有限元程序.docx

有限元程序.docx

《有限元程序.docx》由会员分享，可在线阅读，更多相关《有限元程序.docx（55页珍藏版）》请在冰豆网上搜索。

有限元程序

有限单元法程序设计

学生姓名：

汪菊

学号：

094812284

专业：

结构工程

指导老师：

周文伟

目录

程序1：

平面三节点有限元程序………………………………………..1

程序2：

四节点矩形薄板单元程序……………………………………..8

程序3：

平面刚架静力分析程序（PF.FOR）……………………………20

程序一：

平面三节点有限元程序

如图1所示的悬臂梁，受均布荷载q=1N/mm2作用。

E=2.1×105N/mm2，

，厚度h=10mm，现用有限元法分析其位移及应力。

图1

梁可视为平面应力状态，按图1尺寸划分为均匀的三角形网格，共80个单元，55个节点，坐标轴及单元与节点的编号如图。

将均不荷载分配到相应的节点上，把有约束的节点51、52、53、54、55视作固定铰链，建立如图1所示的离散化计算模型。

计算源程序如下：

PROGRAMMAIN

DIMENSIONSK（300,30）,EK（12,12）,Q（300）,MC（55）,XY（3,100）,XYE（3,4）,&

QE（12）,NX（4,100）

OPEN（7,FILE='INP.DAT'）

REWIND7

READ（7,10）NF,NE,NN,MB,ND,LE,LS

READ（7,12）E,UM,T

10FORMAT（7I5）

12FORMAT（3F15.2）

WRITE（*,600）NF,NE,NN,MB,ND,LE,LS,E,UM,T

ME=NE*NF

MS=NN*NF

CALLINPUT（XY,Q,NX,MC,LS,NN,MS,NE,LE,ND）

WRITE（*,102）（（XY（I,J）,I=1,LS）,J=1,NN）

102FORMAT（10X,'XY'/,（2X,6F12.3））

WRITE（*,101）（Q（I）,I=1,MS）

101FORMAT（10X,'Q'/,（2X,6F12.3））

WRITE（*,500）（（NX（I,J）,I=1,NE）,J=1,LE）

WRITE（*,400）（MC（I）,I=1,ND）

500FORMAT（10X,'NX'/,（2X,12I6））

600FORMAT（10X,'NFNENNMBNDLELSEUMT'/7（2X,I4）,3（2X,F8.4））

400FORMAT（10X,'MC'/,（2X,10I6））

CALLSTIFS（SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,&

UM,T）

CALLSOLVE（SK,Q,MS,MB）

OPEN（9,FILE='OUT.DAT'）

REWIND9

WRITE（9,200）

WRITE（9,250）（Q（I）,I=1,MS）

200FORMAT（5X,'DISPLACEMENT'）

250FORMAT（2X,6E14.5）

CALLSTRES（Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T）

STOP1000

END

SUBROUTINEINPUT（XY,Q,NX,MC,LS,NN,MS,NE,LE,ND）

DIMENSIONXY（LS,NN）,Q（MS）,NX（NE,LE）,MC（ND）

READ（7,*）XY

READ（7,*）Q

READ（7,*）NX

READ（7,*）MC

CLOSE（7）

10FORMAT（6F11.2）

20FORMAT（12I5）

RETURN

END

SUBROUTINESTIFS（SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,UM,T）

DIMENSIONSK（MS,MB）,EK（ME,ME）,Q（MS）,NX（NE,LE）,MC（ND）,XY（LS,NN）,XYE（LS,NE）

DO35I=1,MS

DO35J=1,MB

35SK（I,J）=0

DO200L=1,LE

DO40J=1,NE

LJ=NX（J,L）

DO40I=1,LS

40XYE（I,J）=XY（I,LJ）

DO50I=1,ME

DO50J=1,ME

50EK（I,J）=0.0

CALLSTIFE（EK,XYE,ME,NE,NF,LS,E,UM,T）

IF（L.EQ.1）WRITE（*,70）EK

70FORMAT（10X,'EK'/,（6E14.5））

DO200I=1,NE

DO200II=1,NF

M=NF*（I-1）+II

M1=NF*（NX（I,L）-1）+II

DO200J=1,NE

DO200JJ=1,NF

N=NF*（J-1）+JJ

N1=NF*（NX（J,L）-1）+JJ

MN=N1-M1+1

IF（MN）200,200,150

150SK（M1,MN）=SK（M1,MN）+EK（M,N）

200CONTINUE

DO220I=1,ND

M=MC（I）

220SK（M,1）=SK（M,1）*1E8

RETURN

END

SUBROUTINESOLVE（SK,Q,MS,MB）

DIMENSIONSK（MS,MB）,Q（MS）

K1=MS-1

DO125K=1,K1

IF（K+MB-1-MS）105,106,106

