计算力学课程设计

前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法，由于其通用性和有效性，受到工程技术界的高度重视。

伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计和计算机辅助制造的重要组成部分。

由于有限元法是通过计算机实现的，因此有限元程序的编制及相关软件的研发就变得尤为重要。

从二十世纪五十年代以来，有限元软件的发展按目的和用途可分为专用软件和大型通用商业软件，而且软件往往集成了网格自动划分，结果分析和显示等前后处理功能，而且随着时间的发展，大型通用商业软件的功能由线性扩展到非线性，由结构扩展到非结构等等，这一系列强大功能的实现与运用都要求我们对有限元法的基础理论知识有较为清楚的认识以及对程序编写的基本能力有较好掌握。

一、有限元分析程序的理论基础作为弹性力学微分方程的等效积分形式，虚位移原理和虚应原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。

将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理和最小余能原理，它们本质上和等效积分的伽辽金“弱”形式相一致，这是建立弹性力学有限原方程一般表达格式的理论基础。

对于通过弹性力学变分原理建立的弹性力学问题有限元方法，其未知场变量是位移，以节点位移为基本未知量，并以最小位能原理为基础建立的有限单元为位移元。

弹性力学平面问题有限元分析表达格式的建立步骤一般为：1、对所给的模型进行单元离散，一般采用三角形单元或四边形单元，不过由于四边形单元具有更高精度，应用更为普遍，一般而言一次或二次单元已足以满足精度要求。

离散后对单元进行编号，并且给定单元各节点的整体编码以及局部坐标系下按逆时针局部编码。

2、在局部坐标系下对各单元构造形函得到由单元各节点位移表示的单元位移形式，进而得到单元刚度矩阵，利用等参元性质和雅阁比矩阵进行组集，建立整体刚度矩阵。

3、建立单元等效结点载荷列阵并组集成结构节点等效载荷列阵4、得到单元各节点位移与结构节点等效载荷列阵的线代方程组，求解得到位移，进而可得应力应变。

根据以上一般步骤，编制相应的计算机程序并采用数值积分方法处理有关数学计算就可以得到完整的弹性力学问题有限元分析程序。

二、平面问题有限元分析教学程序本程序可对二维弹性力学问题行有限元分析计算。

计算采用的单元形式为四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽存储方法，最终输出结果为单元各节点的位移。

