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

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