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

1、平面四边形四节点等参单元Fortran源程序平面四边形四节点等参单元Fortran源程序30/30平面四边形四节点等参单元Fortran源程序*C*C*FINITEELEMENTPROGRAM*C*FORTwoDIMENSIONALELASticityPROBLEM*C*WITH4NODE*PROGRAMELASTICITYcharacter*32dat,cchDIMENSIONSK(80000),COOR(2,300),AE(4,11),MEL(5,200),COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8

2、),NN(8)WRITE(*,*)PLEASEENTERINPUTFILENAMEREAD(*,(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)NUMBEROFMATERIALNM=,nmWRITE(7,(A,I6)NUMBEROFsurportingNC=,NrCALLINPUT(JR,COOR,AE,MEL)*C

3、ALLCBAND(MA,JR,MEL)DOI=1,NHSK(I)=0.0enddoCALLSK0(SK,MEL,COOR,JR,MA,AE)doI=1,NR(I)=0.0enddopauseaaastopREAD(4,*)NCP,NBE,izWRITE(*,(5i8)NCP,NBE,izWRITE(7,(5i8)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

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)WRITE(7,(A)PROGRAMSAFFHASBEENENDEDWRITE(*,(A)PROGRAMSAFFHASBEENENDEDSTOPcRETURNEND*SUBROUTINEINPUT(JR,COOR,AE,MEL)DIMENSIONJR(2,*),COOR(2,*),AE(4,*),

5、MEL(5,*)CONTINUEDO11J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME*11CONTINUEDO10I=1,NPDO10J=1,2JR(J,I)=1DO20I=1,NRREAD(4,*)IP,IX,IYCONTINUEN=0DO30I=1,NPDO30J=1,2IF(JR(J,I)30,30,25N=N+1JR(J,I)=NCONTINUEDO55J=1,NMREAD(4,*)JJ,(AE(I,JJ),I=1,4)55CONTINUEFORMAT(/20X,MATERIALPROPERTIES/(3X,I5,4(1

6、x,E8.3)RETURN*END*SUBROUTINECBAND(MA,JR,MEL)DIMENSIONMA(*),JR(2,*),MEL(5,*),NN(8)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHDO65I=1,N65MA(I)=0CONTINUECONTINUEL=NDO80I=1,2*4NNI=NN(I)IF(NNI.EQ.0)GOTO80IF(NNI.LT.L)L=NNI80CONTINUE*DO85M=1,2*4JP=NN(M)IF(JP.EQ.0)GOTO85JPL=JP-L+1IF(JPL.GT.MA(JP)MA(JP)=JPL85

7、CONTINUE90CONTINUEMX=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,(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*SUBROUTINESK

8、0(SK,MEL,COOR,JR,MA,AE)DIMENSIONSK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*),*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

9、(2)=0.00RSTG(3)=-RSTG(1)DO10IE=1,NENEE=IENME=MEL(5,IE)DO75K=1,4IEK=MEL(K,IE)iven(k)=IEKDO95M=1,2*JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)XYZ(M,K)=COOR(M,IEK)CONTINUECALLSTIF(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)60CONTINUE70CONTINUEwr

10、ite(7,1111)(ske(i,j),j=1,8),i=1,8)1111format(2x,8f12.2)10CONTINUERETURNENDC*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

11、,8SKE(I,J)=0.00CONTINUECONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(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+2*DO115J=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)CALLF

12、DNX(XYZ,DNX,DET,R,S,RJAC,iven,NEE)DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=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*SHCONTINUECONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3CONTINUECONTINUER

13、ETURNENDC*SUBROUTINECONCR(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)*CONTINUECONTINUERETURNENDC*SUBROUTINEBODYR(NBE,R,MEL,COOR,JR,AE)DIMENSIONR(*),MEL(5,*),COOR(2,*),JR(2,*),&AE(4,*),XYZ(2,4),iven(4)COMMON/C

14、MN1/NP,NE,NM,NRCOMMON/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)=i

15、ekDO95M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)XYZ(M,K)=COOR(M,IEK)CONTINUEDO99IS=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*GAMA30CONTINUE*CONTINUECONTINUECALLASLOAD(R)CONTINUERETURNENDC*SUBROUTINEFACER(iew,NSE,R,MEL,COOR,JR,WG)DIMENSI

16、ONR(*),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=1DO510I

17、E=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)XYZ(M,K)=COOR(M,IEK)CONTINUEIF(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,Z0,GAMA,

18、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)COMMON/GAUSS/R

19、STG(3),H(3)DATAKCRD/1,1,2,2/DATAKFACE/1,4,*2,3,*1,2,*4,3/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*GAMAENDDOENDIFML=KCRD(NSU)*IF(ML.EQ.

20、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)*FACTNUSCONTINUECONTINUE*RETURNEND*SUBROUTINEASLOAD(R)CONTINUERETURN*

21、SUBROUTINEDECOP(SK,MA)DIMENSIONSK(*),MA(*)COMMON/CMN2/N,MX,NHDO50I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1*IF(L1.GT.K)GOTO30DO20J=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)CONTINUECONTINUEIF(L.GT.K)GOTO50DO40LP=L,

22、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)CONTINUECONTINUE*RETURNEND*SUBROUTINEFOBA(SK,MA,R)DIMENSIONSK(*),MA(*),R(*)COMMON/CMN2/N,MX,NHDO10I=2,NL=I-MA(I)+MA(I-1)+1CONTINUECONTINUEDO20I=1,NII=MA(I)R(I)=R(I)/SK(II)CONTINUEDO30J1=2,NI=2+N-J1L=I-MA(I)+MA(I-

23、1)+1*K=I-1IF(L.GT.K)GOTO30DO25J=L,KIJ=MA(I)-I+JR(J)=R(J)-SK(IJ)*R(I)CONTINUECONTINUERETURNENDC*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)COM

24、MON/CMN5/FUN(4),PN(2,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)+MNN(JRR)=JR(M,IEK)CONTINUEE=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

25、,RR,SS,RJAC,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)=CI*B(1,J2)=0.B(2,J2)=CIB(3,J2)=BICONTINUEDO55II=1,3SIG(II)=0.00CONTINUEDO70K=1,8CONTINUECONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106C

26、ONTINUECALLOUTSTRE(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,*)(iven(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJ

27、AC(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)*CONTINUECONTINUECONTINUE600FORMAT(1X,ERR0R*NEGTIVEORZERO*JACOBIANDETERMINANTFOR*ELEMENT/ELE.=,I5,R=,F10.5,6X,S=,F10.5,det=,f12.5)RETURN*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.00DO