C---------------------------------------------------------------------C |C FEA2DP----A FINITE ELEMENT ANAYLYSIS PROGRAM FOR |C 2D ELASTIC PROBLEMS |C TANGENT MATRIX IS STORED WITH VARIOUS BAND METHOH |C THIS PROGRAM IS USED TO DEMONSTRATE THE USAGE OFVRIOUS BAND |C STRORAGE SCHEM OF SYMMETRIC AND UNSYMMETRIC TANGENTSMATRIX |C |C XIAOBAO HE |C AT CHONGQING UNIVERSITY,CHINA.(12/6/2012) |C |C---------------------------------------------------------------------CPROGRAM FEA2DPCC A(1)-A(N1-1):X(NDM,NUMAP); A(N1)-A(N2-1):F(NDF,NUMAP)C A(N2)-A(N3-1):B(NEQ); A(N3)-A(N4-1):AD(NEQ)C A(N4)-A(N5-1):AL(NED); A(N5)-A(N6-1):AU(NAD)CC IA(1)-IA(N1-1):IX(NEN1,NUMEL); IA(N1)-IA(N2-1):ID(NDF,NUMNP)C IA(N2)-IA(N3-1):JD(NDF*NUMNP); IA(N3)-IA(N4-1): IDL(NEN*NUMEL*NDF)CIMPLICIT REAL*8(A-H,O-Z)DIMENSION A(10000000),IA(100000)CHARACTER*80 HEADCOMMON /CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON /SDATA/NDF,NDM,NEN1,NSTCOMMON /IOFILE/IOR,IOWNMAXM=10000000IMAXM=100000IOR=1IOW=2CC Open files for data input and outputCOPEN(IOR,FILE='INPUT.DAT',FORM='FORMATTED') OPEN(IOW,FILE='OUTPUT.DAT')CC....READ TITLECREAD(IOR,'(A)') HEADWRITE(IOW,'(A)') HEADCC....READ AND PRINT CONTROL INFORMATIONCC NUMNP:number of nodesC NUMEL:number of elementsC NUMMAT:number of material typesC NLOAD:number of loadsC NDM:number of coordinates of each nodeC NDF:number of degrees of freedomC NEN:number of nodes in each elementCREAD(IOR,'(7I5)')NUMNP,NUMEL,NUMMAT,NLOAD,NDM,NDF,NENWRITE(IOW,2000) NUMNP,NUMEL,NUMMAT,NLOAD,NDM,NDF,NENCC....SET POITERS FOR ALLOCATION OF DATA ARRAYSCNEN1=NEN+4NST=NEN*NDFNNEQ=NDF*NUMNPCN1=NDM*NUMNP+1N2=N1+NDF*NUMNP+1CI1=NEN1*NUMEL+1I2=I1+NDF*NUMNP+1I3=I2+NDF*NUMNP+1I4=I3+NUMEL*NEN*NDF+1CC....CALL MESH INPUT SUBROUTINE TO READ ALL MESH DATACCALL PMESH (A(1),A(N1),IA(1),IA(I1),NDF,NDM,NEN1,NLOAD)CPUTE PROFILECCALLPROFIL(IA(I2),IA(I3),IA(I1),IA(1),NDF,NEN1,NAD)CN3=N2+NEQ+1N4=N3+NEQ+1N5=N4+NAD+1N6=N5+NAD+1CC The lengths of Real and Integer arraysCWRITE (IOW,2222) N6,I4CC The length o9f exceeds the limitationCIF (N6>NMAXM.OR.I4>IMAXM) THENIF (N6>NMAXM) WRITE (IOW,3333) N6,NMAXMIF (I4>IMAXM) WRITE (IOW,4444) I4,IMAXMSTOPEND IFCPUTE AND ASSEMBLE ELEMENT ARRAYSCCALLASSEM(NAD,IA(1),IA(I1),IA(I2),A(1),A(N2),A(N3),A(N4),A(N5)) CC FORM LOAD VECTORCCALL PLOAD(IA(I1),A(N1),A(N2),NNEQ,NEQ)CC....TRIANGULAR DECOMPOSITION OF A MATRIX SYORED IN PROFILE FORMCCALLDATRI(NDF,NUMNP,IA(I2),NEQ,NAD,.FALSE.,A(N3),A(N5),A(N5)) CCC FOR UNSYMMETRIC TANGENT MATRIXC CALL DATRI(NDF,NUMNP,IA(I2),NEQ,NAD,.TRYE.,A(N3),A(N4),A(N5)) CC SOLVE EQUATIONSCCALLDASOL(NDF,NUMNP,A(N2),IA(I2),NEQ,NAD,AENGY,A(N3),A(N5),A(N5 ))CC FOR UNSYMMETRIC TANGENT MATRIXC CALL DASOL(NDF,NUMNP,A(N2),IA(I2),NEQ,NAD,AENGY,A(N3),A(N4),A(N5 ))CC OUTPUT NODAL DISPLACEMENTSCCALL PRTDIS(IA(I1),A(N2),NDF,NUMNP,NEQ)CC....CLOSE INPUT AND OUTPUT FILES; DESTROY TEMPORARY DISK FILESCCLOSE(IOR)CLOSE(IOW)CC....INPUT/OUTPUT FORMATSC1000 FORMAT (20A4)2000 FORMAT(//X 5X,'NUMBER OF NODAL POINTS =',I6/1 5X,'NUMBER OF ELEMENTS =',I6/2 5X,'NUMBER OF MATERIAL SETS =',I6/2 5X,'NUMBER OF NODAL LOADS =',I6/3 5X,'DIMENSION OF COORDINATE SPACE =',I6/4 5X,'DEGREE OF FREEDOMS/NODE =',I6/5 5X,'NODES PER ELEMENT (MAXIMUM) =',I6/)2222 FORMAT (//,10X,'THE LENGTH OF REAL ARRAY IS',I10,/,1 10X,'THE LENGTH OF INTEGER ARRAY IS',I10)3333 FORMAT (//,10X,'THE LENGTH OF REAL ARRAY',I10,'EXCEED THE',1 'MAXIMUN VALUE',I10)4444 FORMAT(//,10X,'THE LENGTH OF INTEGER ARRAY',I10,'EXCEED THE',1 'MAXIMUN VALUE',I10)CSTOPENDCCSUBROUTINE PMESH(X,F,IX,ID,NDF,NDM,NEN1,NLOAD)CC....DATA INPUT ROUTINE FOR MESH DESCRIPTIONCIMPLICIT REAL*8(A-H,O-Z)DIMENSIONX(NDM,NUMNP),F(NDF,NUMNP),ID(NDF,NUMNP),IX(NEN1,NUMEL)1 ,J0(NLOAD)COMMON/BDATA/HEAD(20)COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/MATER/EE,XNU,ITYPECOMMON/IOFILE/IOR,IOWCC....INPUT CONSTRAIN CODES ANDNODAL COORDINATE DATACC ID(K,J):CONSTRAIN CODE OF Kth degree of freedom of node J,=0:free,=1:fixedC X(K,J):Kth coordinate of node JcDO I=1,NUMNPREAD(IOR,'(3I5,2F10.4)')J,(ID(K,J),K=1,NDM),(X(K,J),K=1,NDM )END DOCWRITE(IOW,'(//,17HNODAL COORDINATES,/)')DO I=1,NUMNPWRITE(IOW,'(3I5,2F10.4)')I,(ID(K,I),K=1,NDM),(X(K,I),K=1,ND M)END DOCC.....ELEMENT DATA INPUTCC IX(K,J):global node number of Kth node in element JCDO I=1,NUMELREAD(IOR,'(9I5)')J,(IX(K,J),K=1,NEN)END DOCCWRITE(IOW,'(//,18HELEMENT DEFINITION,/)')DO I=1,NUMELWRITE(IOW,'(9I5)')I,(IX(K,I),K=1,NEN)END DOCC.....MATERIAL DATA INPUTCC EE:Young's nodulus,XNU:Poisson ratioC ITYPE:Type of problem,=1:Plane stress,=2:Plane strain,=3:Axi-symmetricCREAD(IOR,'(2F10.4,I5)')EE,XNU,ITYPEWRITE(IOW,'(//,19HMATERIAL PROPERTIES,/)')WRITE(IOW,'(2(E10.4,5X),I5)')EE,XNU,ITYPECC.....FORCE/DISP DATA INPUTCC F(K,J):Concentrate load at node J in K directionCF=0.0D0DO I=1,NLOADREAD(IOR,'(I5,2F10.4)')J, (F(K,J),K=1,NDF)J0(I)=JEND DOCWRITE(IOW,'(//,20HAPPLIED NODAL FORCES,/)') DO I=1,NLOADWRITE(IOW,'(I5,2F10.4)')J0(I),(F(K,J0(I)),K=1,NDF)END DOCRETURNCC FORMAT STATEMENTSC2000 FORMAT(' Mesh 1>',$)3000 FORMAT(1X,'**WARNING**ELEMENT CONNECTIONS NECESSARY'1'TO USE BLOCK IN MACRO PROGRAM')4000 FORMAT('**CURRENT PROBLEM VALIES**'/I6,'NODES,', 1I5,'ELMTS,',I3,'MATLS,',I2,'DIMS,',I2,'DOF/NODE,',2I3,'NODES/ELMT')ENDCC SUBROUTINE ASSEM(NAD,IX,ID,JD,X,B,AD,AL,AU)CC Call element subuoutine and assemble global tangent matrixCIMPLICIT REAL*8(A-H,O-Z)DIMENSIONILX(NEN),XL(NDF,NEN),LD(NDF,NEN),S(NST,NST),P(NST)DIMENSION IX(NEN1,NUMEL),ID(NDF,NUMNP),JD(NDF*NUMNP)DIMENSION X(NDM,NUMNP),B(NEQ),AD(NEQ),AL(NAD),AU(NAD)COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/SDATA/NDF,NDM,NEN1,NSTCNEL=NENCC Element loopCDO 320 N=1,NUMELS=0.0D0 !element stiffness matrixP=0.0D0 !nodal forceNE=NDO 310 I=1,NENILX(I)=IX(I,NE) !current element difinitionDO K=1,NDMXL(K,I)=X(K,ILX(I)) !nodal coords in current elementEND DOKK=ILX(I)DO K=1,NDFLD(K,I)=ID(K,KK) !equation numbersEND DO310 CONTINUECC Call element libCCALL ELMT01(XL,ILX,S,P,NDF,NDM,NST)CC Assemble tangent matrix and load vector if neededCCALLDASBLY(NDF,NAD,S,P,LD,JD,NST,B,AD,AL,AU)C320 CONTINUE !end element loopCRETURNENDCC SUBROUTINE DASBLY(NDF,NAD,S,P,LD,JP,NS,B,AD,AL,AU)CC.....ASSUMBLE THE SYMMETRIC OR UNSYMMETRIC ARRAYS FOR 'DASOL'CIMPLICIT REAL*8(A-H,O-Z)C LOGICAL ALFL,AUFL,BFLDIMENSION AD(NEQ),AL(NAD),AU(NAD)DIMENSIONLD(NS),JP(NDF*NUMNP),B(NEQ),S(NS,NS),P(NS)COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/IOFILE/IOR,IOWCC ALFL=TRUE:FOR UNSUMMETRIC MATRIX ASSEMBLEC ALFL=FALSE:FOR SUMMETRIC MATRIX ASSUMBLEC S:ELEMENT STIFFNESS MATRIXC P:LOAD OR INTERNAL FORCE VECTORC AD:DIAGONAL ELEMENTSC AU:UPPER TRIANGLE ELEMENTSC AL:LOWER TRIANGLE ELEMENTSC JP:POINTER TO LAST ELEMENT IN EACH ROW/COLUMN OF AL/AU RESPECTIVELC LD:EQUATION NUMBERS OF EACH FREEDOM DEGREE IN AN ELEMENT(GET FROMC ID)CC.....LOOP THROUGH THE ROWS TO PERFORM THE ASSEMBLYCDO 200 I=1,NSII=LD(I)IF(II.GT.0)THENC IF(AUFL)THEN !Assemble stiffness matrixCC.....LOOP THROUGH THE COLUMNS TO PERFORM THE ASSEMBLYCDO 100 J=1,NSIF(LD(J).EQ.II)THENAD(II)=AD(II)+S(I,J)ELSEIF (LD(J).GT.II)THENJC=LD(J)JJ=II+JP(JC)-JC+1AU(JJ)=AU(JJ)+S(I,J)C IF(ALFL) AL(JJ)=AL(JJ)+S(J,I) !UnsymmetricENDIF100 CONTINUEENDIF200 CONTINUECRETURNENDCC SUBROUTINE DASOL(NDF,NUMNP,B,JP,NEQ,NAD,ENERGY,AD,AL,AU)CC.....SOLUTION OF SYMMETRIC EQUATIONS IN PROFILE FORMC.....COEFFICIENT MATRIX MUST BE DECOMPOSED INTO ITS TRIANGULARC.....FACTORS USING DATRI BEFORE USING DASOL.CC JP:POINTER TO LAST ELEMENT IN EACH ROW/COLUMN OF AL/AU RESPECTIVELYCIMPLICIT REAL*8(A-H,O-Z)DIMENSION AD(NEQ),AL(NAD),AU(NAD)DIMENSION B(NEQ),JP(NDF*NUMNP)COMMON/IOFILE/IOR,IOWDATA ZERO/0.0D0/CC.....FIND THE FIRST NONZERO ENTRY IN THE RIGHT HAND SIDECDO IS=1,NEQIF(B(IS).NE.ZERO)GO TO 200END DOWRITE(IOW,2000)RETURNC200 IF(IS.LT.NEQ)THENCC.....REDUCE THE RIGHT HAND SIDECDO 300 J=IS+1,NEQJH=JP(J)-JRIF(JH.GT.0)THENB(J)=B(J)-DOT(AL(JR+1),B(J-JH),JH)ENDIF300 CONTINUEENDIFCC.....MULTIPLY BY INVERSE OF DIAGONAL ELEMENTS CENERGY=ZERODO 400 J=IS,NEQBD=B(J)B(J)=B(J)*AD(J)ENERGY=ENERGY+BD*B(J)400 CONTINUECC.....BACKSUBSTITUTIONCIF(NEQ.GT.1)THENDO 500 J=NEQ,2,-1JR=JP(J-1)IF(JH.GT.0)THENCALL SAXPB(AU(JR+1),B(J-JH),-B(J),JH,B(J-JH))ENDIF500 CONTINUEENDIFCRETURNC2000 FORMAT('**DASOL WARNING 1 ** ZERO RIGHT-HAND-SIDE VECTOR')ENDCC SUBROUTINE DATEST(AU,JH,DAVAL)CC.....TEST FOR RANKCIMPLICIT REAL*8(A-H,O-Z)DIMENSION AU(JH)CDAVAL=0.0D0CDO J=1,JHDAVAL=DAVAL+ABS(AU(J))END DOCRETURNENDC SUBROUTINE DATRI(NDF,NUMNP,JP,NEQ,NAD,FLG,AD,AL,AU)CC.....TRIANGULAR DECOMPOSITION OF A MATRIX STORED IN PROFILE FORMCIMPLICIT REAL*8(A-H,O-Z)LOGICAL FLGDIMENSION JP(NDF*NUMNP),AD(NEQ),AL(NAD),AU(NAD)COMMON/IOFILE/IOR,IOWCC.....N.B. TOL SHOULD BE SET TO APPROXIMATE HALF-WORD PRECISIONCDATA ZERO,ONE/0.0D0,1.0D0/,TOL/0.5D-07/CC.....SET INITIAL VALUES FOR CONTDITIONING CHECKCDIMX=ZERODIMN=ZEROCDO J=1,NEQDIMN=MAX(DIMN,ABS(AD(J)))END DODFIG=ZEROCC.....LOOP THROUGH THE COLUMS TO PERFOM THE TRIANGULAR DECOMPOSITIONCJD=1DO 200 J=1,NEQJR=JD+1JD=JP(J )JH=JD-JRIF(JH.GT.0)THENIS=J-JHIE=J-1CC.....IF DIAGONAL IS ZERO COMPUTE A NORM FOR SINGULARITY TESTCIF(AD(J).EQ.ZERO)CALL DATEST(AU(JR),JH,DAVAL)DO 100 I=IS,IEJR=JR+1ID=JP(I)IH=MIN(ID-JP(I-1),I-IS+1)IF(IH.GT.0)THENJRH=JR-IHIDH=ID-IH+1AU(JR)=AU(JR)-DOT(AU(JRH),AL(IDH),IH)IF(FLG)AL(JR)=AL(JR)-DOT(AL(JRH),AU(IDH),IH)ENDIF100 CONTINUEENDIFCC...REDUCE THE DIAGONALCIF(JH.GE.0)THENDD=AD(J)JR=JD-JHJRH=J-JH-1CALL DREDU(AL(JR),AU(JR),AD(JRH),JH+1,FLG,AD(J)) CC.....CHECK FOR POSSIBLE ERRORS AND PRINT WARNINGSCIF(ABS(AD(J)).LT.TOL*ABS(DD))WRITE(IOW,2000)J IF(DD.LT.ZERO.AND.AD(J).GT.ZERO)WRITE(IOW,2001)JIF(DD.LT.ZERO.AND.AD(J).LT.ZERO)WRITE(IOW,2001)JIF(AD(J).EQ.ZERO)WRITE(IOW,2002)JIF(DD.EQ.ZERO.AND.JH.GT.0)THENIF(ABS(AD(J)).LT.TOL*DAVAL)WRITE(IOW,2003)JENDIFENDIFCC.....STORE RECIPROCAL OF DIAGONAL,COMPUTE CONDITION CHECKSCIF(AD(J).NE.ZERO)THENDIMX=MAX(DIMX,ABS(AD(J)))DIMN=MIN(DIMN,ABS(AD(J)))DIFG=MAX(DIFG,ABS(DD/AD(J)))AD(J)=ONE/AD(J)ENDIF200 CONTINUECC.....PRINT CONDITIONING INFORMATIONCDD=ZEROIF(DIMN.NE.ZERO)DD=DIMX/DIMNIFIG=DLOG10(DIFG)+0.6WRITE(IOW,2004)DIMX,DIMN,DD,IFIGCC.....FORMATSC2000 FORMAT('**DATRI WARNING 1**LOSS OF AT LEAST 7 DIGITS IN',1 'REDUCING DIAGONAL OF EQUATION',I5)2001 FORMAT('**DATRI WARNING 2**SIGN OF DIAGONAL CHANGED WHEN',1 'REDUCING EQUATION',I5)2002 FORMAT('**DATRI WARNING 3**REDUCED DIAGONAL IS ZERO FOR',1 'EQUATION',I5)2003 FORMAT('**DATRI WARNING 4**RANK FAILURE FOR ZEROUNREDUCED',1 'DIAGONAL IN EQUATION',I5)2004 FORMAT(//'CONDITION CHECK:D-MAX',E11.4,';D-MIN',E11.4,1 ';RATIO',E11.4/' MAXIMIM NO.DIAGONAL DIGITS LOST:',I3)2005 FORMAT('COND CK:DMAX',1P1E9.2,';DMIN',1P1E9.2,1 ';RATIO',1P1E9.2)ENDCC SUBROUTINE DREDU(AL,AU,AD,JH,FLG,DJ)CC.....REDUCE DIAGONAL ELEMENT IN TRIANGULAR DECOMPOSITIONCIMPLICIT REAL*8(A-H,O-Z)LOGICAL FLGDIMENSION AL(JH),AU(JH),AD(JH)CDO J=1,JHUD=AU(J)*AD(J)DJ=DJ-AL(J)*UDAU(J)=UDEND DOCC.....FINISH COMPUTATION OF COLUMN OF AL FOR UNSYMMETRIC MATRICESCIF(FLG)THENDO J=1,JHAL(J)=AL(J)*AD(J)END DOENDIFCRETURNENDCC SUBROUTINE PROFIL(JD,IDL,ID,IX,NDF,NEN1,NAD)CPUTE PROFILE OF GLOBAL ARRAYSCIMPLICIT REAL*8(A-H,O-Z)DIMENSIONJD(NDF*NUMNP),IDL(NUMEL*NEN*NDF),ID(NDF,NUMNP),1 IX(NEN1,NUMEL)COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/FRDATA/MAXFCOMMON/IOFILE/IOR,IOWCC JD: COLUMH HIGHT (ADDRESS OF DIAGONAL ELEMENTS)C ID: BOUNDARY CONDITION CODES BEFORE THIS SUBROUTINE'S RUNNINGC ID: EQUATION NUMBERS IN GLOBAL ARRAY (EXCLUDED RESTRAINED NODES)C AFTER RUNNINGC IDL:ELEMENT SEARCH ORDERC NAD:TOTAL NUMBER OF NON-ZERO ELEMENTS EXCEPT DIAGONAL ELEMENTSC IN GLOBAL TANGENT MATRIXCC.....SET UP THE EQUATION NUMBERSCNEQ=0CDO 10 K=1,NUMNPDO 10 N=1,NDFJ=ID(N,K)IF(J.EQ.0)THENNEQ=NEQ+1ID(N,K)=NEQELSEID(N,K)=0ENDIF10 CONTINUECPUTE COLUMN HEIGHTSCCALL PCONSI(JD,NEQ,0)CDO 50 N=1,NUMELMM=0NAD=0DO 30 I=1,NENII=IABS(IX(I,N))IF(II.GT.0)THENDO 20 J=1,NDFJJ=ID(J,II)IF(JJ.GT.0)THENIF(MM.EQ.0) MM=JJMM=MIN(MM,JJ)NAD=NAD+1IDL(NAD)=JJENDIF20 CONTINUEENDIF30 CONTINUEIF(NAD.GT.0)THENDO 40 I=1,NADII=IDL(I)JJ=JD(II)JD(II)=MAX(JJ,II-MM)40 CONTINUEENDIF50 CONTINUECPUTE DIAGONAL POINTERS FOR PROFILE CNAD=0JD(1)=0IF(NEQ.GT.1)THENDO 60 N=2,NEQJD(N)=JD(N)+JD(N-1)60 CONTINUENAD=JD(NEQ)ENDIFCC.....SET ELEMENT SEARCH ORDER TO SEQUENTIALCDO 70 N=1,NUMELIDL(N)=N70 CONTINUECC.....EQUATION SUMMARYCMAXF=0MM=0IF(NEQ.GT.0)MM=(NAD+NEQ)/NEQWRITE(IOW,2001)NEQ,NUMNP,MM,NUMEL,NAD,NUMMAT CRETURNC2001 FORMAT(5X,'NEQ =',I5,5X,'NUMNP=',I5,5X,'MM =',I5,/5X,1 'NUMEL=',I5,5X,'NAD =',I5,5X,'NUMMAT=',I5/)ENDCC SUBROUTINE SAXPB(A,B,X,N,C)CC.....VECTOR TIMES SCALAR ADDED TO SECOND VECTORCIMPLICIT REAL*8(A-H,O-Z)DIMENSION A(N),B(N),C(N)CDO K=1,NC(K)=A(K)*X+B(K)END DOCRETURNENDCC SUBROUTINE PCONSI(IV,NN,IC)CC.....ZERO INTEGER ARRAYCDIMENSION IV(NN)CDO N=1,NNIV(N)=ICEND DOCRETURNENDCC SUBROUTINE ELMT01(XL,ILX,S,P,NDF,NDM,NST)CC.....PLANE LINEAR ELASTIC ELEMENT ROUTINEC ITYP=1:PLANE STRESSC =2:PLANE STRAINC =3:AXISYMMETRICCIMPLICIT REAL*8(A-H,O-Z)DIMENSION XL(NDM,NEN),ILX(NEN),SIGR(6)DIMENSIOND(18),S(NST,NST),P(NST),SHP(3,9),SG(16),TG(16),WG(16) CHARACTER WD(3)*12COMMON/CDATA/NUMNP,NUMEL,NUMMAT,NEN,NEQCOMMON/MATER/EE,XNU,ITYPECOMMON/FCIOFILE/IOR,IOWDATA WD/'PLANE STRESS','PLANE STRAIN','AXISYMMETRIC'/INTEGER::ALP=1CC XL(NDM,NEN):COORDS OF EACH NODE IN CURRENT ELEMENTC ILX(NEN):ELEMENT DEFINITION OF CURRENT ELEMENTC D(18):MATERIAL PROPERTIESC S(NST,NST):ELEMENT STIFFNESS MATRIXC P(NS):NODAL FORCE AND/OR INTERNAL FORCEC SHP(3,9):SHAPE FUCTION AND ITS DERIVATIVESC SG(16),TG(16),WG(16):WEIGHT COEFFICIENTS OF GUASS INTEGRATIONC L,K:INTERGRATION POINTSCL=2K=2E=EENEL=NENCC D(14):THICKNESS;D(11),D(12):BODY FORCEC.....SET MATERIAL PARAMETER TYPE AND FLAGSCITYP=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)=ITYPD(16)=ED(17)=XNUD(18)=-XNU/EL=MIN(4,MAX(1,L))K=MIN(4,MAX(1,K))D(5)=LD(6)=KC D(9)=T0C D(10)=E*ALP/(1.-J*XNU)LINT=0CWRITE(IOW,2000)WD(ITYP),D(16),D(17),D(4),L,K,D(14),1 D(11),D(12),ALPCC.....STIFFNESS/RESIDUAL COMPUTATIONCL=KCC cmopute cordinates and weights of integration pointsC SG,TG:Coords;WG=Wp*WqCIF(L*L.NE.LINT)CALL PGUASS(L,LINT,SG,TG,WG)CPUTE INTEGRALS OF SHAPE FUNCTIONSCDO 340L=1,LINTCC compute shape function and their derivatives to localc and global coordinate systemCCALL SHAPE (SG(L),TG(L),XL,SHP,XSJ,NDM,NEN,ILX,.FALSE.)CC compute global coordinates of integration points CXX=0.0YY=0.0DO J=1,NENXX=XX+SHP(3,J)*XL(1,J)YY=YY+SHP(3,J)*XL(2,J)END DOXSJ=XSJ*WG(L)*D(14) !XSJ=|J|(Sp,Tq)*Wp*Wq*tCPUTE JACOBIAN CORRECTIONC for plane stress and strain problemsCIF(ITYP.LE.2)THENDV=XSJXSJ=0.0ZZ=0.0C SIGR4=-D(110*DV !D(11)BODY FORCEELSECC for axisymmetric problemCDV=XSJ*XX*3.1415926*2ZZ=1./XXC SIGR4=SIGR(4)*XSJ-D(11)*DVEND IFJ1=1CC.....LOOP OVER ROWSCDO 330 J=1,NELW11=SHP(1,J)*DVW12=SHP(2,J)*DVC W22=SHP(3,J)*XSJW22=SHP(3,J)*DV*ZZCPUTE THE INTERNAL FORCESC out of balance forceCC P(J1)-(SHP(1,J)*SIGR(1)+SHP(2,J)*SIGR(2))*DVC 1 -SHP(3,J)*SIGR4CP(J1+1)=P(J1+1)-(SHP(1,J)*SIGR(2)+SHP(2,J)*SIGR(3))*DVC 1 +D(12)*SHP(3,J0*DV !D(12) BODAY FORCECC....LOOP OVER COLUMNS(SYMMETRY NOTED)C compute stiffness matrixCK1=J1A11=D(1)*W11+D(2)*W22A21=D(2)*W11+D(1)*W22A31=D(2)*(W11+W22)A41=D(3)*W12A12=D(2)*W12A32=D(1)*W12A42=D(3)*W11DO 320 K=J,NELW11=SHP(1,K)W12=SHP(2,K)W22=SHP(3,K)*ZZS(J1,K1)=S(J1,K1)+W11*A11+W22*A21+W12*A41S(J1+1,K1)=S(J1+1,K1)+(W11+W22)*A12+W12*A42 S(J1,K1+1)=S(J1,K1+1)+W12*A31+W11*A41S(J1+1,K1+1)=S(J1+1,K1+1)+W12*A32+W11*A42 K1=K1+NDF320 CONTINUEJ1=NDF+J1330 CONTINUE340 CONTINUEC.....MAKE STIFFNESS SYMMETRICCDO 360 J=1,NSTDO 360 K=J,NSTS(K,J)=S(J,K)360 CONTINUECRETURNCC.....FORMATS FOR INPUT-OUTPUTC1000 FORMAT(3F10.0,3I10)1001 FORMAT(8F10.0)2000 FORMAT(/5X,A12,'LINEAR ELASTIC ELEMENT'//110X,'MODULUS',E18.5/10X,'POISSIONRATIO',F8.5/10X,'DENSITY',E18.5/210X,'GUASS PTR/DIR',I3/10X,'STRESS PTS',I6/10X,'THICKNESS',E16.5/310X,'1-GRAVITY',E16.5/10X,'2-GRAVITY',E16.5/10X,'ALPHA',E20.5/410X,'BASE TEMP',E16.5/)2001 FORMAT(5X,'ELEMENT STRESSES'//'ELMT 1-COORD',2X,'11-STRESS',2X,1'12-STRESS',2X'22-STRESS',2X,'33-STRESS',3X,'1-STRESS',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/) 5000 FORMAT('INPUT:E,NU,RHO,PTS/STIFF,PTS/STRE',1',TYPE(1=STRESS,2=STRAIN,3=AXISM)',/3X,'>',$)5001 FORMAT('INPUT:THICKNESS,1-BODY FORCE,1-BODY FORCE,ALPHA,'1 ,'TEMP-BASE'/3X,'>',$)ENDCC SUBROUTINE SHAPE(SS,TT,XL,SHP,XSJ,NDM,NEL,ILX,FLG)CC.....SHAPE FUNCTION ROUTINE FOR TWO DIMENSIONAL ELEMENTS CIMPLICIT REAL*8(A-H,O-Z)LOGICAL FLGDIMENSION XL(NDM,NEL),S(4),T(4),ILX(NEL)DIMENSION SHP(3,NEL),XS(2,2),SX(2,2)DATA S/-0.5D0,0.5D0,0.5D0,-0.5D0/,1 T/-0.5D0,-0.5D0,0.5D0,0.5D0/CC.....FORM 4-NODE QUATRILATERAL SHAPE FUNCTIONSC NEL:number of nodes per elementCDO 100 I=1,4SHP(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)100 CONTINUECC.....FORM TRIANGLE BY ADDING THIRD AND FOURTH TOGETHER C for triangle elementCIF (NEL.EQ.3)THENDO I=1,3SHP(I,3)=SHP(I,3)+SHP(I,4)END DOEND IFCC.....ADD QUADRATIC TERMS IF NECESSARYC for element with more than 4 nodescIF (NEL.GT.4)CALL SHAP2(SS,TT,SHP,ILX,NEL) CC.....CONSTRUCT JACOBIAN AND ITS INVERSECDO 125 I=1,2DO 125 J=1,2XS(I,J)=0.0DO 120 K=1,NELXS(I,J)=XS(I,J)+XL(I,K)*SHP(J,K)120 CONTINUE125 CONTINUECC XSJ:determinate of Jacob matrixcXSJ=XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1)CIF(FLG)RETURNCC FLG=false:form global derivativesCIF(XSJ.LE.0.0D0)XSJ=1.0SX(1,1)=XS(2,2)/XSJSX(2,2)=XS(1,1)/XSJSX(1,2)=-XS(1,2)/XSJSX(2,1)=-XS(2,1)/XSJCC.....FORM GLOBAL DERIVATIVESCDO 130 I=1,NELTP=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)=TP130 CONTINUECRETURNENDC SUBROUTINE SHAP2(S,T,SHP,ILX,NEL) CC.....ADD QUADRATIC FUNCTION AS NECESSARY CIMPLICIT REAL*8(A-H,O-Z)DIMENSION SHP(3,9),ILX(NEL)CS2=(1.-S*S)/2.T2=(1.-T*T)/2.DO 100 I=5,9DO 100 J=1,3SHP(J,I)=0.0100 CONTINUECC.....MIDSIZE NODES (SERENIPITY)CIF(ILX(5).EQ.0)GO TO 101SHP(1,5)=-S*(1.-T)SHP(2,5)=-S2SHP(3,5)=S2*(1.-T)101 IF(NEL.LT.6)GO TO 107IF(ILX(6).EQ.0)GO TO 102SHP(1,6)=T2SHP(2,6)=-T*(1.+S)SHP(3,6)=T2*(1.+S)102 IF(NEL.LT.7)GO TO 107IF(ILX(7).EQ.0)GO TO 103 SHP(1,7)=-S*(1.+T)SHP(2,7)=S2SHP(3,7)=S2*(1.+T)103 IF(NEL.LT.8)GO TO 107IF(ILX(8).EQ.0)GO TO 104 SHP(1,8)=-T2SHP(2,8)=-T*(1.-S)SHP(3,8)=T2*(1.-S)CC.....INTERIOR NODE (LAGRAGIAN) C104 IF(NEL.LT.9)GO TO 107IF(ILX(9).EQ.0)GO TO 107 SHP(1,9)=-4.*S*T2SHP(2,9)=-4.*T*S2SHP(3,9)=4.*S2*T2CC.....CORRECT EDGE NODES FOR INTERIOR NODE(LAGRANGIAN) CDO 106 J=1,3DO 105 I=1,4105 SHP(J,I)=SHP(J,I)-0.25*SHP(J,9)DO 106 I=5,8106 IF(ILX(I).NE.0)SHP(J,I)=SHP(J,I)-.5*SHP(J,9)CC.....CORRECT CORNER NODES FOR PRESENSE OF MIDSIZE NODES C107 DO 108 I=1,4K=MOD(I+2,4)+5L=I+4DO 108 J=1,3108 SHP(J,I)=SHP(J,I)-0.5*(SHP(J,K)+SHP(J,L))RETURNENDCCSUBROUTINE PGUASS(L,LINT,R,Z,W)CC....GUASS POINTS AND WEIGHTS FOR TWO DIMENSIONSCIMPLICIT REAL*8(A-H,O-Z)DIMENSION LR(9),LZ(9),LW(9),R(16),Z(16),W(16)C COMMON/ELDATA/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/ DATA LW/4*25,4*40,64/cc Lint:Number of integration pointsc R,Z:Coordinates of Integration pointsc W:Wp*Wq,product of the two weightscLINT=L*LGO TO(1,2,3),LCC....1X1 INTEGERATIONC1 R(1)=0.Z(1)=0.W(1)=4.CRETURN。