2022年有限元中三角形单元源程序 .pdf

《2022年有限元中三角形单元源程序 .pdf》由会员分享，可在线阅读，更多相关《2022年有限元中三角形单元源程序 .pdf（7页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、56 附录：平面问题三角形单元源程序* * 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

2、) 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

3、,*) 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

4、,NJ2,T,V,LND,X,Y,PJ,P)(P7) 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),

5、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) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 7 页 - - - - - - - - - 57 READ(5,*)(LND(I,J),J=1,3),I=1,NE)（结点编码） WRITE(6,20) 20 FORMAT(/1X,

6、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)（约束信息） W

7、RITE(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,

8、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)

9、,B(3,6),D(3,3) REAL KE(6,6) CALL ATE(IE,NJ,NE,LND,X,Y,AE) CALL DTE(E,PR,D) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 7 页 - - - - - - - - - 58 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

10、,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*YI

11、J) 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

12、,B) DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 7 页 - - - - - - - - - 59 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

13、(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(N

14、J),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 名师资料总结 -

15、 - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 7 页 - - - - - - - - - 60 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(N

16、J2) 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(N

17、S,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 名师资料总结 - - -精品资料欢迎下载 - - -

18、- - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 7 页 - - - - - - - - - 61 *- 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

19、-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,NO

20、DAL 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) 名师资

21、料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 7 页 - - - - - - - - - 62 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,

22、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) C

23、ETA=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名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 7 页 - - - - - - - - -