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

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]