有限元编程算例(fortran).docx

有限元编程算例(Fortran)
本程序通过Fortran语言编写，程序在Intel Parallel Studio XE 2013 with VS2013中成功运行，程序为《计算力学》（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。

源程序为：
!Page149
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)
OPEN(5,FILE='DATAIN')
!OPEN(6,FILE='DATAOUT',STATUS='NEW')
CALL DATA

IF(IND.EQ.0)GOTO 10
EO=EO/(1.0-UN*UN)
UN=UN/(1.0-UN)
10 CALL TOTSTI
CALL LOAD
CALL SUPPOR
CALL SOLVEQ
CALL STRESS
PAUSE
!STOP
END

SUBROUTINE DATA

COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

READ(5,*)NJ,NE,NZ,NDD,NPJ,IND
NJ2=NJ*2
NPJ1=NPJ+1
READ(5,*)EO,UN,GAMA,TE
READ(5,*)JM(I,J),J=1,3),I=1,NE)
READ(5,*)((CJZ(I,J),J = 1,2),I=1,NJ)

!Page150
READ(5,*)(NZC(I),I=1,NZ)
READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)
WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)
10 FORMAT(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7.2,F7.2))

RETURN

END

SUBROUTINE ELEST(MEO,IASK)
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

IE=JM(MEO,1)
JE=JM(MEO,2)
ME=JM(MEO,3)

CM=CJZ(JE,1)-CJZ(IE,1)
BM=CJZ(IE,2)-CJZ(JE,2)
CJ=CJZ(IE,1)-CJZ(ME,1)
BJ=CJZ(ME,2)-CJZ(IE,2)

AE=(BJ*CM-BM*CJ)/2.0
IFIASK.LE.1) GOTO 50

DO 10 I=1,3
DO 10 J=1,6
B(I,J)=0.0
10 CONTINUE

B(1,1)=-BJ-BM
B(1,3)=BJ
B(1,5)=BM
B(2,2)=-CJ-CM
B(2,4)=CJ
B(2,6)=CM
B(3,1)=B(2,2)
B(3,2)=B(1,1)
B(3,3)=B(2,4)
B(3,4)=B(1,3)
B(3,5)=B(2,6)
!Page151
B(3,6)=B(1,5)

DO 20 I=1,3
DO 20 J=1,6
B(I,J)=B(I,J)/(2.0*AE)
20 CONTINUE

D(1,1)=EO/(1.0-UN*UN)
D(1,2)=EO*UN/(1.0-UN*UN)
D(2,1)=D(1,2)
D(2,2)=D(1,1)
D(1,3)=0.0
D(2,3)=0.0
D(3,1)=0.0
D(3,2)=0.0
D(3,3)=EO/(2.0*(1.0+UN))

DO 30 I=1,3
DO 30 J=1,6
S(I,J)=0.0
DO 30 K=1,3
S(I,J)=S(I,J)+D(I,K)*B(K,J)
30 CONTINUE

IF(IASK.LE.2) GOTO 50

DO 40 I=1,6
DO 40 J=1,6
EKE(I,J)=0.0
DO 40 K=1,3
!**********************************Exchange B And S***********************************************
EKE(I,J)=EKE(I,J)+B(K,I)*S(K,J)*AE*TE
40 CONTINUE

50 CONTINUE

RETURN

END

SUBROUTINE TOTSTI

COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE

COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

!Page152
DO 20 I=1,NJ2
DO 20 J=1,NDD

TKZ(I,J)=0.0

20 CONTINUE
!*************Not Understanded***
DO 30 MEO=1,NE
CALL ELEST(MEO,3)
DO 30 I=1,3
DO 30 II=1,2
LH=2*(I-1)+II
LDH=2*(JM(MEO,I)-1)+II
DO 30 J=1,3
DO 30 JJ=1,2
L=2*(J-1)+JJ
LZ=2*(JM(MEO,J)-1)+JJ
LD=LZ-LDH+1
IF(LD.LE.0) GOTO 30
TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)
30 CONTINUE
RETURN
END

SUBROUTINE LOAD
COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

