有限元大作业matlab课程设计例子范文.doc

有限元大作业matlab课程设计例子

28
2020年4月19日
文档仅供参考

有 限 元 大 作 业 程 序 设 计

学校：天津大学
院系：建筑工程与力学学院
专业：01级工程力学
姓名：刘秀
学号：\\\\\\\\\\\
指导老师：

连续体平面问题的有限元程序分析
[题目]：
如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，，同时在沿对角线y轴上受一对集中压力，载荷为2KN，若取板厚，泊松比。

2kN
2kN
1kN/m

[分析过程]：
由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。

1kN/m

[程序原理及实现]：
用FORTRAN程序的实现。由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件为BASIC.IN生成。
该程序的特点如下：
问题类型：可用于计算弹性力学平面问题和平面应变问题
单元类型：采用常应变三角形单元
位移模式：用用线性位移模式
载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷
材料性质：弹性体由单一的均匀材料组成
约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束
方程求解：针对半带宽刚度方程的Gauss消元法
输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN
结果文件：输出一般的结果文件DATA.OUT
程序的原理如框图：
开始

输入数据（子程序READ_IN）
BASIC.IN（基本信息文件）
NODE.IN（节点信息文件）
ELEMENT.IN(单元信息文件)
形成单元刚度矩阵（子程序FORM_KE）
以半带存储方式形成整体刚度矩阵（BAND_K）
形成节点载荷向量（子程序FORM_P）
处理边界条件（子程序DO_BC）
求解方程获得节点位移（子程序SOLVE）
计算单元及节点应力（子程序）
结束
输出方件
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格式如下表
栏目
实际输出的数据
节点位移
I RESULT_N(2*I_ 1) RESULT_N(2*I)
节点号 x方向位移 y方向位移
单元应力的三个分量
IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3)
单元号 x方向应力 y方向应力 剪切应力
节点应力的三个分量
I STS-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的数据为
1 0.0 2.0
2 0.0 1.0
3 1.0 1.0
4 0. 0.
5 1.0 0.
6 2.0 0.
（3）手工准备的单元信息文件ELEMENT.IN的数据为
1 2 3 3 0 0 0 0 1 1 1 1 0 1
2 4 5 5 0 0 0 0 1 1 1 1 0 2
5 3 2 2 0 0 0 0 1 1 1 1 0 3
3 5 6 6 0 0 0 0 1 1 1 1 0 4
（4）源程序文件chengxu.for为：
PROGRAM FEM2D
DIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),
&RESULT_N(500),AK(500,100)
DIMENSION STS_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)
20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')
25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')
CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,
& IJK_ELE,X,Y,IJK_U,P_IJK)
CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,
& IJK_ELE,X,Y,PE,PR,PT,AK)
CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,
& RESULT_N)
CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)
CALL SOLVE(N_NODE,N_DOF,NBAND,AK,RESULT_N)
CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,
& STS_ELE,STS_ND)
c to putout a data file
WRITE(9,70)REAL(N_NODE),REAL(N_ELE)
70 FORMAT(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_NDI,3),I=1,N_NODE)
71 FORMAT(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)
72 FORMAT(7f9.4)
c
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(8)
CLOSE(9)
END
c
c to get the original data in order to model the problem

SUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,
&PT,IJK_ELE,X,Y,IJK_U,P_IJK)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),
& P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)
REAL ND_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)
DO 10 I=1,N_NODE
X(I)=ND_ANSYS(I,2)
Y(I)=ND_ANSYS(I,3)
10 CONTINUE
DO 11 I=1,N_ELE
DO 11 J=1,3
IJK_ELE(I,J)=NE_ANSYS(I,J)
11 CONTINUE
N_BAND=0
DO 20 IE=1,N_ELE
DO 20 I=1,3
DO 20 J=1,3
IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))
IF(N_BAND.LT.IW)N_BAND=IW
20 CONTINUE
N_BAND=(N_BAND+1)*2
IF(ID.EQ.1) THEN
ELSE
PE=PE/(1.0-PR*PR)
PR=PR/(1.0-PR)
END IF
RETURN
END
c
C to form the stiffness matrix of element
SUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),
& AKE(6,6), SS(6,6)
CALL CAL_DD(PE,PR,DD)
CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DO 10 I=1,3
DO 10 J=1,6
SS(I,J)=0.0
DO 10 K=1,3
10 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)
DO 20 I=1,6
DO 20 J=1,6
AKE(I,J)=0.0
DO 20 K=1,3
20 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT
RETURN
END
c
c to form banded global stiffness matrix
SUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,
& PR,PT,AK)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)
N_DOF=2*N_NODE
DO 40 I=1,N_DOF
DO 40 J=1,N_BAND
40 AK(I,J)=0
DO 50 IE=1,N_ELE
CALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PRPT,AKE)
DO 50 I=1,3
DO 50 II=1,2
IH=2*(I-1)+II
IDH=2*(IJK_ELE(IE,I)-1)+II
DO 50 J=1,3
DO 50 JJ=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)
END IF
50 CONTINUE
RETURN
END
c
c to calculate the area of element
SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)
DIMENSION IJK_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
c to calculate the elastic matrix of element
SUBROUTINE CAL_DD(PE,PR,DD)
DIMENSION DD(3,3)
DO 10 I=1,3
DO 10 J=1,3
10 DD(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
c to calculate the strain-displacement matrix of element
SUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DIMENSION IJK_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)
DO 10 II=1,3
DO 10 JJ=1,3
10 BB(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)
CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)
DO 20 I1=1,3
DO 20 J1=1,6
20 BB(I1,J1)=BB(I1,J1)/(2.0*AE)
RETURN
END
c
c to form the global load matrix
SUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,
& RESULT_N)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),
& RESULT_N(N_DOF)
DO 10 I=1,N_DOF
10 RESULT_N(I)=0.0
DO 20 I=1,N_LOAD
II=P_IJK(I,1)
RESULT_N(2*II-1)=P_IJK(I,2)
20 RESULT_N(2*II)=P_IJK(I,3)
RETURN
END
c
c to deal with BC(u) (here only for fixed displacement) using "1-0" method
SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)
DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)
DO 30 I=1,N_BC
IR=IJK_U(I,1)
DO 30 J=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
DO 10 JJ=2,N_BAND
10 AK(II,JJ)=0.0
DO 20 JJ=2,II
20 AK(II-JJ+1,JJ)=0.0
END IF
30 CONTINUE
RETURN
END
c
c to solve the banded FEM equation by GAUSS elimination
SUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)
DIMENSION RESULT_N(N_DOF),AK(500,100)
DO 20 K=1,N_DOF-1
IF(N_DOF.GT.K+N_BAND-1)IM=K+NBAND-1
IF(N_DOF.LE.K+N_BAND-1)IM=N_DOF
DO 20 I=K+1,IM
L=I-K+1
C=AK(K,L)/AK(K,1)
IW=N_BAND-L+1
DO 10 J=1,IW
M=J+I-K
10 AK(I,J)=AK(I,J)-C*AK(K,M)
20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)
RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)
DO 40 I1=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
DO 30 J=2,JQ
K=J+I-1
30 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)
40 RESULT_N(I)=RESULT_N(I)/AK(I,1)
WRITE(8,50)
50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,
&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP')
DO 60 I=1,N_NODE
60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)
70 FORMAT(8X,I5,7X,2E15.6)
RETURN
END
c
c calculate the stress components of element and node
SUBROUTINE CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,YRESULT_N,
&STS_ELE,STS_ND)
DIMENSION IJK_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)
DIMENSION STS_ELE(500,3),STS_ND(500,3)
WRITE(8,10)
10 FORMAT(//8X,'--STRESSES OF ELEMENT--')
CALL CAL_DD(PE,PR,DD)
DO 50 IE=1,N_ELE
CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DO 20 I=1,3
DO 20 J=1,6
SS(I,J)=0.0
DO 20 K=1,3
20 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)
DO 30 I=1,3
DO 30 J=1,2
IH=2*(I-1)+J
IW=2*(IJK_ELE(IE,I)-1)+J
30 DISP_E(IH)=RESULT_N(IW)
STX=0
STY=0
TXY=0
DO 40 J=1,6
STX=STX+SS(1,J)*DISP_E(J)
STY=STY+SS(2,J)*DISP_E(J)
40 TXY=TXY+SS(3,J)*DISP_E(J)
STS_ELE(IE,1)=STX
STS_ELE(IE,2)=STY
STS_ELE(IE,3)=TXY
50 WRITE(8,60)IE,STX,STY,TXY
60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',
&E12.6,2X,'TXY=',E12.6)
c the following part is to calculate stress components of node
WRITE(8,55)
55 FORMAT(//8X,'--STRESSES OF NODE--')
DO 90 I=1,N_NODE
A=0.
B=0.
C=0.
II=0
DO 70 K=1,N_ELE
DO 70 J=1,3
IF(IJK_ELE(K,J).EQ.I) THEN
II=II+1
A=A+STS_ELE(K,1)
B=B+STS_ELE(K,2)
C=C+STS_ELE(K,3)
END IF
70 CONTINUE
STS_ND(I,1)=A/II
STS_ND(I,2)=B/II
STS_ND(I,3)=C/II
WRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)
75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',
&E12.6,2X,'TXY=',E12.6)
90 CONTINUE
RETURN
END
c FEM2D programm end
[算例结果]：
chengxu.for所输出的数据文件DATA.OUT数据内容如下：

=========PLANE STRESS PROBLEM========

* * * * * RESULTS BY FEM2D * * * * *

--DISPLACEMENT OF NODE--

NODE NO X-DISP Y-DISP
1 .000000E00 -.525275E+01
2 .000000E+00 -.225275E+01
3 -.108791E+01 -.137363E+01
4 .000000E+00 .000000E+00
5 -.824176E+00 .000000E+00
6 -.182418E+01 .000000E+00

--STRESSES OF ELEMENT--
ELEMENT NO.= 1
STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00
ELEMENT NO.= 2
STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00
ELEMENT NO.= 3
STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00
ELEMENT NO.= 4
STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00

--STRESSES OF NODE--
NODE NO.= 1
STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00
NODE NO.= 2
STX=-.100000E+01 STY=-.220879E+01 TXY .249084E+00
NODE NO.= 3
STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00
NODE NO.= 4
STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00
NODE NO.= 5
STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01
NODE NO. 6
STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00

[结论与体会]:
经过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。经过编程及对chengxu.for的调试进一步巩固了FORTRAN语言知识，也锻炼了自我，使自己在学习工作中更具耐心、细心。
[参考数目]：
《FORTRAN语言》清华大学出版社 谭浩强，田淑清
《弹性力学》 高等教育出版社 徐芝纶
《有限元分析及其应用》清华大学 曾攀