平面三角形3节点有限元程序.docx
《平面三角形3节点有限元程序.docx》由会员分享，可在线阅读，更多相关《平面三角形3节点有限元程序.docx（27页珍藏版）》请在冰豆网上搜索。
平面三角形3节点有限元程序
平面三角形３结点有限元程序
1、程序名：
FEMT3.FOR，FEMT3.EXE
2、程序功能
本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种
荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：
结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结
点的支座反力。
程序采用VisualFortran5.0编制而成，输入数据全部采用自由格式。
3、程序流程及框图
图１-1程序流程图
图１-2程序框图
其中，各子程序的功能如下：
INPUT——输入结点坐标、单元信息和材料参数；
MR——形成结点自由度序号矩阵；
FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；
DIV——取出单元的3个结点号码和该单元的材料号并计算单元的bi，ci等；
MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；
LOAD——形成整体结点荷载列阵F；
OUTPUT——输出结点位移或结点荷载；
TREAT——由于有非零已知位移，对K和F进行处理；
DECOMP——整体劲度矩阵的分解运算；
FOBA——前代、回代求出未知结点位移δ；
ERFAC——计算约束结点的支座反力；
KRS——计算单元劲度矩阵中的子块Krs。
4、输入数据及变量说明
当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个
文件的名字由少于8个字符或数字组成。
数据文件包括如下内容：
⑴ 总控信息，共一条，9个数据
NP，NE，NM，NR，NI，NL，NG，ND，NC
NP——结点总数；
NE——单元总数；
NM——材料类型总数；
NR——约束结点总数；
NI——问题类型标识，0为平面应力问题，1为平面应变问题；
NL——受荷载作用的结点的数目；
NG——考虑自重作用为1，不计自重为0；
ND——非零已知位移结点的数目；
NC——要计算支座约束反力的结点数目。
⑵ 材料信息，共NM条，每条依次输入
EO，VO，W，t
EO——弹性模量（kN/m2）；
VO——泊松比；
W——材料容重（kN/m3）；
t——单元厚度（m）。
这些信息都存放在数组AE（4，NM）中。
⑶ 坐标信息，共NP条，每条依次输入
IP，X，Y
IP——结点号；
X，Y——分别为结点的x坐标和y坐标。
坐标信息存放在数组X（2，NP）中。
⑷ 单元信息，共NE条，每条依次输入
JE，L，Io，Jo，Mo
JE——单元号；
L——为该单元的材料类型号。
Io，Jo，Mo——分别为该单元i，j，m的整体编码。
单元信息存放在数组MEO（4,NE）中。
⑸ 约束信息共NR条，每条依次输入一个数
××
IP
IP——结点号；
Ix，Iy——分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。
⑹ 荷载信息，共NL条，每条依次输入y
IP，Fx，Fy
IP——结点号；
Fx，Fy——分别为该结点的x，y方向的荷载分量（kN）。
结点号存放在数组NF（NL）中，结点荷载分量存放在数组FV（2，NL）中。
⑺ 若ND＞0，输入非零已知位移信息，共ND条，每条依次输入
IP，ux，uy
IP——结点号；
ux，uy——分别为该结点x，y方向已知位移分量（m），若其中某方向为自由，则其相应分量为0。
结点号存放在数NDI（ND）中，已知位移分量存放在数组DV（2，ND）中。
⑻ 支座反力信息，共NC条，每条依次输入
IP，M1，M2，M3，M4
IP——支座结点号；
M1，M2，M3，M4——为与该支座结点相关的单元号，若不足4个，则用0补充。
支座结点号存放在数组NCI（NC）中，相关单元号存放在数组NCE（4，NC）中。
以上数据须按如上顺序存放在数据文件中。
除此之外，程序中还用到其他一些主要变量和数组，说明如下：
N——结构自由度总数；
NH——按一维存贮的整体劲度矩阵的总容量；
NX——最大半带宽；
SK（10000）——维存贮的劲度矩阵；
R（1000）——开始存放等效结点荷载，求解方程以后，用来存放结点位移；
B（6）——存放单元应力σx，σy，τxy，σ1，σ2，α；
MA（1000）——主元素序号指标矩阵；
JR（2，500）——结点自由度序号矩阵；
ME（3）——存放单元结点i，j，m的整体编码；
NN（6）——单元结点自由度序号；
BI（3），CI（3）——单元劲度矩阵计算公式中的bi，bj，bm和ci，cj，cm；
S——三角形单元的面积；
H11，H12，H21，H22——单元劲度矩阵中子块Krs的4个元素。
5、算例
一个正方形弹性体，厚度为1m，四边受单位均布法向力作用，由于对称性，取其1/4进行计算，其有限元网格如图2-3所示，设
，
，不考虑自重。
该问题的精确解应力为
＝1
，
＝1
，
＝0。
图１-3有限元网格
（1）输入文件数据
641503005
2000.00.00.01.0
10.02.0
20.01.0
31.01.0
40.00.0
51.00.0
62.00.0
11312
21245
31325
41563
101
201
400
510
610
1-0.5-0.5
3-1.0-1.0
6-0.5-0.5
11000
21230
42000
52340
64000
（2）输出文件结果
NODALDISPLACEMENTS
NODEX-COMP.Y-COMP.
10.00000E+00-0.10000E-02
20.00000E+00-0.50000E-03
3-0.50000E-03-0.50000E-03
40.00000E+000.00000E+00
5-0.50000E-030.00000E+00
6-0.10000E-020.00000E+00
ELEMENTSTRESSES
ELEMENTX-STRESSY-STRESSXY-STRESSMAX-STRESSMIN-STRESSANGLE
1-1.000-1.0000.000-1.000-1.00090.000
2-1.000-1.0000.000-1.000-1.00090.000
3-1.000-1.0000.000-1.000-1.00090.000
4-1.000-1.0000.000-1.000-1.00090.000
NODESTRESSES
NODEX-STRESSY-STRESSXY-STRESSMAX-STRESSMIN-STRESSANGLE
1-1.000-1.0000.000-1.000-1.00090.000
2-1.000-1.0000.000-1.000-1.00090.000
3-1.000-1.0000.000-1.000-1.00090.000
4-1.000-1.0000.000-1.000-1.00090.000
5-1.000-1.0000.000-1.000-1.00090.000
6-1.000-1.0000.000-1.000-1.00090.000
NODALREACTIONS
NODEX-COMPY-COMP
10.0000.000
21.0000.000
40.5000.500
50.0001.000
60.0000.000
6、源程序
CFINITEELEMENTPROGRAMFORTWODIMENSIONAL
CTRIANGLEELEMENT
C
DIMENSIONK（800000）,COOR（2,3000）,AE（4,11）,
*MEL（5,2000）,MA（6000）
CHARACTER*32dat
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
WRITE（*,300）
300FORMAT（///''
*':
****'/'+PLEASEINPUTFILENAMEOFDATA'）
READ（*,*）data
OPEN（4,FILE=data,STATUS='OLD'）
OPEN（7,FILE='OUT',STATUS='UNKNOWN'）
READ（4,*）NP,NE,NM,NR,NI,NL,NG,ND,NC
CWRITE（*,400）NP,NE,NM,NR,NI,NL,NG,ND,NC
CWRITE（7,400）NP,NE,NM,NR,NI,NL,NG,ND,NC
CALLINPUT（JR,COOR,MEL,AE）
CALLCBAND（MA,JR,MEL）
CALLSK0（SK,R,COOR,MEL,MA,JR,AE）
CALLLOAD（COOR,MEL,R,JR,AE）
IF（ND.GT.0）CALLTREAT（SK,MA,JR,R）
CALLDECOMP（SK,MA）
CALLFOBA（SK,MA,R）
WRITE（*,650）
WRITE（7,650）
CALLOUTPUT（JR,R）
WRITE（*,700）
WRITE（7,700）
CALLCES（COOR,MEL,JR,R,AE）
IF（NC.GT.0）callERFAC（COOR,MEL,JR,R,AE）
400FORMAT（/2X,'NP=',I3,2X,'NE=',I3,2X,'NM='
*,I3,2X,'NR=',I3,2X,'NI='I3,2X,'NL=',I3,2X,
*'NG=',I3,2X,'ND=',I3,2X,'NC=',I3）
500FORMAT（1X,'TOTALDEGREESOFFREEDOMN=',
*I4,1X,'TOTAL-STORAGE','NH=',I5,1X,
*'MAX-SEMI-BANDWIDTHMX=',I3）
550FORMAT（/20X,'TOTALSTORAGEIS
*GREATERTHAN50000'）
600FORMAT（30X,'NODALFORCES'/8X,'NODE',
*11X,'X-COMP.',14X,'Y-COMP.'）
650FORMAT（/30X,'NODALDISPLACEMENTS'/8X,
*'NODE',13X,'X-COMP.',12X,'Y-COMP.'）
700FORMAT（/30X,'ELEMENTSTRESSES'/5X,
*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',
*2X,'XY-STRESS',1X,'MAX-STRESS',1X,
*'MIN-STRESS',6X,'ANGLE'/）
STOP
END
C*********************************************
SUBROUTINEKRS（BR,BS,CR,CS）
COMMON/CB/EO,VO,W,T,A,H11,H12,H21,H22
*,ME（3）,BI（3）,CI（3）
ET=EO*T/（1.0-VO*VO）/A/4.0
V=（1.0-VO）/2.0
H11=ET*（BR*BS+V*CR*CS）
H12=ET*（VO*BR*CS+V*BS*CR）
H21=ET*（VO*CR*BS+V*BR*CS）
H22=ET*（CR*CS+V*BR*BS）
RETURN
END
C*********************************************
SUBROUTINEINPUT（JR,COOR,MEL,AE）
DIMENSIONJR（2,*）,COOR（2,*）,AE（4,*）,MEL（3,*）
COMMON/CA/NP,NE,NM,NR
COMMON/CC/N,MX,NH
DO70I=1,NP
READ（4,*）IP,X,Y
COOR（1,IP）=X
COOR（2,IP）=Y
70CONTINUE
DO11J=1,NE
READ（4,*）NEE,NME,（MEL（I,NEE）,I=1,3）
MEL（3,NEE）=NME
11CONTINUE
DO10I=1,NP
DO10J=1,2
10JR（J,I）=1
DO20I=1,NR
READ（4,*）IP,IX,IY
JR（1,IP）=IX
JR（2,IP）=IY
20CONTINUE
N=0
DO30I=1,NP
DO30J=1,2
IF（JR（J,I））30,30,25
25N=N+1
JR（J,I）=N
30CONTINUE
DO55J=1,NM
READ（4,*）jj,（AE（I,jj）,I=1,4）
CWRITE（*,910）jj,（AE（I,jj）,I=1,4）
IF（NI.eq.1）then
AE（1,jj）=AE（1,jj）/（1.0-AE（2,jj）*AE（2,jj））
AE（2,jj）=AE（2,jj）/（1.0-AE（2,jj））
endif
55CONTINUE
910FORMAT（/20X,'MATERIALPROPERTIES'/（3X,I5,4（1x,E8.3）））
RETURN
END
C**********************************************
SUBROUTINECBAND（MA,JR,MEL）
DIMENSIONMA（*）,JR（2,*）,MEL（3,*）,NN（6）
COMMON/CA/NP,NE,NM,NR
COMMON/CC/N,MX,NH
DO65I=1,N
65MA（I）=0
DO90IE=1,NE
DO75K=1,3
IEK=MEL（K,IE）
DO95M=1,2
JJ=2*（K-1）+M
NN（JJ）=JR（M,IEK）
95CONTINUE
75CONTINUE
L=N
DO80I=1,6
NNI=NN（I）
IF（NNI.EQ.0）GOTO80
IF（NNI.LT.L）L=NNI
80CONTINUE
DO85M=1,6
JP=NN（M）
IF（JP.EQ.0）GOTO85
JPL=JP-L+1
IF（JPL.GT.MA（JP））MA（JP）=JPL
85CONTINUE
90CONTINUE
MX=0
MA
（1）=1
DO10I=2,N
IF（MA（I）.GT.MX）MX=MA（I）
MA（I）=MA（I）+MA（I-1）
10CONTINUE
NH=MA（N）
WRITE（*,500）N,MX,NH
WRITE（7,500）N,MX,NH
500FORMAT（/5X,'FREEDOMN='
*,I5,3X,'SEMI-BANDWI.MX=',I5,3X,
*'STORAGENH=',I7）
RETURN
END
C***********************************************
SUBROUTINESK0（SK,R,COOR,MEL,MA,JR,AE）
DIMENSIONAE（4,*）,COOR（2,*）,MEL（3,*）,JR（2,*）,R（*）,
*MA（*）,SK（*）,SKE（6,6）,NN（6）
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CB/EO,VO,W,T,A,H11,H12,H21,H22,
*ME（3）,BI（3）,CI（3）
COMMON/CC/N,NH
DO10I=1,NH
10SK（I）=0.0
DO70IE=1,NE
CALLDIV（IE,COOR,MEL,AE）
DO30I=1,3
DO30J=1,3
CALLKRS（BI（I）,BI（J）,CI（I）,CI（J））
SKE（2*I-1,2*J-1）=H11
SKE（2*I-1,2*J）=H12
SKE（2*I,2*J-1）=H21
SKE（2*I,2*J）=H22
30CONTINUE
DO40I=1,3
J2=ME（I）
DO40J=1,2
J3=2*（I-1）+J
NN（J3）=JR（J,J2）
40CONTINUE
DO60I=1,6
DO60J=1,6
IF（NN（J）.EQ.0.OR.NN（I）.LT.NN（J））GOTO60
JJ=NN（I）
JK=NN（J）
JL=MA（JJ）
JM=JJ-JK
JN=JL-JM
SK（JN）=SK（JN）+SKE（I,J）
60CONTINUE
70CONTINUE
cWRITE（0,500）（SK（I）,I=1,20）
500FORMAT（/10X,'SK='/（6F12.5））
RETURN
END
C**********************************************
SUBROUTINELOAD（COOR,MEL,R,JR,AE）
DIMENSIONAE（4,*）,COOR（2,*）,MEL（3,*）,R（*）,JR（2,*）,
*NF（50）,FV（2,50）
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CB/EO,VO,W,T,A,A1（4）,ME（3）,BB（6）
COMMON/CC/N,NH
DO10I=1,N
10R（I）=0.0
IF（NG）70,70,30
30DO60IE=1,NE
CALLDIV（IE,COOR,MEL,AE）
DO50I=1,3
J2=ME（I）
J3=JR（2,J2）
IF（J3）50,50,40
40R（J3）=R（J3）-T*W*A/3.0
50CONTINUE
60CONTINUE
70IF（NL）110,110,80
80READ（4,*）（NF（I）,I=1,NL）
READ（4,*）（（FV（I,J）,I=1,2）,J=1,NL）
CWRITE（*,500）（NF（I）,I=1,NL）
CWRITE（7,500）（NF（I）,I=1,NL）
CWRITE（*,600）（（FV（I,J）,I=1,2）,J=1,NL）
CWRITE（7,600）（（FV（I,J）,I=1,2）,J=1,NL）
DO100I=1,NL
JJ=NF（I）
J=JR（1,JJ）
M=JR（2,JJ）
IF（J.GT.0）R（J）=R（J）+FV（1,I）
IF（M.GT.0）R（M）=R（M）+FV（2,I）
100CONTINUE
110RETURN
500FORMAT（/20X,'NODESOFAPPLIEDLOAD***NF='
*/（1X,10I8））
600FORMAT（/30X,'LUMPED-LOADS***FV='
*/（5X,5F15.3））
END
C*********************************************
SUBROUTINETREAT（SK,MA,JR,R）
DIMENSIONSK（*）,MA（*）,NDI（75）,DV（2,75）,JR（2,*）,R（*）
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CC/N,NH
READ（4,*）（NDI（J）,J=1,ND）
READ（4,*）（（DV（I,J）,I=1,2）,J=1,ND）
CWRITE（*,500）（NDI（J）,J=1,ND）
CWRITE（7,500）（NDI（J）,J=1,ND）
CWRITE（*,550）（（DV（I,J）,I=1,2）,J=1,ND）
CWRITE（7,550）（（DV（I,J）,I=1,2）,J=1,ND）
DO20I=1,ND
DO20J=1,2
IF（DV（J,I））10,20,10
10JJ=NDI（I）
L=JR（J,JJ）
JN=MA（L）
SK（JN）=1.0E30
R（L）=DV（J,I）*1.0E30
20CONTINUE
500FORMAT（/25X,'NODENO.**NDI='/（1X,10I8））
550FORMAT（/25X,'DISPLACEMENT-VALUES**DV='/
*（10X,6F10.6））
RETURN
END
C*********************************************
SUBROUTINEDECOMP（SK,MA）
DIMENSIONSK（*）,MA（*）
COMMON/CC/N,NH
DO50I=2,N
L=I-MA（I）+MA（I-1）+1
K=I-1
L1=L+1
IF（L1.GT.K）GOTO30
DO20J=L1,K
IJ=MA（I）-I+J
M=J-MA（J）+MA（J-1）+1
IF（L.GT.M）M=L
MP=J-1
IF（M.GT.MP）GOTO20
DO10LP=M,MP
IP=MA（I）-I+LP
JP=MA（J）-J+LP
SK（IJ）=SK（IJ）-SK（IP）*SK（JP）
10CONTINUE
20CONTINUE
30IF（L.GT.K）GOTO50
DO40LP=L,K
升级会员 