连续体平面问题的有限元分析.docx
《连续体平面问题的有限元分析.docx》由会员分享，可在线阅读，更多相关《连续体平面问题的有限元分析.docx（19页珍藏版）》请在冰豆网上搜索。
连续体平面问题的有限元分析
连续体平面问题的有限元分析
【题目】：
正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，P=1kN/m，同时在沿对角线y轴上受一对集中压力，载荷为2kN，板厚t=1，泊松比ν=0，见下图：
Y
2kN
1kN/m
x
o
2kN
摘要：
有限单元法作为一门课程在现实工程中的应用已经十分的广泛，利本文用计算机，结合FORTRAN语言和有限单元法的理论课知识，对平面应力问题的薄板进行内力值的计算和相应的结构分析，并可以将由计算机计算的实际结果与由用弹性力学理论知识计算出来的解答相比较。
关键词：
连续体；平面应力；应力；节点；单元
一、理论依据与分析：
此问题，为弹性力学里的平面应力问题，在板的内部，到处都有【1】：
σz=0,τyz=0,τxz=0;σx=f1（x,y）,σy=f2（x,y）,τxy=f3（x,y），应力具有这种性质的问题，称为平面应力问题。
弹性薄板在工程中应用很广泛，对于一些简单的情况，如等厚、单跨、无大孔口，外形规则（如矩形，圆形等）的薄板，已有一些解答和表格可资利用盘【2】。
由于连续平板的连续性，仅需要取其在第一象限的四分之一部分研究计算，然后做出一些辅助线将平板分成若干部分，在为每个部分选择分子单元，采用此模型化为4个全等的三角形单元，利用其对称性，四分之一的边界约束，荷载可等效如图所示。
1kN/m
二、程序原理及实现：
用FORTRAN程序的实现。
有节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个一般性的文件DATA.OUT。
模型基本信息由文件BASIC.IN生成。
该程序的特点如下：
问题类型：
可用于计算弹性力学平面应力问题和平面应变问题。
单元类型：
采用常应变三角形单元。
位移模式：
用线性位移模式。
载荷类型：
节点载荷，非节点载荷应先转换为等效节点载荷。
材料性质：
弹性体由单一的均匀材料组成。
约束方式：
为'“0”位移固定约束，为保证无刚体位移，弹性体至少应有三个自由度的独立约束。
方程求解：
针对半带宽刚度方程的GUASS消去法。
输出文件：
由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT,IN。
结果文件：
输出一般的结果文件DATA.OUT。
三、程序的原理如框图：
开始
输入数据（子程序READ_IN）
BASIC.IN（基本信息文件）
NODE.IN（节点信息文件）
ELEMEENT.IN（单元信息文件）
↓
形成单元刚度矩阵（子程序FORM_KE）
↓
以半带存储方式形成整体刚度矩阵
（BANDK）
↓
形成节点载荷向量（子程序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坐标值。
LJK_U（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_FE:
计算单元刚度矩阵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的格式如下表
栏目
格式说明
实际需输入的数据
节点信息
每行为一个节点的信息（每行三个数，每两个数之间用空格或“，”分开
LJK_U（N_BC,3）
节点号，该节点的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,100）,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,'=========PLANESTRAIN
PROBLEM========='）
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）
CALLFROM_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（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（5,*）（（ND_ANSYS（I,J）,J=1,3）,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（IJKK_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）=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
Ctoformbandedglobalstiffnessmatrix
SUBROUTINEBAND_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
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
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）
YJJ=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（IIJJ）=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,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）
DIMENTIONRESULT_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,'*****RESULTBYFEM2D*****',//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）
RESULT
END
C
Ccalculatethestresscomponentofelementandnode
SUBROUTINE
CAL_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,'--STRESSOFELEMENT--'）
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,L）
DO30I=1,3
DO30I=1,3
DO30
升级会员 