有限元中三角形单元源程序

56附录：平面问题三角形单元源程序********************************************************************ANALYSISPROGTAMOFFINITEELEMENTMETHOD**FORPLANESTRESS/STRAINOFTRIANGULARELEMENT**-----FEMT3.FOR-----**-------------------------------------------------------------**Subroutines:1-SDATA,2-STE,3-ATE,4-DTE,5-BTE,6-STIFF**7-EQUPE,8-INSCD,9-BGSMT,10-SIGME********************************************************************DIMENSIONLND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REALKS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*)NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)'FINITEELEMENTANALYSISINPLANEPROBLEM'WRITE(6,*)'SOURCEDATAOUTPUT'WRITE(6,20)NJ,NE,NS,NPJ,IPS20FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0)WRITE(6,*)'PLANESTRESSPROBLEM'IF(IPS.EQ.1)WRITE(6,*)'PLANESTRAINPROBLEM'CALLSDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50)NJ250FORMAT(/1X,'DEGREESOFFREEDOM=',I5)WRITE(6,60)NW60FORMAT(1X,'BANDWIDTH=',I5)CALLSTIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALLEQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALLINSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALLBGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALLSIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------CSUBPROGRAM-1CINPUTSTRUCTURALDATASUBROUTINESDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,*T,V,LND,X,Y,JR,PJ)DIMENSIONLND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*)E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10)E,PR,T,V10FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)57READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20FORMAT(/1X,'ELEMENTINFORMATION'/3X,'ELEM',3X,*'IJK'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40FORMAT(/1X,'COORDINATESOFNODES'/3X,'NODES',*8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60FORMAT(/1X,'CONSTRAINEDNODES'/3X,'NODE',3X,'X',4X,'Y')WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80FORMAT(/1X,'LOADCASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90FORMAT(1X,F5.0,2E15.6)100NW=0（半带宽）DO110IE=1,NEDO110I=1,3DO110J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW)THENNW=IWENDIF110CONTINUENW=(NW+1)*2IF(IPS.NE.0)THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------CSUBPROGRAM-2CCALCULATEELEMENTSTIFFNESSMATRIXSUBROUTINESTE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSIONLND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REALKE(6,6)CALLATE(IE,NJ,NE,LND,X,Y,AE)CALLDTE(E,PR,D)58CALLBTE(IE,NJ,NE,LND,X,Y,AE,B)DO10I=1,6DO10J=1,6KE(I,J)=0.DO10K=1,3DO10K1=1,310KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO30I=1,6DO30J=1,630KE(I,J)=KE(I,J)*CEND*------------------------------------------------CSUBPROGRAM-3CCALCULATEELEMENTAREASUBROUTINEATE(IE,NJ,NE,LND,X,Y,AE)DIMENSIONLND(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*----------------------------------------------CSUBPROGRAM-4CCALCULATEELASTICITYMATRIXSUBROUTINEDTE(E,PR,D)DIMENSIOND(3,3)DO10I=1,3DO10J=1,310D(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*------------------------------------------------CSUBPROGRAM-5CCALCULATEMATRIX[B]SUBROUTINEBTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSIONLND(NE,3),X(NJ),Y(NJ),B(3,6)59I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO10II=1,3DO10JJ=1,610B(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)DO20I1=1,3DO20J1=1,620B(I1,J1)=.5/AE*B(I1,J1)END*-------------------------------------------------------CSUBPROGRAM-6CCALCULATEGLOBALSTIFFNESSMATRIXSUBROUTINESTIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)DIMENSIONLND(NE,3),X(NJ),Y(NJ)REALKS(NJ2,NW),KE(6,6)DO5I=1,NJ2DO5J=1,NW5KS(I,J)=0.DO10IE=1,NECALLSTE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO10I=1,3IZ=LND(IE,I)DO10II=1,2IH=2*(I-1)+IIIDH=2*(IZ-1)+IIDO10J=1,3JZ=LND(IE,J)DO10JJ=1,2L=2*(J-1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH)THENIDL=IL-IDH+160KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10CONTINUEEND*--------------------------------------------------------CSUBPROGRAM-7CCALCULATENODALLOADVECTORSUBROUTINEEQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)DIMENSIONLND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2)DO10I=1,NJ210P(I)=0.DO20I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20P(2*II)=PJ(I,3)30IF(V.EQ.0.)GOTO50DO40IE=1,NECALLATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO40I=1,3II=LND(IE,I)40P(2*II)=P(2*II)+PE50RETURNEND*---------------------------------------------CSUBPROGRAM-8CINTRODUCEBOUNDARYCONDITIONSUBROUTINEINSCD(NS,NW,NJ2,JR,KS,P)DIMENSIONP(NJ2),JR(NS,3)REALKS(NJ2,NW)DO30I=1,NSIR=JR(I,1)DO30J=2,3IF(JR(I,J).EQ.0)GOTO30II=2*IR+J-3KS(II,1)=1.DO10JJ=2,NW10KS(II,JJ)=0.IF(II.GT.NW)JO=NWIF(II.LE.NW)JO=IIDO20JJ=2,JO20KS(II-JJ+1,JJ)=0.P(II)=0.30CONTINUEEND61*-------------------------------------------CSUBPROGRAM-9CSOLVEEQUATIONSBYGSMETHODSUBROUTINEBGSMT(NJ,NJ2,NW,KS,P)DIMENSIONP(NJ2)REALKS(NJ2,NW)NJ1