有限元程序設(shè)計(jì)

前言有限單元法是在當(dāng)今工程分析中獲得最廣泛應(yīng)用的數(shù)值計(jì)算方法。由于它的通用性和有效性，受到工程技術(shù)界的高度重視。伴隨著計(jì)算機(jī)科學(xué)與技術(shù)的快速開展，現(xiàn)已成為計(jì)算機(jī)輔助設(shè)計(jì)與計(jì)算機(jī)輔助制造的重要組成局部。由于有限元法是通過(guò)計(jì)算機(jī)實(shí)現(xiàn)的，因此它的計(jì)算機(jī)程序的研制和開發(fā)是其理論和方法應(yīng)用于生產(chǎn)和科研實(shí)際的前提和根底。同時(shí)所研制和開發(fā)的計(jì)算機(jī)程序又是有限元理論的和方法的研究必要平臺(tái)。本次程序設(shè)計(jì)是將Fortran程序設(shè)計(jì)與有限元理論結(jié)合。根據(jù)有限元理論知識(shí)，進(jìn)行設(shè)計(jì)程序，從而獲得簡(jiǎn)單平面問(wèn)題的計(jì)算方法。平面4、8節(jié)點(diǎn)有限元公式及計(jì)算原理〔1〕通過(guò)Serendipity四邊形單元格式構(gòu)造插值函數(shù)。對(duì)于4節(jié)點(diǎn)單元，插值函數(shù)為：〔i=1,2,3,4〕對(duì)于8節(jié)點(diǎn)單元，插值函數(shù)為：〔i=1,2,3,4〕〔2〕通過(guò)，計(jì)算出D和B矩陣?！?〕運(yùn)用高斯積分公式求積分公式近似數(shù)值積分解：其中權(quán)函數(shù)〔4〕單元矩陣的變換：，因此所以面積微元可表示成?！?〕單元?jiǎng)偠染仃嚺c載荷列向量可表示為：，將其代入整體剛度矩陣與載荷列向量可得整體剛度矩陣K和載荷列向量P?！?〕模型泛函總勢(shì)能為：。工程力學(xué)系根據(jù)最小勢(shì)能原理可得：工程力學(xué)系運(yùn)用上式即可求得求得各節(jié)點(diǎn)位移。程序結(jié)構(gòu)〔1〕有限元程序系統(tǒng)的組成及分析過(guò)程前處理幾何模型、材料參數(shù)、邊界條件、分析類型定義、網(wǎng)格劃分等。前處理幾何模型、材料參數(shù)、邊界條件、分析類型定義、網(wǎng)格劃分等。核心：各種計(jì)算方法的實(shí)現(xiàn)。有限元分析本體程序核心：各種計(jì)算方法的實(shí)現(xiàn)。有限元分析本體程序以圖形、曲線和表格等方式表達(dá)、分析計(jì)算結(jié)果。后處理以圖形、曲線和表格等方式表達(dá)、分析計(jì)算結(jié)果。后處理〔2〕程序框圖輸入離散模型數(shù)據(jù)型數(shù)據(jù)計(jì)算單元?jiǎng)偠汝嚱M集結(jié)構(gòu)剛度矩陣計(jì)算單元等效結(jié)點(diǎn)載荷組集結(jié)構(gòu)結(jié)點(diǎn)載荷列陣引入位移邊界條件求解線性方程組其它輔助計(jì)算輸出結(jié)果結(jié)束單元循環(huán)形成K形成P消除K的奇異求解Ka=P，得結(jié)點(diǎn)位移a計(jì)算應(yīng)力、應(yīng)變等〔3〕變量說(shuō)明NUMNP：節(jié)點(diǎn)數(shù)；NUMEL：?jiǎn)卧獢?shù)；NUMMAT：材料種類；NLOAD：載荷數(shù)；NDM：每個(gè)節(jié)點(diǎn)坐標(biāo)方向數(shù)；NDF：節(jié)點(diǎn)自由度；NEN：每個(gè)單元節(jié)點(diǎn)數(shù)。ID〔K，J〕：節(jié)點(diǎn)J的第K個(gè)自由度約束編號(hào)；X〔K，J〕：節(jié)點(diǎn)J的第K個(gè)坐標(biāo)。IX〔K，J〕：J單元中K節(jié)點(diǎn)的全局節(jié)點(diǎn)編號(hào)。EE：楊氏彈性模量；XNU：泊松比；ITYPE：?jiǎn)栴}類型。F(K,J)：節(jié)點(diǎn)J在K方向上集中載荷大小。ALFL=TRUE：非對(duì)稱矩陣集合；ALFL=FALSE：對(duì)稱矩陣集合：Ｓ：?jiǎn)卧獎(jiǎng)偠染仃?；Ｐ：載荷和內(nèi)力向量；AD：對(duì)角元素；AU：上三角元素；AL：下三角元素；JP：上〔下〕三角元素每行〔列〕最后元素；LD：一個(gè)單元中各自由度方程數(shù)。JD：列高；ID：邊界條件?！?〕數(shù)據(jù)輸入格式整數(shù)：I5格式；小數(shù)：F10.4格式?！?〕程序本體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單元四節(jié)點(diǎn)形式，計(jì)算其數(shù)值解。網(wǎng)格劃分及受力情況如下列圖：程序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調(diào)試程序得到結(jié)果為：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-