105N=K+MB-1

GOTO110

106N=MS

110I1=K+1

DO125I=I1,N

L=I-K+1

C=SK（K,L）/SK（K,1）

J1=MB-L+1

DO122J=1,J1

122SK（I,J）=SK（I,J）-C*SK（K,M）

125Q（I）=Q（I）-C*Q（K）

Q（MS）=Q（MS）/SK（MS,1）

M=MS-1

DO145I1=1,M

I=MS-I1

IF（MS-I+1-MB）135,136,136

135N=MS-I+1

GOTO140

136N=MB

140DO142J=2,N

L=J+I-1

142Q（I）=Q（I）-SK（I,J）*Q（L）

145Q（I）=Q（I）/SK（I,1）

WRITE（*,147）

147FORMAT（5X,'DISPLACEMENT'）

WRITE（*,150）（Q（I）,I=1,MS）

150FORMAT（2X,6E14.5）

RETURN

END

SUBROUTINESTRES（Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T）

DIMENSIONQ（MS）,QE（ME）,NX（NE,LE）,XY（LS,NN）,XYE（LS,NE）

DO400L=1,LE

DO160I=1,NE

DO160J=1,NF

N=NF*（I-1）+J

N1=NF*（NX（I,L）-1）+J

160QE（N）=Q（N1）

WRITE（*,165）L

WRITE（*,170）（QE（I）,I=1,ME）

165FORMAT（4X,'L=',I4）

170FORMAT（6E14.5）

DO200J=1,NE

LJ=NX（J,L）

DO200I=1,LS

200XYE（I,J）=XY（I,LJ）

CALLSTE（XYE,QE,NE,LS,ME,E,UM,T）

400CONTINUE

RETURN

END

SUBROUTINESTIFE（EK,XYE,ME,NE,NF,LS,E,UM,T）

DIMENSIONEK（ME,ME）,XYE（LS,NE）,B（3）,C（3）

B

（1）=XYE（2,2）-XYE（2,3）

B

（2）=XYE（2,3）-XYE（2,1）

B（3）=XYE（2,1）-XYE（2,2）

C

（1）=XYE（1,3）-XYE（1,2）

C

（2）=XYE（1,1）-XYE（1,3）

C（3）=XYE（1,2）-XYE（1,1）

AE=（B

（2）*C（3）-B（3）*C

（2））/2

D=E*T/4/（1-UM**2）/AE

DO30I=1,3

DO30J=1,3

M=2*（I-1）+1

N=2*（J-1）+1

EK（M,N）=D*（B（I）*B（J）+C（I）*C（J）*（1-UM）/2）

EK（M,N+1）=D*（UM*B（I）*C（J）+C（I）*B（J）*（1-UM）/2）

EK（M+1,N）=D*（UM*C（I）*B（J）+B（I）*C（J）*（1-UM）/2）

30EK（M+1,N+1）=D*（C（I）*C（J）+B（I）*B（J）*（1-UM）/2）

RETURN

END

SUBROUTINESTE（XYE,QE,NE,LS,ME,E,UM,T）

DIMENSIONXYE（LS,NE）,QE（ME）,SG（3）,B（3）,C（3）,S（3,6）

B

（1）=XYE（2,2）-XYE（2,3）

B

（2）=XYE（2,3）-XYE（2,1）

B（3）=XYE（2,1）-XYE（2,2）

C

（1）=XYE（1,3）-XYE（1,2）

C

（2）=XYE（1,1）-XYE（1,3）

C（3）=XYE（1,2）-XYE（1,1）

AE=（B

（2）*C（3）-B（3）*C

（2））/2

D=E*T/2/（1-UM**2）/AE

DO220I=1,3

S（1,2*I-1）=D*B（I）

S（2,2*I-1）=D*UM*B（I）

S（3,2*I-1）=D*（1-UM）*C（I）/2

S（1,2*I）=D*UM*C（I）

S（2,2*I）=D*C（I）

220S（3,2*I）=D*（1-UM）*B（I）/2

DO250I=1,3

SG（I）=0

DO250K=1,ME

250SG（I）=SG（I）+S（I,K）*QE（K）

WRITE（*,300）SG

300FORMAT（5X,'SEGMA',3E20.4）

RETURN

END

悬臂梁计算的输入数据文件：

2,3,55,14,10,80,2,

2.1e5,0.3,10,

400,80,400,60,400,40,400,20,400,0,

360,80,360,60,360,40,360,20,360,0,

320,80,320,60,320,40,320,20,320,0,

280,80,280,60,280,40,280,20,280,0,

240,80,240,60,240,40,240,20,240,0,

200,80,200,60,200,40,200,20,200,0,

160,80,160,60,160,40,160,20,160,0,

120,80,120,60,120,40,120,20,120,0,

