有限元程序设计

前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法。由于它的通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学与技术的快速开展，现已成为计算机辅助设计与计算机辅助制造的重要组成局部。由于有限元法是通过计算机实现的，因此它的计算机程序的研制和开发是其理论和方法应用于生产和科研实际的前提和根底。同时所研制和开发的计算机程序又是有限元理论的和方法的研究必要平台。本次程序设计是将Fortran程序设计与有限元理论结合。根据有限元理论知识，进行设计程序，从而获得简单平面问题的计算方法。平面4、8节点有限元公式及计算原理〔1〕通过Serendipity四边形单元格式构造插值函数。对于4节点单元，插值函数为：〔i=1,2,3,4〕对于8节点单元，插值函数为：〔i=1,2,3,4〕〔2〕通过，计算出D和B矩阵。〔3〕运用高斯积分公式求积分公式近似数值积分解：其中权函数〔4〕单元矩阵的变换：，因此所以面积微元可表示成。〔5〕单元刚度矩阵与载荷列向量可表示为：，将其代入整体刚度矩阵与载荷列向量可得整体刚度矩阵K和载荷列向量P。〔6〕模型泛函总势能为：。工程力学系根据最小势能原理可得：工程力学系运用上式即可求得求得各节点位移。程序结构〔1〕有限元程序系统的组成及分析过程前处理几何模型、材料参数、边界条件、分析类型定义、网格划分等。前处理几何模型、材料参数、边界条件、分析类型定义、网格划分等。核心：各种计算方法的实现。有限元分析本体程序核心：各种计算方法的实现。有限元分析本体程序以图形、曲线和表格等方式表达、分析计算结果。后处理以图形、曲线和表格等方式表达、分析计算结果。后处理〔2〕程序框图输入离散模型数据型数据计算单元刚度阵组集结构刚度矩阵计算单元等效结点载荷组集结构结点载荷列阵引入位移边界条件求解线性方程组其它辅助计算输出结果结束单元循环形成K形成P消除K的奇异求解Ka=P，得结点位移a计算应力、应变等〔3〕变量说明NUMNP：节点数；NUMEL：单元数；NUMMAT：材料种类；NLOAD：载荷数；NDM：每个节点坐标方向数；NDF：节点自由度；NEN：每个单元节点数。ID〔K，J〕：节点J的第K个自由度约束编号；X〔K，J〕：节点J的第K个坐标。IX〔K，J〕：J单元中K节点的全局节点编号。EE：杨氏弹性模量；XNU：泊松比；ITYPE：问题类型。F(K,J)：节点J在K方向上集中载荷大小。ALFL=TRUE：非对称矩阵集合；ALFL=FALSE：对称矩阵集合：Ｓ：单元刚度矩阵；Ｐ：载荷和内力向量；AD：对角元素；AU：上三角元素；AL：下三角元素；JP：上〔下〕三角元素每行〔列〕最后元素；LD：一个单元中各自由度方程数。JD：列高；ID：边界条件。〔4〕数据输入格式整数：I5格式；小数：F10.4格式。〔5〕程序本体C----------------------------------------CC.....FEA2DP---AFINITEELEMENTANALYSISPROGRAMFORC2DELASTICPROBLEMSCCTANGENTMATRIXISSTOREDWITHVARIOUDBANDMETHODCTHISPROGRAMISUSEDTODEMONSTRTETHEUSAGEOFVRIOUSBANDCSTORAGESCHEMOFSYMMETRICANDUNSYMMETRICTANGENTMATRIXCCChenYangCATCHONGQINGVNIVERSITY(15/06/2011)CC------------------------------------------------------PROGRAMFEA2DPCCA(1)-A(N1-1):X(NDM,NUMMNP);A(N1)-A(N2-1):F(NDF,NUMNP)CA(N2)-A(N3-1):B(NEQ);A(N3)-A(N4-1):AD(NEQ)CA(N4)-A(N5-1):AL(NAD);A(N5)-A(N6-1):NU(NAD)CCIA(1)-IA(N1-1):IX(NEN1,NUMEL);IA(N1)-IA(N2-1):ID(NDF,NUMNP)CIA(N2)-IA(N3-1):JD((NDF*NUMNP);IA(N3)-IA(N4-1):IDL(NEN*NUMEL*NDF)CIMPLICITREAL*8(A-H,O-Z) DIMENSIONA(100000),IA(1000) CHARACTER*80HEAD COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQ COMMON/SDATA/NDF,NDM,NEN1,NST COMMON/IOFILE/IOR,IOWCNMAXM=100000 IMAXM=1000 IOR=1 IOW=2CCOpenfilesfordatainputandoutputCOPEN(IOR,FILE='INPUT3.DAT',FORM='FORMATTED') OPEN(IOW,FILE='OUTPUT3.DAT')CC.....READTITLECREAD(IOR,'(A)')HEAD WRITE(IOW,'(A)')HEADCC.....READANDPRINTCONTROLINFORMATIONCCNUMNP:numberofnodesCNUMEL:numberofelementsCNUMMAT:numberofmaterialtypesCNLOAD:numberofloadsCNDM:numberofcoordinatsofeachnodeCNDF:numberofdegreesoffreedomCNEN:numberofnodesineachelementCREAD(IOR,'(7I5)')NUMNP,NUMEL,NUMMAT,NLOAD,NDM,NDF,NEN WRITE(IOW,2000)NUMNP,NUMEL,NUMMAT,NLOAD,NDM,NDF,NENCC.....SETPOITERSFORALLOCATIONOFDATAARRAYSCNEN1=NEN+4 NST=NEN*NDF NNEQ=NDF*NUMNPCN1=NDM*NUMNP+1 N2=N1+NDF*NUMNP+1CI1=NEN1*NUMEL+1 I2=I1+NDF*NUMNP+1 I3=I2+NDF*NUMNP+1 I4=I3+NUMEL*NEN*NDF+1CC.....CALLMESHINPUTSUBROUTINETOREADALLMESHDATACCALLPMESH(A(1),A(N1),IA(1),IA(I1),NDF,NDM,NEN1,NLOAD)CC.....COMPUTEPROFILECCALLPROFIL(IA(I2),IA(I3),IA(I1),IA(1),NDF,NEN1,NAD)CN3=N2+NEQ+1 N4=N3+NEQ+1N5=N4+NAD+1 N6=N5+NAD+1CCThelengthesofRealandIntegerarraysCWRITE(IOW,2222)N6,I4CCThelengthesofarrayexceedsthelimitationCIF(N6>NMAXM.OR.I4>IMAXM)THEN IF(N6>NMAXM)WRITE(IOW,3333)N6,NMAXM IF(I4>NMAXM)WRITE(IOW,4444)I4,IMAXM STOP ENDIFCC.....COMTUTEANDASEEMBLEELEMENTARRAYSCCALLASSEM(NAD,IA(1),IA(I1),IA(I2),A(1),A(N2),A(N3), 1A(N4),A(N5))CCFORMLOADVECTORCCALLPLOAD(IA(I1),A(N1),A(N2),NNEQ,NEQ)CC.....TRIANGULARDECOMPOSITIONOFAMATRIXSTOREDINPROFILEFORMCCALLDATRI(NDF,NUMNP,IA(I2),NEQ,NAD,.FALSE.,A(N3),A(N5),A(N5))CCFORUNSYMMTRICTANGENTMATIRXCCALLDATRI(NDF,NUMNP,IA(I2),NEQ,NAD,.TRUE.,A(N3),A(N4),A(N5))CCSOLVEEQUATIONSCCALLDASOL(NDF,NUMNP,A(N2),IA(I2),NEQ,NAD,AENGY,A(N3),A(N5),A(N5))CCFORUNSYMMETRICTANGENTMATRIXCCALLDASOL(NDF,NUMNP,A(N2),IA(I2),NEQ,NAD,AENGY,A(N3),A(N5),A(N5))CCOUTPUTNODALDISPLACEMENTSCCALLPRTDIS(IA(I1),A(N2),NDF,NUMNP,NEQ)CC.....CLOSEINPUTANDOUTPUTFILES;DESTROYTEMPORARYDISKFILESCCLOSE(IOR) CLOSE(IOW)CC.....INPUT/OUTPUTFORMATSC1000FORMAT(20A4)2000FORMAT(//X5X,'NUMBEROFNODALPOINTS=',I6/ 15X,'NUMBEROFELEMENTS=',I6/ 25X,'NUMBEROFMATERIALSETS=',I6/ 35X,'NUMBEROFNODALLOADS=',I6/45X,'DIMENSIONOFCOORDINATESPACE=',I6/ 55X,'DEGREEOFFREEDOMS/NODE=',I6/ 65X,'NODESPERELEMENT(MAXIMUM)=',I6)2222FORMAT(//,10X,'THELENGTHEOFREALARRAYIS',I10,/,110X,'THELENGTHEOFINTEGERARRAYIS',I10)3333FORMAT(//,10X,'THELENGTHEOFREALARRAY',I10,'EXCEEDTHE',1'MAXIMUNVALUE',I10)4444FORMAT(//,10X,'THELENGTHEOFINTEGERARRAY',I10,'EXCEEDTHE',1'MAXIMUNVALUE',I10)CSTOP ENDCCSUBROUTINEPMESH(X,F,IX,ID,NDF,NDM,NEN1,NLOAD)CC......DATAINPUTROUTINEFORMESHDESCRIPRIONCIMPLICITREAL*8(A-H,O-Z) DIMENSIONX(NDM,NUMNP),F(NDF,NUMNP),ID(NDF,NUMNP),IX(NEN1,NUMEL)COMMON/BDATA/HEAD(20) COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQ COMMON/MATER/EE,XNU,ITYPE COMMON/IOFILE/IOR,IOWCC.....INPUTCONSTRAINCODESANDNODALCOORDINATEDATACCID(K,J):constraincodeofKthdegreeoffreedomofnodeJ,=0:free,=1:fixedCX(K,J):KthcoordinateofnodeJCDOI=1,NUMNP READ(IOR,'(3I5,2F10.4)')J,(ID(K,J),K=1,NDM),(X(K,J),K=1,NDM) ENDDOCWRITE(IOW,'(//17HNODALCOORDINATES,/)') DOI=1,NUMNP WRITE(IOW,'(3I5,2F10.4)')I,(ID(K,I),K=1,NDM),(X(K,I),K=1,NDM) ENDDOCC.....ELEMENTDATAINPUTCCIX(K,J):globalnodenumberofKthnodeinelementJCDOI=1,NUMEL READ(IOR,'(9I5)')J,(IX(K,J),K=1,NEN) ENDDOCWRITE(IOW,'(//,18HELEMENTDEFINITION,/)') DOI=1,NUMEL WRITE(IOW,'(9I5)')J,(IX(K,J),K=1,NEN) ENDDOCC.....MATERIALDATAINPUTCCEE:Young'smodulus,XNU:PoissonratioCITYPE:Typeofproblem,=1,:Planestress,=2:Planestrain,=3:Axi-symmetricCREAD(IOR,'(2F10.4,I5)')EE,XNU,ITYPE WRITE(IOW,'(//,19HMATEIALPROPERTIES,/)') WRITE(IOW,'(2(E10.4,5X),I5)')EE,XNU,ITYPECC.....FORCE/DISPDATAINPUTCCF(K,J):ConcentrateloadatnodeJinKdirectionCF=0.0D0 DOI=1,NLOAD READ(IOR,'(I5,2F10.4)')J,(F(K,J),K=1,NDF) ENDDOCWRITE(IOW,'(//,20HAPPLIEDNODALFORCES,/)') DOI=1,NLOAD WRITE(IOW,'(I5,2F10.4)')J,(F(K,J),K=1,NDF) ENDDOCRETURNCCFORMATSTATEMENTSC2000FORMAT('Mesh1>',$)3000FORMAT(1X,'**WARNING**ELEMENTCONNECTIONSNECESSARY'1'TOUSEBLOCKINMACROPROGRAM')4000FORMAT('**CURRENTPROBLEMVALIES**'/I6,'NODES,',1I5,'ELMTS,',I3,'MATLS,',I2,'DIMS,',I2,'DOF/NODE,', 2I3,'NODES/ELMT') ENDCCSUBROUTINEASSEM(NAD,IX,ID,JD,X,B,AD,AL,AU)CCCALLelementsubroutineandassembleglobaltangentmatrixCIMPLICITREAL*8(A-H,O-Z) DIMENSIONILX(NEN),XL(NDF,NEN),LD(NDF,NEN),S(NST,NST),P(NST)DIMENSIONIX(NEN1,NUMEL),ID(NDF,NUMNP),JD(NDF*NUMNP)DIMENSIONX(NDM,NUMNP),B(NEQ),AD(NEQ),AL(NAD),AU(NAD) COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQ COMMON/SDATA/NDF,NDM,NEN1,NSTCNEL=NENCCElenmentloopCDO320N=1,NUMEL S=0.0D0!elementstiffnessmatrix P=0.0D0!nodalforce NE=N DO310I=1,NEN ILX(I)=IX(I,NE)!currentelementdefinition DOK=1,NDM XL(K,I)=X(K,ILX(I))!nodalcoordsincurrentelement ENDDO KK=ILX(I) DOK=1,NDF LD(K,I)=ID(K,KK)!equationnumbers ENDDO310CONTINUECCCallelementlibCCALLELMT01(XL,ILX,S,P,NDF,NDM,NST)CCAsemmbletangentmatrixandloadvectorifneededC CALLDASBLY(NDF,NAD,S,P,LD,JD,NST,B,AD,AL,AU)C320CONTINUE!endelementloopCRETURN ENDCCSUBROUTINEDASBLY(NDF,NAD,S,P,LD,JP,NS,B,AD,AL,AU)CC.....ASSEMBLETHESYMMETRICORUNSYMMETRICARRAYSFOR'DASOL'CIMPLICITREAL*8(A-H,O-Z)CLOGICALALFL,AUFL,BFLDIMENSIONAD(NEQ),AL(NAD),AU(NAD)DIMENSIONLD(NS),JP(NDF*NUMNP),B(NEQ),S(NS,NS),P(NS) COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQ COMMON/IOFILE/IOR,IOWCCALFL=TRUE:FORUNSYMMETRICMATIRXASSEMBLECALFL=FALSE:FORSYMMETRICMATIRXASSEMBLECS:ELEMENTSTIFFNESSMATRIXCP:LOADORINTERNALFORCEVECTORCAD:DIAGONALELEMENTSCAU:UPPERTRIANGLEELEMENTSCAL:LOWERTRIANGLEELEMENTSCJP:POINTERTOLASTELEMENTINEACHROW/COLUMNOFAL/AURESPECTIVECLD:EQUATIONNUMBERSOFEACHFREEDOMDEGREEINANELEMENT(GETFROMID)CC.....LOOPTHROUGHTHEROWSTOPERFORMTHEASSEMBLYCDO200I=1,NS II=LD(I) IF(II.GT.0)THENC IF(AUFL)THEN!AssemblestiffnessmatrixCC.....LOOPTHROUGHTHECOLUMNSTOPERFORMTHEASSEMBLYCDO100J=1,NS IF(LD(J).EQ.II)THEN AD(II)=AD(II)+S(I,J) ELSEIF(LD(J).GT.II)THEN JC=LD(J) JJ=II+JP(JC)-JC+1AU(JJ)=AU(JJ)+S(I,J)CIF(ALFL)AL(JJ)=AL(JJ)+S(J,I)!UnsymmetricENDIF100CONTINUE ENDIFCIF(BFL)B(II)=B(II)+P(I)!AssemblenodalforceCENDIF200CONTINUECRETURN ENDCCSUBROUTINEDASOL(NDF,NUMNP,B,JP,NEQ,NAD,ENERGY,AD,AL,AU)CC.....SOLUTIONOFSYMMETRICEQUATIONSINPROFILEFORMC.....COEFICIENTMATRIXMUSTBEDECOMPOSEDINTOITSTRIANGULARC.....FACTORUSINGDATRIBEFORCEUSINGDASOL.CCJP:POINTERTOLASTELEMENTINEACHROW/COLUMNOFAL/AURESPECIVECCIMPLICITREAL*8(A-H,O-Z) DIMENSIONAD(NEQ),AL(NAD),AU(NAD) DIMENSIONB(NEQ),JP(NDF*NUMNP) COMMON/IOFILE/IOR,IOW DATAZERO/0.0D0/CC.....FINDTHEFIRSTNONZEROENTRYINTHERINGHANDSIDECDOIS=1,NEQ IF(B(IS).NE.ZERO)GOTO200 ENDDO WRITE(IOW,2000) RETURNC200IF(IS.LT.NEQ)THENCC.....REDUCETHERIGHTHANDSIDECDO300J=IS+1,NEQ JR=JP(J-1) JH=JP(J)-JR IF(JH.GT.0)THEN B(J)=B(J)-DOT(AL(JR+1),B(J-JH),JH) ENDIF300CONTINUEENDIFCC.....MULTIPLYINVERSEOFDIAGONALELEMENTSCENERGY=ZERO DO400J=IS,NEQ BD=B(J) B(J)=B(J)*AD(J) ENERGY=ENERGY+BD*B(J)400CONTINUECC.....BACKSUBSTITUTIONCIF(NEQ.GT.1)THEN DO500J=NEQ,2,-1 JR=JP(J-1) JH=JP(J)-JR IF(JH.GT.0)THEN CALLSAXPB(AU(JR+1),B(J-JH),-B(J),JH,B(J-JH)) ENDIF500CONTINUE ENDIFCRETURNC2000FORMAT('**DASOLWARNING1**ZERORIGHT-HAND-SIDEVECTOR') ENDCCSUBROUTINEDATEST(AU,JH,DAVAL)CC.....TESTFORRANKCIMPLICITREAL*8(A-H,O-Z) DIMENSIONAU(JH)CDAVAL=0.0D0CDOJ=1,JH DAVAL=DAVAL+ABS(AU(J)) ENDDOCRETURN ENDCCSUBROUTINEDATRI(NDF,NUMNP,JP,NEQ,NAD,FLG,AD,AL,AU)CC.....TRIANGULARDECOMPOSIONTIONOFAMATRIXSTOREDINPROFILEFORMCIMPLICITREAL*8(A-H,O-Z) LOGICALFLG DIMENSIONJP(NDF*NUMNP),AD(NEQ),AL(NAD),AU(NAD) COMMON/IOFILE/IOR,IOWCC.....N.B.TOLSHOULDBESETTOAPPROXIMATEHALF-WORDPRECISION.CDATAZERO,ONE/0.0D0,1.0D0/,TOL/0.5D-07/CC.....SETINITIALVALUESFORCONTDITIONINGCHECKCDIMX=ZERO DIMN=ZEROCDOJ=1,NEQ DIMN=MAX(DIMN,ABS(AD(J))) ENDDO DFIG=ZEROCC.....LOOPTHROUGHTHECOLUMNSTOPERFORMTHETRIANGULARDECOMPOSITIONCJD=1 DO200J=1,NEQ JR=JD+1 JD=JP(J) JH=JD-JR IF(JH.GT.0)THEN IS=J-JH IE=J-1CC.....IFDIAGONALISZEORCOMPUTEANORMFORSINGULARITYTESTCIF(AD(J).EQ.ZERO)CALLDATEST(AU(JR),JH,DAVAL) DO100I=IS,IE JR=JR+1 ID=JP(I) IH=MIN(ID-JP(I-1),I-IS+1) IF(IH.GT.0)THEN JRH=JR-IH IDH=ID-IH+1 AU(JR)=AU(JR)-DOT(AU(JRH),AL(IDH),IH) IF(FLG)AL(JR)=AL(JR)-DOT(AL(JRH),AU(IDH),IH) ENDIF100CONTINUE ENDIFCC.....REDUCETHEDIAGONALCIF(JH.GE.0)THEN DD=AD(J) JR=JD-JH JRH=J-JH-1 CALLDREDU(AL(JR),AU(JR),AD(JRH),JH+1,FLG,AD(J))CC.....CHECKFORPOSSIBLEERRORSANDPRINTWARNINGSCIF(ABS(AD(J)).LT.TOL*ABS(DD))WRITE(IOW,2000)J IF(DD.LT.ZERO.AND.AD(J).GT.ZERO)WRITE(IOW,2001)J IF(DD.GT.ZERO.AND.AD(J).LT.ZERO)WRITE(IOW,2001)J IF(AD(J).EQ.ZERO)WRITE(IOW,2002)J IF(DD.EQ.ZERO.AND.JH.GT.0)THEN IF(ABS(AD(J)).LT.TOL*DAVAL)WRITE(IOW,2003)J ENDIF ENDIFCC.....STROERECIPROCALOFDIAGONAL,COMPUTECONDITIONCHECKSCIF(AD(J).NE.ZERO)THEN DIMX=MAX(DIMX,ABS(AD(J))) DIMN=MIN(DIMN,ABS(AD(J))) DFIG=MAX(DFIG,ABS(DD/AD(J))) AD(J)=ONE/AD(J) ENDIF200CONTINUECC.....PRINTCONDITIONINGINFORMATIONCDD=ZERO IF(DIMN.NE.ZERO)DD=DIMX/DIMN IFIG=DLOG10(DFIG)+0.6 WRITE(IOW,2004)DIMX,DIMN,DD,IFIGCRETURNCC.....FORMATSC2000FORMAT('**DATRIWARNING1**LOSSOFATLEAST7DIGITSIN',1'REDUCINGDIAGONALOFEQUATION',I5)2001FORMAT('**DATRIWARNING2**SIGNOFCHANGEDWHEN',1'REDUCINGEQUATION',I5)2002FORMAT('**DATRIWARNING3**REDUCEDDIAGONALISZEROZERIFOR',1'EQUATION',I5)2003FORMAT('**DATRIWARNING4**RANKFAILUREFFOZEROUNREDUCED',1'DIAGONALINEQUATION',I5)2004FORMAT(//'CONDITONCHECK:D-MAX',E11.4,';D-MIN',E11.4,1';RATIO',E11.4/'MAXIMIMNO.DIAGONALDIGITSLOST:',I3)2005FORMAT('CONDCK:DMAX',1P1E9.2,';DMIN',1P1E9.2,1';RATIO',1P1E9.2)ENDCCSUBROUTINEDREDU(AL,AU,AD,JH,FLG,DJ)CC.....REDUCEDIAGONALELEMENTINTRIANGULARDECOMPOSITIONCIMPLICITREAL*8(A-H,O-Z) LOGICALFLG DIMENSIONAL(JH),AU(JH),AD(JH)CDOJ=1,JH UD=AU(J)*AD(J) DJ=DJ-AL(J)*UD AU(J)=UD ENDDOCC.....FINISHCOMPUTATIONOFCOLUMNOFALFORUNSYMMETRICMATRICESCIF(FLG)THEN DOJ=1,JH AL(J)=AL(J)*AD(J) ENDDO ENDIFCRETURN ENDCCSUBROUTINEPROFIL(JD,IDL,ID,IX,NDF,NEN1,NAD)CC.....COMPUTEPROFILEOFGLOBALARRAYSCIMPLICITREAL*8(A-H,O-Z) DIMENSIONJD(NDF*NUMNP),IDL(NUMEL*NEN*NDF),ID(NDF,NUMNP),1IX(NEN1,NUMEL)COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQ COMMON/FRDATA/MAXF COMMON/IOFILE/IOR,IOWCCJD:COLUMNHIGHT(ADDRESSOFDIAGONALELEMENTS)CID:BOUDARYCONDITIONCODESBEFORETHISBUBROUTINE'SRUNNING CID:EQUATIONNUMBERSINGLOBALARRAY(EXCLUDEDRESTRAINEDNODES)AFTERRUNNINGCIDL:ELEMENTSTRECHORDERCNAD:TOTALNUMBEROFNON-ZEROELEMENTSEXCEPTDIAGONALELEMENTSCINGLOBALTANGENTMATRIXCC.....SETUPTHEEQUATIONNUMBERSCNEQ=0CDO10K=1,NUMNP DO10N=1,NDF J=ID(N,K) IF(J.EQ.0)THEN NEQ=NEQ+1 ID(N,K)=NEQ ELSE ID(N,K)=0 ENDIF10CONTINUECC.....COMPUTECOLUMNHEIGHTSCCALLPCONSI(JD,NEQ,0)CDO50N=1,NUMEL MM=0 NAD=0 DO30I=1,NEN II=IABS(IX(I,N)) IF(II.GT.0)THEN DO20J=1,NDF JJ=ID(J,II) IF(JJ.GT.0)THEN IF(MM.EQ.0)MM=JJ MM=MIN(MM,JJ) NAD=NAD+1 IDL(NAD)=JJ ENDIF20CONTINUEENDIF30CONTINUEIF(NAD.GT.0)THEN DO40I=1,NAD II=IDL(I) JJ=JD(II) JD(II)=MAX(JJ,II-MM)40 CONTINUEENDIF50CONTINUECC.....COMPUTEDIAGONGALPOINTERSFORPROFILECNAD=0 JD(1)=0 IF(NEQ.GT.1)THEN DO60N=2,NEQ JD(N)=JD(N)+JD(N-1)60CONTINUENAD=JD(NEQ) ENDIFCC.....SETELEMENTSEARCHORDERTOSEQUENTIALCDO70N=1,NUMEL IDL(N)=N70CONTINUECC.....EQUATIONSUMMARYCMAXF=0 MM=0 IF(NEQ.GT.0)MM=(NAD+NEQ)/NEQ WRITE(IOW,2001)NEQ,NUMNP,MM,NUMEL,NAD,NUMMATCRETURNC2001FORMAT(5X,'NEQ=',I5,5X,'NUMNP=',I5,5X,'MM=',I5,/5X,1'NUMEL=',I5,5X,'NAD=',I5,5X,'NUMMAT=',I5/) ENDCCSUBROUTINESAXPB(A,B,X,N,C)CC.....VECTORTIMESSCALARADDEDTOSECONDVECTORCIMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N),B(N),C(N)CDOK=1,N C(K)=A(K)*X+B(K) ENDDOCRETURN ENDCCSUBROUTINEPCONSI(IV,NN,IC)CC.....ZEROINTEGERARRAYCDIMENSIONIV(NN)CDON=1,NN IV(N)=IC ENDDOCRETURN ENDCCSUBROUTINEELMT01(XL,ILX,S,P,NDF,NDM,NST)CC.....PLANELINEARELASTICELEMENTROUTINECityp=1:planestressC=2:planestrainC=3:axisymmetricCIMPLICITREAL*8(A-H,O-Z) DIMENSIONXL(NDM,NEN),ILX(NEN),SIGR(6) DIMENSIOND(18),S(NST,NST),P(NST),SHP(3,9),SG(16),TG(16),WG(16) CHARACTERWD(3)*12 COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/MATER/EE,XNU,ITYPE COMMON/IOFILE/IOR,IOW DATAWD/'PLANESTRESS','PLANESTRAIN','AXISYMMETRIC'/CCXL(NDM,NEN):COORDSOFEACHNODEINCURRENTELEMENTCILX(NEN):ELEMENTDEFINITIONOFCURRENTELEMENTCD(18):MATERIALSPROPERTIESCS(NST,NST):ELEMENTSTIFFNESSMATRIXCp(NS):NODALFORCEANDINTERNALFORCECSHP(3,9):SHAPEFUNCTIONANDITSDERIVATIVESCSG(16),TG(16),WG(16):WEIGHTCOEFFICIENTSOFGUASSINTERGTATIONCL,K:integrationpointsCL=2 K=2 E=EE NEL=NENCCD(14):thickness;D(11),D(12):bodyforcesC.....SETMATERIALPATAMETERTYPEANDFLAGSCITYP=MAX(1,MIN(ITYP,3)) J=MIN(ITYP,2)CD(1)=E*(1.+(1-J)*XNU)/(1.+XNU)/(1.-J*XNU) D(2)=XNU*D(1)/(1.+(1-J)*XNU) D(3)=E/2./(1.+XNU) D(13)=D(2)*(J-1) IF((D(14).LE.0.0D0).OR.ITYP.GE.2)D(14)=1.0 D(15)=ITYP D(16)=E D(17)=XNU D(18)=-XNU/E L=MIN(4,MAX(1,L)) K=MIN(4,MAX(1,K)) D(5)=L D(6)=KCD(9)=T0CD(10)=E*ALP/(1.-J*XNU)LINT=0CWRITE(IOW,2000)WD(ITYP),D(16),D(17),D(4),L,K,D(14), 1D(11),D(12)CC.....STIFFNESS/RESIDUALCOMPUTATIONCL=KCCcomputecordinatesandweightsofintegtationpointC`SG,TG:Cootds;WG=Wp*WqCIF(L*L.NE.LINT)CALLPGUASS(L,LINT,SG,TG,WG)CC.....COMPUTEINTEGRALSOFSHAPEFUNCTIONSCDO340L=1,LINTCCcomputeshapefunctionandtheirderivativestolocalCandglobalcoordinatesystemCCALLSHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEN,ILX,.FALSE.)CCcomputeglobalcoordinatesofintegrationpointsCXX=0.0 YY=0.0 DOJ=1,NEN XX=XX+SHP(3,J)*XL(1,J)YY=YY+SHP(3,J)*XL(2,J) ENDDO XSJ=XSJ*WG(L)*D(14)!XSJ+|J|(Sp,Tq)*Wp*Wq*tCC.....COMPUTEJACOBIANCORRECTIONCforplanestressandstrainproblemsCIF(ITYP.LE.2)THEN DV=XSJ XSJ=0.0 ZZ=0.0CSIGR4=-D(11)*DV!D(11)BODYFORCEELSECCforanisymmetricproblemCDV=XSJ*XX*3.1415926*2. ZZ=1./XXC SIGR4=SIGR(4)*XSJ-D(11)*DV ENDIF J1=1CC.....LOOPOVERROWSCDO330J=1,NEL W11=SHP(1,J)*DV W12=SHP(2,J)*DV W22=SHP(3,J)*XSJ W22=SHP(3,J)*DV*ZZCCC.....COMPUTETHEINTERNALFORCESCoutofbalanceCCP(J1)=P(J1)-(SHP(1,J)*SIGR(1)+SHP(2,J)*SIGR(2))*DVC1-SHP(3,J)*SIGR4CP(J1+1)=P(J1+1)-(SHP(1,J)*SIGR(2)+SHP(2,J)*SIGR(3))*DVC1+D(12)*SHP(3,J)*DV!D(12)BODYFORCECC.....LOOPOVERCOLUMNS(SYMMETRYNOTED)CcomputestiffnessmatrixCK1=J1 A11=D(1)*W11+D(2)*W22 A21=D(2)*W11+D(1)*W22 A31=D(2)*(W11+W22) A41=D(3)*W12 A12=D(2)*W12 A32=D(1)*W12 A42=D(3)*W11 DO320K=J,NEL W11=SHP(1,K) W12=SHP(2,K) W22=SHP(3,K)*ZZ S(J1,K1)=S(J1,K1)+W11*A11+W22*A21+W12*A41 S(J1+1,K1)=S(J1+1,K1)+(W11+W22)*A12+W12*A42 S(J1,K1+1)=S(J1,K1+1)+W12*A31+W11*A41 S(J1+1,K1+1)=S(J1+1,K1+1)+W12*A32+W11*A42 K1=K1+NDF320CONTINUEJ1=NDF+J1330CONTINUE340CONTINUECC.....MAKESTIFFNESSSYMMETRICCDO360J=1,NST DO360K=J,NST S(K,J)=S(J,K)360CONTINUECRETURNCC.....FORMATSFORINPUT-OUTPUTC1000FORMAT(3F10.0,3I10)1001FORMAT(8F10.0)2000FORMAT(/5X,A12,'LINEARELASTICELEMENT'//110X,'MODULUS',E18.5/10X,'POISSIONRATIO',F8.5/10X,'DENSITY',E18.5/210X,'GUASSPTR/DIR',I3/10X,'STRESSPTS',I6/10X,'THICKNESS',E16.5/310X,'1-GRAVITY',E16.5/10X,'2-GTAVITY',E16.5/10X,'ALPHA',E20.5/410X,'BASETEMP',E16.5/)2001FORMAT(5X,'ELEMENTSTRESSES'//'ELMT1-COORD',2X,'11-STRESS',2X,1'12-STRESS',2X,'22-STRESS',2X,'33-STRESS',3X,'1-COORD',2X,3X,2'2-STRESS'/'MATL2-COORD',2X,'11-STRAIN',2X,'12-STRAIN'2X,3'22-STRAIN',2X,'33-STRAIN',6X,'ANGLE'/39('-'))2002FORMAT(I4,0P1F9.3,1P6E11.3/I4,0P1F9.3,1P4E11.3,0P1F11.2/)5000FORMAT('INPUT:E,NU,RHO,PTS/STIFF,PTS/STRE',1',TYPE(1=STRESS,2=STRAIN,3=AXISM)',/3X,'>',$)5001FORMAT('INPUT:THICKNESS,1-BODYFORCE,1-BODYFORCE,ALPHA,'1,'TEMP-BASE'/3X,'>',$) ENDCCSUBROUTINESHAPE(SS,TT,XL,SHP,XSJ,NDM,NEL,ILX,FLG)CC.....SHAPEFUNCTIONROUTINEFORTWODIMENSIONELEMENTSCIMPLICITREAL*8(A-H,O-Z) LOGICALFLG DIMENSIONXL(NDM,NEL),S(4),T(4),X(NEL) DIMENSIONSHP(3,NEL),XS(2,2),SX(2,2) DATAS/-0.5D0,0.5D0,0.5D0,-0.5D0/,1T/-0.5D0,-0.5D0,0.5D0,0.5D0/CC.....FORM4-NODEQUATRILATERALSHAPEFUNCTIONSCCNEL:nuberofnodesperelementCDO100I=1,4 SHP(3,I)=(0.5+S(I)*SS)*(0.5+T(I)*TT) SHP(1,I)=S(I)*(0.5+T(I)*TT) SHP(2,I)=T(I)*(0.5+S(I)*SS)100CONTINUECC.....FORMTRIANGGEBUADDINGTHEIRANDFOURTHTOGETHERCfortriangleelementCIF(NEL.EQ.3)THEN DOI=1,3 SHP(I,3)=SHP(I,3)+SHP(I,4) ENDDO ENDIFCC.....ADDQUATRATICTERMSIFNECESSARYCforelementwithmorethan4nodesCIF(NEL.GT.4)CALLSHAP2(SS,TT,SHP,ILX,NEL)CC.....CONSTRUCTJACOBIANANDITSINVERSECDO125I=1,2DO125J=1,2 XS(I,J)=0.0 DO120K=1,NEL XS(I,J)=XS(I,J)+XL(I,K)*SHP(J,K)120CONTINUE125CONTINUECCXSJ:determinateofJacobmatrixCXSJ=XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1)CIF(FLG)RETURNCFLG=false:formglobalderivativesCIF(XSJ.LE.0.0D0)XSJ=1.0 SX(1,1)=XS(2,2)/XSJ SX(2,2)=XS(1,1)/XSJSX(1,2)=-XS(1,2)/XSJSX(2,1)=-XS(2,1)/XSJCC....FORMGLOBALDERIVATIVESCDO130I=1,NEL TP=SHP(1,I)*SX(1,1)+SHP(2,I)*SX(2,1) SHP(2,I)=SHP(1,I)*SX(1,2)+SHP(2,I)*SX(2,2) SHP(1,I)=TP130CONTINUECRETURN ENDCCSUBROUTINESHAP2(S,T,SHP,ILX,NEL)CC....ADDQUADTATICFUNCTIONASNECESSARYCIMPLICITREAL*8(A-H,O-Z) DIMENSIONSHP(3,9),ILX(NEL)CS2=(1.-S*S)/2. T2=(1.-T*T)/2. DO100I=5,9 DO100J=1,3 SHP(J,I)=0.0100CONTINUECC.....MIDSIZENODES(SERENIPITY)CIF(ILX(5).EQ.0)GOTO101 SHP(1,5)=-S*(1.-T) SHP(2,5)=-S2 SHP(3,5)=S2*(1.-T)101IF(NEL.LT.6)GOTO107IF(ILX(6).EQ.0)GOTO102 SHP(1,6)=T2 SHP(2,6)=-T*(1.+S) SHP(3,6)=T2*(1.+S)102IF(NEL.LT.7)GOTO107IF(ILX(7).EQ.0)GOTO103 SHP(1,7)=-S*(1.+T) SHP(2,7)=S2 SHP(3,7)=S2*(1.+T)103IF(NEL.LT.8)GOTO107IF(ILX(8).EQ.0)GOTO104 SHP(1,8)=-T2 SHP(2,8)=-T*(1.-S) SHP(3,8)=T2*(1.-S)CC.....INTERIORNODE(LAGRAGIAN)C104IF(NEL.LT.9)GOTO107IF(ILX(9).EQ.0)GOTO107SHP(1,9)=-4.*S*T2 SHP(2,9)=-4.*T*S2 SHP(3,9)=4.*S2*T2CC.....CORRECTEDGENODESFORINTERIORNODE(LAGRANGIAN)CDO106J=1,3 DO105I=1,4105SHP(J,I)=SHP(J,I)-0.25*SHP(J,9)DO106I=5,8106IF(ILX(I).NE.0)SHP(J,I)=SHP(J,I)-.5*SHP(J,9)CC.....CORRECTCORNERNODESFORPRESENSEOFMIDSIZENODESC107DO108I=1,4K=MOD(I+2,4)+5 L=I+4 DO108J=1,3108SHP(J,I)=SHP(J,I)-0.5*(SHP(J,K)+SHP(J,L))RETURN ENDCCSUBROUTINEPGUASS(L,LINT,R,Z,W)CC.....GUASSPOINTSANDWEIGHTSFORTWODIMENSIONSCIMPLICITREAL*8(A-H,O-Z) DIMENSIONLR(9),LZ(9),LW(9),R(16),Z(16),W(16)CCOMMON/ELDTAT/DM,N,MA,MCT,IEL,NELDATALR/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/ DATALW/4*25,4*40,64/CCLint:NumberofintegrationpointsCR,Z:CoordinatesofIntegrationpointsCW:Wp*Wq,productofthetwoweightsCLINT=L*L GOTO(1,2,3),LCC.....1X1INTEGERATIONC1R(1)=0.Z(1)=0. W(1)=4.CRETURNCC.....2X2INTEGERATIONC2G=1.0/SQRT(3.D0)DOI=1,4 R(I)=G*LR(I) Z(I)=G*LZ(I) W(I)=1. ENDDOCRETURNCC.....3X3INTEGERATIONC3G=SQRT(0.60D0)H=1.0/81.0D0CDOI=1,9 R(I)=G*LR(I) Z(I)=G*LZ(I) W(I)=H*LW(I) ENDDOCRETURNCENDCCSUBROUTINEPLOAD(ID,F,B,NNEQ,NEQ)CC.....FORMLOADVECTORINCOMPACTFORMCIMPLICITREAL*8(A-H,O-Z) DIMENSIONF(NNEQ),B(NEQ),ID(NNEQ) COMMON/IOFILE/IOR,IOWCB=0.0D0CDON=1,NNEQ J=ID(N) IF(J.GT.0)THEN B(J)=F(N) ENDIF ENDDOCRETURNENDCCSUBROUTINEPRTDIS(ID,B,NDF,NUMNP,NEQ)CCPrintoutnodaldisplacementsCIMPLICITREAL*8(A-H,O-Z) DIMENSIONID(NDF,NUMNP),B(NEQ),U(NDF,NUMNP) COMMON/IOFILE/IOR,IOWCU=0.0D0 DO100I=1,NUMNP DOJ=1,NDF N=ID(J,I) IF(N>0)U(J,I)=B(N) ENDDO100CONTINUECCOutnodaldisplacementsCWRITE(IOW,'(//,19HNODALDISPLACEMENTS,/)') DOI=1,NUMNP WRITE(IOW,'(5X,I5,2X,3(E12.4,3X))')I,(U(K,I),K=1,NDF) ENDDOCRETURN ENDCCDOUBLEPRECISIONFUNCTIONDOT(A,B,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N),B(N)CC.....DOTPRODUCTFUNCTIONCDOT=0.0D0 DO10K=1,N DOT=DOT+A(K)*B(K)10CONTINUERETURN END四、算例及分析算例1：悬臂梁的有限元解；〔1〕将悬臂梁划分为20单元四节点形式，计算其数值解。网格划分及受力情况如下列图：程序input文件中输入：example14220112241110.00000.00002000.05000.00003000.05000.20004110.00000.20005000.10000.00006000.10000.20007000.15000.00008000.15000.20009000.20000.000010000.20000.200011000.25000.000012000.25000.200013000.30000.000014000.30000.200015000.35000.000016000.35000.200017000.40000.000018000.40000.200019000.45000.000020000.45000.200021000.50000.000022000.50000.200023000.55000.000024000.55000.200025000.60000.000026000.60000.200027000.65000.000028000.65000.200029000.70000.000030000.70000.200031000.75000.000032000.75000.200033000.80000.000034000.80000.200035000.85000.000036000.85000.200037000.90000.000038000.90000.200039000.95000.000040000.95000.200041001.00000.000042001.00000.20001123422563357864791085911121061113141271315161481517181691719201810192122201121232422122325262413252728261427293028152931323016313334321733353634183537383619373940382039414240200.0E+90.30002调试程序得到结果为：example1NODALDISPLACEMENTS10.0000E+000.0000E+002-0.3256E-07-0.1139E-0730.3256E-07-0.1139E-0740.0000E+000.0000E+005-0.6345E-07-0.3864E-0760.6345E-07-0.3864E-077-0.9267E-07-0.8092E-0780.9267E-07-0.8092E-079-0.1202E-