平面三角形3节点有限元程序Word格式.docx

平面三角形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