用常应变三角形单元解弹性力学平面问题的程序

用常应变三角形单元解弹性力学平面问题

的程序

******************************************************************* * ANALYSIS PROGTAM OF FINITE ELEMENT METHOD * * FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT * * ----- FEMT3.FOR ----- * *------------------------------------------------------------- * * Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF * * 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME * ******************************************************************* DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200) REAL KS(200,100)

OPEN(5,FILE='FEMT3.DAT')

OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')

READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型) WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM' WRITE(6,*)' SOURCE DATA OUTPUT'

WRITE(6,20) NJ,NE,NS,NPJ,IPS

20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5) IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM' IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'

CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ) NJ2=2*NJ

WRITE(6,50) NJ2

50 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)

WRITE(6,60) NW

60 FORMAT(1X,'BAND WIDTH=',I5)

CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6) CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7) CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8) CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）

CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10） CLOSE(5)

CLOSE(6)

END

*--------------------------------------------------------

C SUBPROGRAM-1

C INPUT STRUCTURAL DATA

SUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,

* T,V,LND,X,Y,JR,PJ)

DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3) READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)

WRITE(6,10) E,PR,T,V

10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2) READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码） WRITE(6,20)

20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X, * 'I J K'/)

WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)

30 FORMAT(1X,4I5)

READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）

WRITE(6,40)

40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES', * 8X,'X',13X,'Y')

WRITE(6,50)(I,X(I),Y(I),I=1,NJ)

50 FORMAT(1X,I5,2E15.6)

READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息） WRITE(6,60)

60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)

70 FORMAT(1X,3I5)

READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息） WRITE(6,80)

80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y') WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)

90 FORMAT(1X,F5.0,2E15.6)

100 NW=0（半带宽）

DO 110 IE=1,NE

DO 110 I=1,3

DO 110 J=1,3

IW=IABS(LND(IE,I)-LND(IE,J))

IF(NW.LT.IW) THEN

NW=IW

ENDIF

110 CONTINUE

NW=(NW+1)*2

IF(IPS.NE.0) THEN

E=E/(1.0-PR*PR)

PR=PR/(1.0-PR)

ENDIF

END

*---------------------------------------------------------

C SUBPROGRAM-2

C CALCULATE ELEMENT STIFFNESS MATRIX

SUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)

DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3) REAL KE(6,6)

CALL ATE(IE,NJ,NE,LND,X,Y,AE)

CALL DTE(E,PR,D)

CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)

DO 10 I=1,6

DO 10 J=1,6

KE(I,J)=0.

DO 10 K=1,3

DO 10 K1=1,3

10 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)

C=AE*T

DO 30 I=1,6

DO 30 J=1,6

30 KE(I,J)=KE(I,J)*C

END

*------------------------------------------------

C SUBPROGRAM-3

C CALCULATE ELEMENT AREA

SUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)

DIMENSION LND(NE,3),X(NJ),Y(NJ)

I=LND(IE,1)

J=LND(IE,2)

K=LND(IE,3)

XIJ=X(J)-X(I)

YIJ=Y(J)-Y(I)

XIK=X(K)-X(I)

YIK=Y(K)-Y(I)

AE=.5*(XIJ*YIK-XIK*YIJ)

END

*----------------------------------------------

C SUBPROGRAM-4

C CALCULATE ELASTICITY MATRIX

SUBROUTINE DTE(E,PR,D)

DIMENSION D(3,3)

DO 10 I=1,3

DO 10 J=1,3

10 D(I,J)=0.

D(1,1)=E/(1.-PR*PR)

D(1,2)=E*PR/(1.-PR*PR)

D(2,1)=D(1,2)

D(2,2)=D(1,1)

D(3,3)=.5*E/(1.+PR)

END

*------------------------------------------------

C SUBPROGRAM-5

C CALCULATE MATRIX [B]

SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)

DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)

I=LND(IE,1)

J=LND(IE,2)

K=LND(IE,3)

DO 10 II=1,3

DO 10 JJ=1,6

10 B(II,JJ)=0.

B(1,1)=Y(J)-Y(K)

B(1,3)=Y(K)-Y(I)

B(1,5)=Y(I)-Y(J)

B(2,2)=X(K)-X(J)

B(2,4)=X(I)-X(K)

B(2,6)=X(J)-X(I)

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)

B(3,6)=B(1,5)

DO 20 I1=1,3

DO 20 J1=1,6

20 B(I1,J1)=.5/AE*B(I1,J1)

END

*-------------------------------------------------------

C SUBPROGRAM-6

C CALCULATE GLOBAL STIFFNESS MATRIX

SUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)

REAL KS(NJ2,NW),KE(6,6)

DO 5 I=1,NJ2

DO 5 J=1,NW

5 KS(I,J)=0.

DO 10 IE=1,NE

CALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)

DO 10 I=1,3

IZ=LND(IE,I)

DO 10 II=1,2

IH =2*(I -1)+II

IDH=2*(IZ-1)+II

DO 10 J=1,3

JZ=LND(IE,J)

DO 10 JJ=1,2

L =2*(J -1)+JJ

IL=2*(JZ-1)+JJ

IF(IL.GE.IDH) THEN

IDL=IL-IDH+1

KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)

ENDIF

10 CONTINUE

END

*--------------------------------------------------------

C SUBPROGRAM-7

C CALCULATE NODAL LOAD VECTOR

SUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ2

10 P(I)=0.

DO 20 I=1,NPJ

II=PJ(I,1)

P(2*II-1)=PJ(I,2)

20 P(2*II)=PJ(I,3)

30 IF(V.EQ.0.) GOTO 50

DO 40 IE=1,NE

CALL ATE(IE,NJ,NE,LND,X,Y,AE)

PE=-V*AE*T/3.

DO 40 I=1,3

II=LND(IE,I)

40 P(2*II)=P(2*II)+PE

50 RETURN

END

*---------------------------------------------

C SUBPROGRAM-8

C INTRODUCE BOUNDARY CONDITION

SUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)

DIMENSION P(NJ2),JR(NS,3)

REAL KS(NJ2,NW)

DO 30 I=1,NS

IR=JR(I,1)

DO 30 J=2,3

IF(JR(I,J).EQ.0) GOTO 30

II=2*IR+J-3

KS(II,1)=1.

DO 10 JJ=2,NW

10 KS(II,JJ)=0.

IF(II.GT.NW) JO=NW

IF(II.LE.NW) JO=II

DO 20 JJ=2,JO

20 KS(II-JJ+1,JJ)=0.

P(II)=0.

30 CONTINUE

END

*-------------------------------------------

C SUBPROGRAM-9

C SOLVE EQUATIONS BY GS METHOD SUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)

DIMENSION P(NJ2)

REAL KS(NJ2,NW)

NJ1=NJ2-1

DO 20 K=1,NJ1

IF(NJ2.GT.K+NW-1) IM=K+NW-1

IF(NJ2.LE.K+NW-1) IM=NJ2

K1=K+1

DO 20 I=K1,IM

L=I-K+1

C=KS(K,L)/KS(K,1)

IW=NW-L+1

DO 10 J=1,IW

M=J+I-K

10 KS(I,J)=KS(I,J)-C*KS(K,M)

20 P(I)=P(I)-C*P(K)

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

DO 40 I1=1,NJ1

I=NJ2-I1

IF(NW.GT.NJ2-I+1) JO=NJ2-I+1

IF(NW.LE.NJ2-I+1) JO=NW

DO 30 J=2,JO

K=J+I-1

30 P(I)=P(I)-KS(I,J)*P(K)

40 P(I)=P(I)/KS(I,1)

WRITE(6,50)

50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X, * 'NODE',5X,'X-DISP.',8X,'Y-DISP.')

DO 60 I=1,NJ

60 WRITE(6,70) I,P(2*I-1),P(2*I)

70 FORMAT(1X,I5,2E15.6)

END

*---------------------------------------------------

C SUBPROGRAM-10

C CALCULATE ELEMENT STRESS MATRIX SUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)

DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)

WRITE(6,*)

WRITE(6,*)' ELEMENT STRESSES'

CALL DTE(E,PR,D)

DO 50 IE=1,NE

CALL ATE(IE,NJ,NE,LND,X,Y,AE)

CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)

DO 10 I=1,3

DO 10 J=1,6

S(I,J)=0.

DO 10 K=1,3

10 S(I,J)=S(I,J)+D(I,K)*B(K,J)

DO 20 I=1,3

DO 20 J=1,2

IH=2*(I-1)+J

IW=2*(LND(IE,I)-1)+J

20 DE(IH)=P(IW)

DO 30 I=1,3

ST(I)=0.

DO 30 J=1,6

30 ST(I)=ST(I)+S(I,J)*DE(J)

SGX=ST(1)

SGY=ST(2)

TXY=ST(3)

ASG=(SGX+SGY)*.5

RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)

SGMA=ASG+RSG

SGMI=ASG-RSG

IF(SGY.EQ.SGMI) CETA=0.

IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN * (TXY/(SGY-SGMI))

50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA 60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4) END

本文来源：https://www.bwwdw.com/article/ugs1.html