有限元大作业matlab课程设计例子.docx
《有限元大作业matlab课程设计例子.docx》由会员分享，可在线阅读，更多相关《有限元大作业matlab课程设计例子.docx（19页珍藏版）》请在冰豆网上搜索。
有限元大作业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
DO40J=1,6
S
升级会员 