有限元大作业matlab课程设计例子.docx
- 文档编号:12540980
- 上传时间:2023-04-20
- 格式:DOCX
- 页数:22
- 大小:93.35KB
有限元大作业matlab课程设计例子.docx
《有限元大作业matlab课程设计例子.docx》由会员分享,可在线阅读,更多相关《有限元大作业matlab课程设计例子.docx(22页珍藏版)》请在冰豆网上搜索。
有限元大作业matlab课程设计例子
有限元大作业matlab---课程设计例子
有限元大作业程序设计
学校:
天津大学
院系:
建筑工程与力学学院
专业:
01级工程力学
姓名:
刘秀
学号:
\\\\\\\\\\\
指导老师:
连续体平面问题的有限元程序分析
[题目]:
如图所示的正方形薄板四周受均匀载荷的作用,该结构在边界上受正向分布压力,
,同时在沿对角线y轴上受一对集中压力,载荷为2KN,若取板厚
,泊松比
。
[分析过程]:
由于连续平板的对称性,只需要取其在第一象限的四分之一部分参加分析,然后人为作出一些辅助线将平板“分割”成若干部分,再为每个部分选择分析单元。
采用将此模型化分为4个全等的直角三角型单元。
利用其对称性,四分之一部分的边界约束,载荷可等效如图所示。
[程序原理及实现]:
用FORTRAN程序的实现。
由节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个一般性的文件DATA.OUT。
模型基本信息由文件为BASIC.IN生成。
该程序的特点如下:
问题类型:
可用于计算弹性力学平面问题和平面应变问题
单元类型:
采用常应变三角形单元
位移模式:
用用线性位移模式
载荷类型:
节点载荷,非节点载荷应先换算为等效节点载荷
材料性质:
弹性体由单一的均匀材料组成
约束方式:
为“0”位移固定约束,为保证无刚体位移,弹性体至少应有对三个自由度的独立约束
方程求解:
针对半带宽刚度方程的Gauss消元法
输入文件:
由手工生成节点信息文件NODE.IN,和单元信息文件ELEMENT.IN
结果文件:
输出一般的结果文件DATA.OUT
程序的原理如框图:
(1)主要变量:
ID:
问题类型码,ID=1时为平面应力问题,ID=2时为平面应变问题
N_NODE:
节点个数
N_LOAD:
节点载荷个数
N_DOF:
自由度,N_DOF=N_NODE*2(平面问题)
N_ELE:
单元个数
N_BAND:
矩阵半带宽
N_BC:
有约束的节点个数
PE:
弹性模量
PR:
泊松比
PT:
厚度
LJK_ELE(I,3):
单元节点编号数组,LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I的三个节点的整体编号
X(N_NODE),Y(N_NODE):
节点坐标数组,X(I),Y(I)分别存放节点I的x,y坐标值
P_LJK(N_BC,3):
节点载荷数组,P_LJK(I,1)表示第I个作用有节点载荷的节点的编号,P_LJK(I,2),P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值
AK(N_DOF,N_BAND):
整体刚度矩阵
AKE(6,6):
单元刚度矩阵
BB(3,6):
位移……应变转换矩阵(三节点单元的几何矩阵)
DD(3,3):
弹性矩阵
SS(3,6);应力矩阵
RESULT_N(N_NOF):
节点载荷数组,存放节点载荷向量,解方程后该矩阵存放节点位移
DISP_E(6):
:
单元的节点位移向量
STS_ELE(N_ELE,3):
单元的应力分量
STS_ND(N_NODE,3):
节点的应力分量
(2)子程序说明:
READ_IN:
读入数据BAND_K:
形成半带宽的整体刚度矩阵
FORM_KE:
计算单元刚度矩阵FORM_P:
计算节点载荷
CAL_AREA:
计算单元面积DO_BC:
处理边界条件
CLA_DD:
计算单元弹性矩阵SOLVE:
计算节点位移
CLA_BB:
计算单元位移……应变关系矩阵
CAL_STS:
计算单元和节点应力
(3)文件管理:
源程序文件:
chengxu.for
程序需读入的数据文件:
BASIC.IN,NODE.IN,ELEMENT.IN(需要手工生成)
程序输出的数据文件:
DATA.OUT
(4)数据文件格式:
需读入的模型基本信息文件BASIC.IN的格式如下表
栏目
格式说明
实际需输入的数据
基本模型数据
第1行,每两个数之间用“,”号隔开
问题类型,单元个数,节点个数,有约束的节点数,有载何的节点数
材料性质
第2行,每两个数之间用“,”号隔开
弹性模量,泊松比,单元厚度
节点约束信息
在材料性质输入行之后另起行,每两个数之间用“,”号隔开
LJK_U(N_BC,3)
位移约束的节点编号,该节点x方向约束代码,该节点y方向代码,
节点荷载信息
在节点约束信息输入行之后另起行,每两个数之间用“,”号隔开
P_IJK(N_LOAD,3)
载荷作用的节点编号,该节点x主向载荷,该节点y方向载荷,……
需读入的节点信息文件NODE.IN的格式如下表
栏目
格式说明
实际需输入的数据
节点信息
每行为一个节点的信息(每行三个数,每两个数之间用空格或“,”分开)
ND_ANSYS(N_NIDE)
节点号,该节点的x坐标,该节点y方向坐标
需读入的单元信息文件ELEMENT.IN的格式如下表
栏目
格式说明
实际需输入的数据
单元信息
每行为一个单元的信息(每行有14个整型数,前4个为单元节点编号,对于3节点编号,第4个节点编号与第3个节点编号相同,后10个数无用,可输入“0”,每两个整型数之间用至少一个空格分开)
NE_ANSYS(N_ELE,14)
单元的节点号1(空格)单元的节点号2(空格)单元的节点号3(空格)单元的节点号4(空格)
0(空格)0(空格)0(空格)0(空格)0(空格)0(空格)0(空格)0(空格)0(空格)0
输出结果文件DATA.OUT格式如下表
栏目
实际输出的数据
节点位移
IRESULT_N(2*I_1)RESULT_N(2*I)
节点号x方向位移y方向位移
单元应力的三个分量
IESTE_ELE(IE,1)STE_ELE(IE,2)STE_ELE(IE,3)
单元号x方向应力y方向应力剪切应力
节点应力的三个分量
ISTS-ND(I,1)STS-ND(I,2)STS-ND(I,3)
节点号x方向应力y方向应力剪切应力
[算例原始数据和程序分析]:
(1)模型基本信息文件BASIC.IN的数据为
1,4,6,5,3
1.,0.,1.
1,1,0,2,1,0,4,1,1,5,0,1,6,0,1
1,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5
(2)手工准备的节点信息文件NODE.IN的数据为
10.02.0
20.01.0
31.01.0
40.0.
51.00.
62.00.
(3)手工准备的单元信息文件ELEMENT.IN的数据为
12330000111101
24550000111102
53220000111103
35660000111104
(4)源程序文件chengxu.for为:
PROGRAMFEM2D
DIMENSIONIJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),
&RESULT_N(500),AK(500,100)
DIMENSIONSTS_ELE(500,3),STS_ND(500,3)
OPEN(4,FILE='BASIC.IN')
OPEN(5,FILE='NODE.IN')
OPEN(6,FILE='ELEMENT.IN')
OPEN(8,FILE='DATA.OUT')
OPEN(9,FILE='FOR_POST.DAT')
READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD
IF(ID.EQ.1)WRITE(8,20)
IF(ID.EQ.2)WRITE(8,25)
20FORMAT(/5X,'=========PLANESTRESSPROBLEM========')
25FORMAT(/5X,'=========PLANESTRAINPROBLEM========')
CALLREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,
&IJK_ELE,X,Y,IJK_U,P_IJK)
CALLBAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,
&IJK_ELE,X,Y,PE,PR,PT,AK)
CALLFORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,
&RESULT_N)
CALLDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)
CALLSOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)
CALLCAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,
&STS_ELE,STS_ND)
ctoputoutadatafile
WRITE(9,70)REAL(N_NODE),REAL(N_ELE)
70FORMAT(2f9.4)
WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),
&STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)
71FORMAT(7F9.4)
WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),
&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),
&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1,N_ELE)
72FORMAT(7f9.4)
c
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(8)
CLOSE(9)
END
c
ctogettheoriginaldatainordertomodeltheproblem
SUBROUTINEREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,
&PT,IJK_ELE,X,Y,IJK_U,P_IJK)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),
&P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)
REALND_ANSYS(N_NODE,3)
READ(4,*)PE,PR,PT
READ(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)
READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)
READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)
READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)
DO10I=1,N_NODE
X(I)=ND_ANSYS(I,2)
Y(I)=ND_ANSYS(I,3)
10CONTINUE
DO11I=1,N_ELE
DO11J=1,3
IJK_ELE(I,J)=NE_ANSYS(I,J)
11CONTINUE
N_BAND=0
DO20IE=1,N_ELE
DO20I=1,3
DO20J=1,3
IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))
IF(N_BAND.LT.IW)N_BAND=IW
20CONTINUE
N_BAND=(N_BAND+1)*2
IF(ID.EQ.1)THEN
ELSE
PE=PE/(1.0-PR*PR)
PR=PR/(1.0-PR)
ENDIF
RETURN
END
c
Ctoformthestiffnessmatrixofelement
SUBROUTINEFORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),
&AKE(6,6),SS(6,6)
CALLCAL_DD(PE,PR,DD)
CALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DO10I=1,3
DO10J=1,6
SS(I,J)=0.0
DO10K=1,3
10SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)
DO20I=1,6
DO20J=1,6
AKE(I,J)=0.0
DO20K=1,3
20AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT
RETURN
END
c
ctoformbandedglobalstiffnessmatrix
SUBROUTINEBAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,
&PR,PT,AK)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)
N_DOF=2*N_NODE
DO40I=1,N_DOF
DO40J=1,N_BAND
40AK(I,J)=0
DO50IE=1,N_ELE
CALLFORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)
DO50I=1,3
DO50II=1,2
IH=2*(I-1)+II
IDH=2*(IJK_ELE(IE,I)-1)+II
DO50J=1,3
DO50JJ=1,2
IL=2*(J-1)+JJ
IZL=2*(IJK_ELE(IE,J)-1)+JJ
IDL=IZL-IDH+1
IF(IDL.LE.0)THEN
ELSE
AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)
ENDIF
50CONTINUE
RETURN
END
c
ctocalculatetheareaofelement
SUBROUTINECAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE)
I=IJK_ELE(IE,1)
J=IJK_ELE(IE,2)
K=IJK_ELE(IE,3)
XIJ=X(J)-X(I)
YIJ=Y(J)-Y(I)
XIK=X(K)-X(I)
YIK=Y(K)-Y(I)
AE=(XIJ*YIK-XIK*YIJ)/2.0
RETURN
END
c
ctocalculatetheelasticmatrixofelement
SUBROUTINECAL_DD(PE,PR,DD)
DIMENSIONDD(3,3)
DO10I=1,3
DO10J=1,3
10DD(I,J)=0.0
DD(1,1)=PE/(1.0-PR*PR)
DD(1,2)=PE*PR/(1.0-PR*PR)
DD(2,1)=DD(1,2)
DD(2,2)=DD(1,1)
DD(3,3)=PE/((1.0+PR)*2.0)
RETURN
END
c
ctocalculatethestrain-displacementmatrixofelement
SUBROUTINECAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)
I=IJK_ELE(IE,1)
J=IJK_ELE(IE,2)
K=IJK_ELE(IE,3)
DO10II=1,3
DO10JJ=1,3
10BB(II,JJ)=0.0
BB(1,1)=Y(J)-Y(K)
BB(1,3)=Y(K)-Y(I)
BB(1,5)=Y(I)-Y(J)
BB(2,2)=X(K)-X(J)
BB(2,4)=X(I)-X(K)
BB(2,6)=X(J)-X(I)
BB(3,1)=BB(2,2)
BB(3,2)=BB(1,1)
BB(3,3)=BB(2,4)
BB(3,4)=BB(1,3)
BB(3,5)=BB(2,6)
BB(3,6)=BB(1,5)
CALLCAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)
DO20I1=1,3
DO20J1=1,6
20BB(I1,J1)=BB(I1,J1)/(2.0*AE)
RETURN
END
c
ctoformthegloballoadmatrix
SUBROUTINEFORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,
&RESULT_N)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),
&RESULT_N(N_DOF)
DO10I=1,N_DOF
10RESULT_N(I)=0.0
DO20I=1,N_LOAD
II=P_IJK(I,1)
RESULT_N(2*II-1)=P_IJK(I,2)
20RESULT_N(2*II)=P_IJK(I,3)
RETURN
END
c
ctodealwithBC(u)(hereonlyforfixeddisplacement)using"1-0"method
SUBROUTINEDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)
DIMENSIONRESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)
DO30I=1,N_BC
IR=IJK_U(I,1)
DO30J=2,3
IF(IJK_U(I,J).EQ.0)THEN
ELSE
II=2*IR+J-3
AK(II,1)=1.0
RESULT_N(II)=0.0
DO10JJ=2,N_BAND
10AK(II,JJ)=0.0
DO20JJ=2,II
20AK(II-JJ+1,JJ)=0.0
ENDIF
30CONTINUE
RETURN
END
c
ctosolvethebandedFEMequationbyGAUSSelimination
SUBROUTINESOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)
DIMENSIONRESULT_N(N_DOF),AK(500,100)
DO20K=1,N_DOF-1
IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1
IF(N_DOF.LE.K+N_BAND-1)IM=N_DOF
DO20I=K+1,IM
L=I-K+1
C=AK(K,L)/AK(K,1)
IW=N_BAND-L+1
DO10J=1,IW
M=J+I-K
10AK(I,J)=AK(I,J)-C*AK(K,M)
20RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)
RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)
DO40I1=1,N_DOF-1
I=N_DOF-I1
IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1
IF(N_BAND.LE.N_DOF-I-1)JQ=N_BAND
DO30J=2,JQ
K=J+I-1
30RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)
40RESULT_N(I)=RESULT_N(I)/AK(I,1)
WRITE(8,50)
50FORMAT(/12X,'*****RESULTSBYFEM2D*****',//8X,
&'--DISPLACEMENTOFNODE--'//5X,'NODENO',8X,'X-DISP',8X,'Y-DISP')
DO60I=1,N_NODE
60WRITE(8,70)I,RESULT_N(2*I-1),RESULT_N(2*I)
70FORMAT(8X,I5,7X,2E15.6)
RETURN
END
c
ccalculatethestresscomponentsofelementandnode
SUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,
&STS_ELE,STS_ND)
DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6),
&SS(3,6),RESULT_N(N_DOF),DISP_E(6)
DIMENSIONSTS_ELE(500,3),STS_ND(500,3)
WRITE(8,10)
10FORMAT(//8X,'--STRESSESOFELEMENT--')
CALLCAL_DD(PE,PR,DD)
DO50IE=1,N_ELE
CALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DO20I=1,3
DO20J=1,6
SS(I,J)=0.0
DO20K=1,3
20SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)
DO30I=1,3
DO30J=1,2
IH=2*(I-1)+J
IW=2*(IJK_ELE(IE,I)-1)+J
30DISP_E(IH)=RESULT_N(IW)
STX=0
STY=0
TXY=0
DO40
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 作业 matlab 课程设计 例子