平面三角形3节点有限元程序Word格式.docx
《平面三角形3节点有限元程序Word格式.docx》由会员分享，可在线阅读，更多相关《平面三角形3节点有限元程序Word格式.docx（23页珍藏版）》请在冰豆网上搜索。
TREAT——由于有非零位移，对K和F进展处理；
DEP——整体劲度矩阵的分解运算；
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
Ix，Iy——分别为该结点的约束情况，如果方向受约束时填0，如果自由如此填1。
⑹荷载信息，共NL条，每条依次输入y
IP，Fx，Fy
Fx，Fy——分别为该结点的x，y方向的荷载分量〔kN〕。
结点号存放在数组NF〔NL〕中，结点荷载分量存放在数组FV〔2，NL〕中。
⑺假如ND＞0，输入非零位移信息，共ND条，每条依次输入
IP，ux，uy
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-P.Y-P.
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
5-1.000-1.0000.000-1.000-1.00090.000
6-1.000-1.0000.000-1.000-1.00090.000
NODALREACTIONS
NODEX-PY-P
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
MON/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）
CALLDEP（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）
400FORMAT（/2X,'
NP='
I3,2X,'
NE='
NM='
*,I3,2X,'
NR='
NI='
I3,2X,'
NL='
I3,2X,
*'
NG='
ND='
NC='
I3）
500FORMAT（1X,'
TOTALDEGREESOFFREEDOMN='
*I4,1X,'
TOTAL-STORAGE'
'
NH='
I5,1X,
MAX-SEMI-BANDWIDTHMX='
550FORMAT（/20X,'
TOTALSTORAGEIS
*GREATERTHAN50000'
600FORMAT（30X,'
NODALFORCES'
/8X,'
NODE'
*11X,'
X-P.'
14X,'
Y-P.'
650FORMAT（/30X,'
NODALDISPLACEMENTS'
/8X,
13X,'
12X,'
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）
MON/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
升级会员 