80,80,80,60,80,40,80,20,80,0,

40,80,40,60,40,40,40,20,40,0,

0,80,0,60,0,40,0,20,0,0,

0,-200,0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-400,

0,0,0,0,0,0,0,0,0,-200,0,0,0,0,0,0,0,0,

1,7,2,2,8,3,3,9,4,4,10,5,1,6,7,2,7,8,3,8,9,4,9,10,

6,12,7,7,13,8,8,14,9,9,15,10,6,11,12,7,12,13,8,13,14,

9,14,15,11,17,12,12,18,13,13,19,14,14,20,15,11,16,17,

12,17,18,13,18,19,14,19,20,16,22,17,17,23,18,18,24,19,

19,25,20,16,21,22,17,22,23,18,23,24,19,24,25,21,27,22,

22,28,23,23,29,24,24,30,25,21,26,27,22,27,28,23,28,29,

24,29,30,26,32,27,27,33,28,28,34,29,29,35,30,26,31,32,

27,32,33,28,33,34,29,34,35,31,37,32,32,38,33,33,39,34,

34,40,35,31,36,37,32,37,38,33,38,39,34,39,40,36,42,37,

37,43,38,38,44,39,39,45,40,36,41,42,37,42,43,38,43,44,

39,44,45,41,47,42,42,48,43,43,49,44,44,50,45,41,46,47,

42,47,48,43,48,49,44,49,50,46,52,47,47,53,48,48,54,49,

49,55,50,46,51,52,47,52,53,48,53,54,49,54,55,

101,102,103,104,105,106,107,108,109,110,

程序输出结果文件：

2.58E+02

-1.37E+03

1.30E+02

-5.13E+02

1.05E+02

-2.42E+02

6.41E+01

-1.20E+02

4.37E+01

-7.32E+01

-4.13E+01

-1.35E+03

6.76E+01

-5.32E+02

6.46E+01

-2.54E+02

4.65E+01

-1.28E+02

4.21E+01

-8.12E+01

-6.18E+00

-1.36E+03

8.22E+01

-5.27E+02

6.54E+01

-2.53E+02

4.48E+01

-1.29E+02

4.06E+01

-8.21E+01

2.48E+00

-1.36E+03

8.64E+01

-5.26E+02

6.74E+01

-2.52E+02

4.51E+01

-1.28E+02

4.04E+01

-8.17E+01

4.53E+00

-1.36E+03

8.75E+01

-5.26E+02

6.81E+01

-2.51E+02

4.55E+01

-1.28E+02

4.06E+01

-8.14E+01

4.73E+00

-1.36E+03

8.76E+01

-5.25E+02

6.82E+01

-2.51E+02

4.56E+01

-1.28E+02

4.07E+01

-8.13E+01

3.65E+00

-1.36E+03

8.70E+01

-5.25E+02

6.78E+01

-2.51E+02

4.53E+01

-1.28E+02

4.06E+01

-8.12E+01

-6.27E-01

-1.35E+03

8.39E+01

-5.24E+02

6.61E+01

-2.51E+02

4.43E+01

-1.27E+02

3.98E+01

-8.10E+01

-1.53E+01

-1.35E+03

7.21E+01

-5.20E+02

5.94E+01

-2.48E+02

4.01E+01

-1.26E+02

3.65E+01

-7.98E+01

-4.41E+01

-1.27E+03

2.95E+01

-4.92E+02

3.53E+01

-2.34E+02

2.63E+01

-1.18E+02

2.21E+01

-7.52E+01

3.01E-06

-8.84E-06

-6.76E-07

-2.63E-07

-2.43E-07

-7.39E-08

-8.58E-08

-1.60E-08

-7.74E-08

-4.72E-09

程序二：

四节点矩形薄板单元程序

如图2为两对边简支、另两对边固支的矩形薄板。

其上受有均布压力q的作用，根据结构与荷载的对称性，可取其1/4为计算对象。

已知E=2.1×105N/mm2，

，厚度h=20mm。

将它划分为4×5的矩形网格，其节点及单元编号如图2。

建立如图2所示的离散化计算模型。

图2

计算源程序如下：

PROGRAMMAIN

DIMENSIONSK（90,21）,EK（12,12）,Q（90）,MC（35）,XY（2,100）,XYE（2,4）,QE（12）,NX（4,30）

OPEN（7,FILE='INP.DAT'）

REWIND7

READ（7,*）NF,NE,NN,MB,ND,LE,LS

READ（7,*）E,UM,T

10FORMAT（7I5）

12FORMAT（6F15.2）

WRITE（*,600）NF,NE,NN,MB,ND,LE,LS,E,UM,T

ME=NE*NF

MS=NN*NF

CALLINPUT（XY,Q,NX,MC,LS,NN,MS,NE,LE,ND）