DO 10 I=1,NJ2
P(I)=0.0
10 CONTINUE
IF(NPJ.EQ.0) GOTO 30

DO 20 I=1,NPJ
I1=I+1
J=IFIX(PJ(I1,2))
P(J)=PJ(I1,1)
20 CONTINUE

30 IF(GAMA.LE.0.0) GOTO 50

!Page153
DO 40 MEO=1,NE
CALL ELEST(MEO,1)
PE=-GAMA*AE*TE/3.0
IE=JM(MEO,1)
JE=JM(MEO,2)
ME=JM(MEO,3)
P(2*IE)=P(2*IE)+PE
P(2*JE)=P(2*JE)+PE
P(2*ME)=P(2*ME)+PE
40 CONTINUE

50 CONTINUE

RETURN

END

SUBROUTINE SUPPOR

COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

DO 60 I=1,NZ
MZ=NZC(I)
TKZ(MZ,1)=1.0
DO 10 J=2,NDD
TKZ(MZ,J)=0.0
10 CONTINUE
IF(MZ-NDD)20,20,30
20 JO=MZ
GOTO 40
30 JO=NDD
40 DO 50 J = 2,JO
J1=MZ-J
TKZ(J1+1,J)=0.0
50 CONTINUE
P(MZ)=0.0
60 CONTINUE
RETURN

END

!Page154
SUBROUTINE SOLVEQ

COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

NJ1=NJ2-1
DO 50 K=1,NJ1
IF(NJ2-K-NDD+1)10,10,20
10 IM=NJ2
GOTO 30
20 IM=K+NDD-1
30 K1=K+1
DO 50 I=K1,IM
L=I-K+1
C=TKZ(K,L)/TKZ(K,1)
LD1=NDD-L+1
DO 40 J=1,LD1
M=J+I-K
TKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)
40 CONTINUE
P(I)=P(I)-C*P(K)

50 CONTINUE

P(NJ2)=P(NJ2)/TKZ(NJ2,1)

DO 100 I1 = 1,NJ1

I=NJ2-I1
!************************************************************************下面一行可能出错
IF(NDD-NJ2+I-1)60,60,70

60 JO=NDD

GOTO 80
70 JO=NJ2-I+1
80 DO 90 J=2,JO

LH=J+I-1

P(I)=P(I)-TKZ(I,J)*P(LH)
90 CONTINUE

P(I)=P(I)/TKZ(I,1)

100 CONTINUE
!Page155
WRITE(6,110)(I,P(2*I-1),P(2*I),I=1,NJ)
!************************************************************************************
110 FORMAT(2X,3HJD=,3X,2HU=,12X,2HV=/(I4,3X,F16.7,3X,F16.7))
RETURN

END

SUBROUTINE STRESS

COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EOUN,GAMA,TE,AE
COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)

DIMENSION WY(6),YL(3)

DO 60 MEO=1,NE

CALL ELEST(MEO,2)

DO 10 I=1,3

DO 10 J=1,2

LH=2*(I-1)+J

LDH=2*(JM(MEO,I)-1)+J

WY(LH)=P(LDH)

10 CONTINUE
DO 20 I=1,3

YL(I)=0.0

DO 20 J=1,6

YL(I)=YL(I)+S(I,J)*WY(J)

20 CONTINUE

SIGX=YL(1)

SIGY=YL(2)

TOXY=YL(3)

PYL=(SIGX+SIGY)/2.0

SIG=(SIGX-SIGY)**2/4.0+TOXY*TOXY

RYL=SQRT(SIG)

SIG1=PYL+RYL

SIG2=PYL-RYL

IF(SIGY.EQ.SIG2) GOTO 30

CETA1=TOXY/(SIGY-SIG2)

CETA=90.0-57.29578*ATAN(CETA1)

GOTO 40

!Page156
30 CETA=0.0

40 WRITE(6,50)MEO,SIGX,SIGY,TOXY,SIG1,SIG2,CETA

50 FORMAT(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)
!50 FORMAT(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)

60 CONTINUE
RETURN
END

