平面四边形四节点等参单元Fortran源程序

1、标准*FINITEELEMENTPROGRAM*FORTwoDIMENSIONALELASticityPROBLEM文案WITH4NODE*PROGRAMELASTICITYcharacter*32dat,cchDIMENSIONSK(80000),COOR(2,300),AE(4,11),MEL(5,200),&WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)WRITE(*,*)'

2、;PLEASEENTERINPUTFILENAME'READ(*,'(A)')DATOPEN(4,FILE=dat,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBEROFNODENP=',npWRITE(7,'(A,I6)')'NUMBEROFELEMENTNE=',neWRITE(7,'(A,I6)')'

3、;NUMBEROFMATERIALNM=',nmWRITE(7,'(A,I6)')'NUMBEROFsurportingNC=',NrCALLINPUT(JR,COOR,AE,MEL)CALLCBAND(MA,JR,MEL)DOI=1,NHSK(I)=0.0enddoCALLSK0(SK,MEL,COOR,JR,MA,AE)doI=1,NR(I)=0.0enddopause'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)'

4、;)NCP,NBE,izIF(NCP.GT.0)CALLCONCR(NCP,R,JR)IF(NBE.GT.0)CALLBODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendojj=1,izREAD(4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALLFACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendif标准CALLDECOP(SK,MA)CALLFOBA(SK,MA,R)CALLOUTDISP(NP,R,JR)CALLSTRESS(COOR,MEL,JR,AE,R,STRE)WRI

5、TE(7,'(A)')'PROGRAMSAFFHASBEENENDED'WRITE(*,'(A)')'PROGRAMSAFFHASBEENENDED'STOPcRETURNENDC*SUBROUTINEINPUT(JR,COOR,AE,MEL)DIMENSIONJR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHDO70I=1,NPREAD(4,*)IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70CONTINUEDO

6、11J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11CONTINUEDO10I=1,NPDO10J=1,210JR(J,I)=1DO20I=1,NRREAD(4,*)IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20CONTINUEN=0DO30I=1,NPDO30J=1,2IF(JR(J,I)30,30,2525N=N+1JR(J,I)=N30CONTINUEDO55J=1,NMREAD(4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910)JJ,(AE(I,JJ),I=1,4)55CONTINU

7、E910FORMAT(/20X,'MATERIALPROPERTIES'/(3X,I5,4(1x,E8.3)RETURN文案标准ENDC*SUBROUTINECBAND(MA,JR,MEL)DIMENSIONMA(*),JR(2,*),MEL(5,*),NN(8)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHDO65I=1,N65MA(I)=0DO90IE=1,NEDO75K=1,4IEK=MEL(K,IE)DO95M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95CONTINUE75CONTINUEL=NDO80I=1

8、,2*4NNI=NN(I)IF(NNI.EQ.0)GOTO80IF(NNI.LT.L)L=NNI80CONTINUEDO85M=1,2*4JP=NN(M)IF(JP.EQ.0)GOTO85JPL=JP-L+1IF(JPL.GT.MA(JP)MA(JP)=JPL85CONTINUE90CONTINUEMX=0MA(1)=1DO10I=2,NIF(MA(I).GT.MX)MX=MA(I)MA(I)=MA(I)+MA(I-1)10CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTALDEGREESOFFREEDOMN=',NWRITE(7,&#

9、39;(A,I8)')'MAX-SEMI-BANDWIDTHMX=',MXWRITE(7,'(A,I8)')'TOTAL-STORAGENH=',NH500FORMAT(/5X,'FREEDOMN='*,I5,3X,'SEMI-BANDWI.MX=',I5,3X,*'STORAGENH=',I7)RETURNEND文案标准C*SUBROUTINESK0(SK,MEL,COOR,JR,MA,AE)DIMENSIONSK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*),*

10、AE(4,*),XYZ(2,4),iven(4)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/GAUSS/RSTG(3),H(3)H(1)=0.5555555555555560H(2)=0.8888888888888890H(3)=H(1)RSTG(1)=-0.7745966692414830RSTG(2)=0.00RSTG(3)=-RSTG(1)DO10IE=1,NENEE=IENME=MEL(5,IE)DO75K=1,4IEK=MEL(K

11、,IE)iven(k)=IEKDO95M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95XYZ(M,K)=COOR(M,IEK)75CONTINUECALLSTIF(XYZ,AE,iven)DO60I=1,8DO60J=1,8II=NN(I)JJ=NN(J)IF(JJ.EQ.0).OR.(II.LT.JJ)GOTO60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60CONTINUE70CONTINUEwrite(7,1111)(ske(i,j),j=1,8),i=1,8)1111format(2x,8f12.2)10CONTINUERETU

12、RNENDC*SUBROUTINESTIF(XYZ,AE,iven)DIMENSIONAE(4,*),DNX(2,4),XYZ(2,*),iven(*),文案标准*RJAC(2,2)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/GAUSS/RSTG(3),H(3)DO40I=1,8RF(I)=0.00DO30J=1,8SKE(I,J)=0.0030CONTINUE40CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(

13、3,NME)D1=E*(1.00-U)/(1.00+U)*(1.00-2.00*U)D2=E*U/(1.00+U)*(1.00-2.00*U)D3=E*0.50/(1.00+U)DO120I=1,4II=2*(I-1)I1=II+1I2=II+2DO115J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO99IS=1,3S=RSTG(IS)SH=H(IS)DO98IR=1,3R=RSTG(IR)RH=H(IR)CALLFDNX(XYZ,DNX,DET,R,S,RJAC,iven,NEE)DNIX=DNX(1,I)DNIY=DNX(2,I)DN

14、JX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SH文案标准DYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115CONTINUE120CONTINUERETURNENDC*SUBROUTINECON

15、CR(NCP,R,JR)DIMENSIONR(*),JR(2,*),XYZ(2)DO100I=1,NCPREAD(4,*)IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO95J=1,2L=JR(J,IP)IF(L.EQ.0)GOTO95R(L)=R(L)+XYZ(J)95CONTINUE100CONTINUERETURNENDC*SUBROUTINEBODYR(NBE,R,MEL,COOR,JR,AE)DIMENSIONR(*),MEL(5,*),COOR(2,*),JR(2,*),&AE(4,*),XYZ(2,4),iven(4)COMMON/CMN1/NP,NE,NM,NR

16、COMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN5/FUN(4),PN(2,4),XJAC(2,2)COMMON/GAUSS/RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)DO10IE=1,NBEDOI=1,8RF(I)=0.00ENDDOcREAD(4,*)NEE文案标准NEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO75K=1,4IEK=MEL(K,NEE)iven(k)=iekDO95M=1,2J

17、J=2*(K-1)+MNN(JJ)=JR(M,IEK)95XYZ(M,K)=COOR(M,IEK)75CONTINUEDO99IS=1,2S=RSTG(IS)SH=H(IS)DO98IR=1,2RR=RSTG(IR)RH=H(IR)CALLFUN8(XYZ,RR,S,DET)DO30I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30CONTINUE98 CONTINUE99 CONTINUECALLASLOAD(R)100 CONTINUERETURNENDC*SUBROUTINEFACER(iew,NSE,R,MEL,COOR,JR,WG)DIMENS

18、IONR(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*)*,XYZ(2,4),iew(*),PR(2)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/GAUSS/RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0.1文案标准if(ir.eq.1)nwf=1if(ir.eq.2)nnf=1DO

19、510IE=1,NSEDOI=1,8RF(I)=0.00ENDDOnee=iew(ie)DO575K=1,4IEK=MEL(K,NEE)DO595M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595XYZ(M,K)=COOR(M,IEK)575CONTINUEIF(NWF.EQ.1)thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0.1CALLSURLOD(NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF.EQ.1)thenq=WG(2)NSU=WG(4)+0.1doj=1,2PR(J)=qenddoCALLSURLOD(NSU,XYZ,PR

20、,Z0,GAMA,2)endifCALLASLOAD(R)510CONTINUERETURNENDC*SUBROUTINESURLOD(NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSIONXYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4),&FVAL(4),NODES(2),FACT(4)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN4/NEE,NMECOMMON/CMN5/FUN(4),PN(2,4),XJAC(2,2)C

21、OMMON/GAUSS/RSTG(3),H(3)DATAKCRD/1,1,2,2/DATAKFACE/1,4,文案2,1,4,3,2,3/标准25306070文案DATAFVAL/-1.00,1.00,-1.00,1.00/FACT(1)=1.0FACT(2)=-1.0FACT(3)=-1.0FACT(4)=1.0FACTNUS=FACT(NSU)DOI=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF(NSI.EQ.1)THENDOI=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0.00IF(Z.GT.0.00)PR(I)=Z*GAMAENDDOEND

22、IFML=KCRD(NSU)IF(ML.EQ.1)MM=2IF(ML.EQ.2)MM=1RST(ML尸FVAL(NSU)DO70LX=1,2RST(MM尸RSTG(LX)CALLFUN8(XYZ,RST(1),RST(2),DET)PXYZ=0.00DO25I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)DO60I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A2CONT

23、INUECONTINUE标准RETURNENDC*SUBROUTINEASLOAD(R)DIMENSIONR(*)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN3/RF(8),SKE(8,8),NN(8)DO20I=1,8L=NN(I)IF(L.EQ.0)GOTO20R(L)=R(L)+RF(I)20CONTINUERETURNENDC*SUBROUTINEDECOP(SK,MA)DIMENSIONSK(*),MA(*)COMMON/CMN2/N,MX,NHDO50I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF(L1.GT.K)GOTO30DO2

24、0J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF(L.GT.M)M=LMP=J-1IF(M.GT.MP)GOTO20DO10LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10CONTINUE20CONTINUE30IF(L.GT.K)GOTO50DO40LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP尸SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40CONTINUE50CONTINUE文案标准RETU

25、RNENDC*SUBROUTINEFOBA(SK,MA,R)DIMENSIONSK(*),MA(*),R(*)COMMON/CMN2/N,MX,NHDO10I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF(L.GT.K)GOTO10DO5LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5CONTINUE10CONTINUEDO20I=1,NII=MA(I)45R(I)=R(I)/SK(II)20CONTINUEDO30J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF(L.GT.K)GOTO30DO25J=L,KIJ

26、=MA(I)-I+J55R(J)=R(J)-SK(IJ)*R(I)25CONTINUE30CONTINUERETURNENDC*SUBROUTINESTRESS(COOR,MEL,JR,AE,R,STRE)DIMENSIONXYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),&COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),&B(3,8),R(*),iven(4)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN3/RF(8),SKE(8,8),NN(8)COMMON/CMN5/FUN(4),PN(2,

27、4),XJAC(2,2)DO106IE=1,NENME=MEL(5,IE)DO300K=1,4IEK=MEL(K,IE)DO310M=1,2文案标准310XYZ(M,K)=COOR(M,IEK)DO320M=1,2JRR=2*(K-1)+M320NN(JRR)=JR(M,IEK)300CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1.00-U)/(1.00+U)*(1.00-2.00*U)D2=E*U/(1.00+U)*(1.00-2.00*U)D3=0.50*E/(1.00+U)SS=0.0RR=0.0CALLFDNX(XYZ,DNX,DET,RR,SS,RJAC,

28、iven,IE)DO30I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0.B(3,J1)=CIB(1,J2)=0.B(2,J2)=CIB(3,J2)=BI30CONTINUEDO55II=1,3SIG(II)=0.0055CONTINUEDO70K=1,8NA=NN(K)IF(NA.EQ.0)GOTO70DO60L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60CONTINUE70CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SX

29、Y=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106CONTINUECALLOUTSTRE(NE,STRE)文案标准RETURNENDC*SUBROUTINEFDNX(XYZ,DNX,DET,R,S,RJAC,iven,NEE)DIMENSIONXYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)COMMON/CMN5/FUN(4),PN(2,4),XJAC(2,2)CALLFUN8(XYZ,R,S,DET)IF(DET.LT.1.0E-5)THENWRITE(7,600)NEE,R,S,detWRITE(7,*)(iv

30、en(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO30K=1,4DO20I=1,2DNX(I,K)=0.DO25M=1,2DNX(I,K尸DNX(I,K)+RJAC(I,M)*PN(M,K)25CONTINUE20CONTINUE30CONTINUE600FORMAT(1X,'ERR0R*NEGTIVEORZERO'* 'JACOBIANDETERMINANTFOR

31、'* 'ELEMENT'/'ELE.=',I5,'R=',F10.5,6X,'S=',F10.5,* 'det=',f12.5)RETURNENDC*SUBROUTINEFUN8(XYZ,R,S,DET)DIMENSIONXYZ(2,*),XI(4),ETA(4)COMMON/CMN5/FUN(4),PN(2,4),XJAC(2,2)DATAXI/-1.0,1.0,1.0,-1.0/DATAETA/-1.0,-1.0,1.0,1.0/DO10I=1,4G1=(1.0+XI(I)*R)G2=(1.0+ETA(I)*S)FUN(I)=0.25*G1*G2PN(1,I)=0.25*XI(I)*G2PN(2,I)=0.25*ETA(I)*G110CONTINUE文案标准DO80I=1,2DO75J=1,2DET=0.00DO70K=1,4DET=DET+PN(I,K)*XYZ(J,K)70CONTINUEXJAC(I,J)=