WRITE（*,102）（（XY（I,J）,I=1,LS）,J=1,NN）

102FORMAT（10X,'XY'/,（2X,6F12.1））

WRITE（*,101）（Q（I）,I=1,MS）

101FORMAT（10X,'Q'/,（2X,3F20.3））

WRITE（*,500）（（NX（I,J）,I=1,NE）,J=1,LE）

WRITE（*,400）（MC（I）,I=1,ND）

500FORMAT（10X,'NX'/,（2X,12I6））

600FORMAT（10X,'NFNENNMBNDLELSEUMT'/7（2X,I4）,3（2X,F8.4））

400FORMAT（10X,'MC'/,（2X,10I6））

!

CALLQS（Q,TL,TH,P,LE,NE,MS,NX）

!

WRITE（*,120）（Q（I）,I=1,MS）

!

120FORMAT（10X,'Q'/（10X,3F20.3））

CALLSTIFS（SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,&

UM,T）

CALLSOLVE（SK,Q,MS,MB）

OPEN（9,FILE='OUT.DAT'）

REWIND9

WRITE（9,200）

WRITE（9,250）（Q（I）,I=1,MS）

200FORMAT（5X,'DISPLACEMENT'）

250FORMAT（2X,6E14.5）

CALLSTRES（Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T）

STOP1000

END

SUBROUTINEINPUT（XY,Q,NX,MC,LS,NN,MS,NE,LE,ND）

DIMENSIONXY（LS,NN）,Q（MS）,NX（NE,LE）,MC（ND）

READ（7,*）XY

READ（7,*）Q

READ（7,*）NX

READ（7,*）MC

CLOSE（7）

10FORMAT（6F11.2）

20FORMAT（12I5）

RETURN

END

!

SUBROUTINEQS（Q,TL,TH,P,LE,NE,MS,NX）

!

DIMENSIONNX（NE,LE）,Q（MS）

!

DOI=1,MS

!

Q（I）=0.0

!

ENDDO

!

DOJ=1,LE

!

I1=3*NX（1,J）-2

!

I2=3*NX（2,J）-2

!

I3=3*NX（3,J）-2

!

I4=3*NX（4,J）-2

!

Q（I1）=Q（I1）+P*TL*TH

!

Q（I1+1）=Q（I1+1）+P*TL*TH*TH/3

!

Q（I1+2）=Q（I1+2）-P*TL*TL*TH/3

!

Q（I2）=Q（I2）+P*TL*TH

!

Q（I2+1）=Q（I2+1）+P*TL*TH*TH/3

!

Q（I2+2）=Q（I2+2）+P*TL*TL*TH/3

!

Q（I3）=Q（I3）+P*TL*TH

!

Q（I3+1）=Q（I3+1）-P*TL*TH*TH/3

!

Q（I3+2）=Q（I3+2）+P*TL*TL*TH/3

!

Q（I4）=Q（I4）+P*TL*TH

!

Q（I4+1）=Q（I4+1）-P*TL*TH*TH/3

!

Q（I4+2）=Q（I4+2）-P*TL*TL*TH/3

!

ENDDO

!

RETURN

!

END

SUBROUTINESTIFS（SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,&

UM,T）

DIMENSIONSK（MS,MB）,EK（ME,ME）,NX（NE,LE）,MC（ND）,XY（LS,NN）,&

XYE（LS,NE）

DO35I=1,MS

DO35J=1,MB

35SK（I,J）=0

DO200L=1,LE

DO40J=1,NE

LJ=NX（J,L）

DO40I=1,LS

40XYE（I,J）=XY（I,LJ）

DO50I=1,ME

DO50J=1,ME

50EK（I,J）=0.0

CALLSTIFE（EK,XYE,ME,NE,LS,E,UM,T）

IF（L.EQ.1）WRITE（*,70）EK

70FORMAT（10X,'EK'/,（6E14.5））

DO200I=1,NE

DO200II=1,NF

M=NF*（I-1）+II

M1=NF*（NX（I,L）-1）+II

DO200J=1,NE

DO200JJ=1,NF

N=NF*（J-1）+JJ

N1=NF*（NX（J,L）-1）+JJ

MN=N1-M1+1

IF（MN）200,200,150

150SK（M1,MN）=SK（M1,MN）+EK（M,N）

200CONTINUE

DO220I=1,ND

M=MC（I）

220SK（M,1）=SK（M,1）*1E8

RETURN

END

SUBROUTINESOLVE（SK,Q,MS,MB）

DIMENSIONSK（MS,MB）,Q（MS）

K1=MS-1

DO125K=1,K1

IF（K+MB-1-MS）105,106,106

105N=K+MB-1

GOTO110

106N=MS

110I1=K+1

DO125I=I1,N

L=I-K+1

C=SK（K,L）/SK