输入文件为datain
28,36,9,10,4,0
1,0.17,0,1
1,5,2
2,5,6
2,6,3
3,6,7
3,7,4
4,7,8
5,9,6
6,9,10
6,10,7
7,10,11
7,11,8
8,11,12
9,13,10
10,13,14
10,14,11
11,14,15
11,15,12
12,15,16
13,17,14
14,17,18
14,18,15
15,18,19
15,19,16
16,19,20
17,21,18
18,21,22
18,22,19
19,22,23
19,23,20
20,23,24
21,25,22
22,25,26
22,26,23
23,26,27
23,27,24
24,27,28
0,6
1,6
2,6
3,6
0,5
1,5
2,5
3,5
0,4
1,4
2,4
3,4
0,3
1,3
2,3
3,3
0,2
1,2
2,2
3,2
0,1
1,1
2,1
3,1
0,0
1,0
2,0
3,0
7,15,23,31,39,47,49,50,55
0,0
-5E4,2
-10E4,4
-10E4,6
-5E4,8
输出结果为：DATAOUT
NO X Y
1 0.00 6.00
2 1.00 6.00
3 2.00 6.00
4 3.00 6.00
5 0.00 5.00
6 1.00 5.00
7 2.00 5.00
8 3.00 5.00
9 0.00 4.00
10 1.00 4.00
11 2.00 4.00
12 3.00 4.00
13 0.00 3.00
14 1.00 3.00
15 2.00 3.00
16 3.00 3.00
17 0.00 2.00
18 1.00 2.00
19 2.00 2.00
20 3.00 2.00
21 0.00 1.00
22 1.00 1.00
23 2.00 1.00
24 3.00 1.00
25 0.00 0.00
26 1.00 0.00
27 2.00 0.00
28 3.00 0.00
JD= U= V=
1 -29766.873 -1173917.750
2 -14003.185 -1174018.875
3 -3753.270 -1179518.125
4 0.000 -1181719.750
5 -26382.471 -1072681.500
6 -10746.993 -1073615.000
7 -2064.593 -1082360.750
8 0.000 -1085873.250
9 -13536.995 -964010.125
10 3372.794 -970055.125
11 7268.415 -989269.125
12 0.000 -998401.812
13 7816.581 -835383.438
14 27176.234 -861713.938
15 22063.230 -905726.125
16 0.000 -927165.188
17 29514.479 -665602.875
18 53419.637 -747340.438
19 34876.832 -839806.812
20 0.000 -881219.125
21 29580.273 -416288.719
22 52944.918 -632601.125
23 17504.195 -803765.688
24 0.000 -859481.938
25 0.000 0.000
26 -120102.820 -583505.375
27 -76202.375 -787347.188
28 0.000 -829170.812
E= 1
SX= -1489.530 SY=-101489.383 TAU= -1489.531
S1= -1467.348 S2=-101511.562 CET= 179.147
E= 2
SX= -1475.844 SY=-100654.875 TAU= -1790.500
S1= -1443.531 S2=-100687.188 CET= 178.966
E= 3
SX= -7021.670 SY=-101597.672 TAU= -3741.688
S1= -6873.875 S2=-101745.469 CET= 177.738
E= 4
SX= -8067.500 SY= -98528.750 TAU= -4459.156
S1= -7848.227 S2= -98748.023 CET= 177.185
E= 5
SX= -13143.328 SY= -99391.750 TAU= -1662.500
S1= -13111.293 S2= -99423.781 CET= 178.896
E= 6
SX= -14652.781 SY= -98337500 TAU= -1501.062
S1= -14625.867 S2= -98364.414 CET= 178.973
E= 7
SX= -2923.122 SY=-109168.297 TAU= -5888.469
S1= -2597.762 S2=-109493.656 CET= 176.837
E= 8
SX= -716.078 SY=-103681.562 TAU= -8617.406
S1= 0.148 S2=-104397.789 CET= 175.249
E= 9
SX= -9188.316 SY=-105121.867 TAU= -9771.594
S1= -8203.125 S2=-106107.062 CET= 174.243
E= 10
SX= -12285.000 SY= -95180.250 TAU= -12199.594
S1= -10526.887 S2= -96938.359 CET= 171.799
E= 11
SX= -14170.516 SY= -95500.750 TAU= -5489.531
S1= -13801.664 S2= -95869.602 CET= 176.156
E= 12
SX= -22797.406 SY= -91347.000 TAU= -3902.844
S1= -22575.914 S2= -91568.492 CET= 176.752
E= 13
SX= -5104.269 SY=-129494.438 TAU= -11708.750
S1= -4011.727 S2=-130586.977 CET= 174.669
E= 14
SX= 969.672 SY=-108176.375 TAU= -21424.750
S1= 5024.582 S2=-112231.281 CET= 169.283
E= 15
SX= -14954.572 SY=-110883.469 TAU= -18383.531
S1= -11552.273 S2=-114285.766 CET= 169.515
E= 16
SX= -19890.141 SY= -86924.312 TAU= -25131.188
S1= 11514.844 S2= -95299.609 CET= 161.569
E= 17
SX= -22109.688 SY= -87301.625 TAU= -10225.406
S1= -20543.453 S2= -88867.859 CET= 171.292
E= 18
SX= -35190.453 SY= -77219.000 TAU= -9162.000
S1= -33280.023 S2= -79129.430 CET= 168.222
E= 19
SX= -9785.850 SY=-171444.172 TAU= -20524.969
S1= -7220.594 S2=-174009.422 CET= 172.876
E= 20
SX= 4594.438 SY=-113592.375 TAU= -46145.688
S1= 20477.398 S2=-129475.336 CET= 161.007
E= 21
SX= -25287.307 SY=-118672.312 TAU= -30023.750
S1= -16467.512 S2127492.109 CET= 163.629
E= 22
SX= -30634.422 SY= -71127.188 TAU= -44991.469
S1= -1543.715 S2=-100217.891 CET= 147.114
E= 23
SX= -34259.609 SY= -71743.438 TAU= -14637.906
S1= -29220.699 S2= -76782.344 CET= 161.005
E= 24
SX= -43958.047 SY= 53418.938 TAU= -17697.562
S1= -30369.627 S2= -67007.359 CET= 142.482
E= 25
SX= -19028.160 SY=-252549.000 TAU= -34958.688
S1= -13907.055 S2=-257670.094 CET= 171.666
E= 26
SX= 3973.812 SY=-114063.750 TAU= -92238.344
S1= 54459.047 S2=-164548.984 CET= 151.307
E= 27
SX= -39180.809 SY=-121400.055 TAU= -39312.688
S1= -23409.074 S2=-137171.781 CET= 158.140
E= 28
SX= -42804.766 SY= -43317.938 TAU= -65723.062
S1= 22662.211 S2=-108784.914 CET= 135.112
E= 29
SX= -42224.094 SY= -43219.188 TAU= -10273.375
S1= -32436.225 S2= -53007.055 CET= 136.386
E= 30
SX= -21830.422 SY= -25448.312 TAU= -23810.344
S1= 239.594 S2= -47518.328 CET= 137.172
E= 31
SX= -48815.199 SY=-424587.344 TAU= -79800.078
S1= -32570.844 S2=-440831.688 CET= 168.494
E= 32
SX=-132271.750 SY= -71582.000 TAU=-175409.250
S1= 76087.781 S2=-279941.531 CET= 130.093
E= 33
SX= -45090.102 SY= -56761.105 TAU= 804.781
S1= -45034.867 S2= -56816.336 CET= 3.926
E= 34
SX= 42332.711 SY= -9221.938 TAU= -47066.344
S1= 70218.328 S2= -37107.555 CET= 149.354
E= 35
SX= -20899.344 SY= -19971.375 TAU= 16235.219
S1= -4193.512 S2= -36677.207 CET= 45.819
E= 36
SX= 73163.914 SY= -17873.250 TAU= -17873.344
S1= 76547.250 S2= -21256.586 CET= 169.281

Welcome To
Download !!!

欢迎您的下载，资料仅供参考！

精品资料