有限元编程算例(fortran)

有限元编程算例

(Fortran)

本程序通 过

Fortran

语言编写，程序在

I ntel Parallel Studio XE 2013 with

VS201 3

中成功运行，程序为《计算力学》

（龙述尧等编）一书中的源 程序，仅

作研究学习使用，省去了敲写的麻烦。

源程序为：

!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(.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(1 00,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 FORMA

T(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,U N,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(20 0,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

IF(.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(.2) GOTO 50

DO 40 I=1,6

DO 40 J=1,6

EKE(I,J)=0.0

DO 40 K=1,3

!**********************************Exchang e

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,N DD,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(20 0)

!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

B

And

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(.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(.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(.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,N DD,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(1 00,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 FORMA

T(2X,3HJD=,3X,2HU=, 12X,2HV=/(I4,3X,F16.7,3X,F16.7))

RETURN

END

SUBROUTINE STRESS