有限元程序设计课件

有限元程序设计有限元程序设计1教学程序：FEP1李明瑞编FEP2李明瑞编ZFEP教学程序：2第一章绪论有限元程序的基本内容有限元求解器第一章绪论有限元程序的基本内容3§1.1有限元程序的基本内容有限元法的解题步骤结构离散化为有限元网格计算单元刚度矩阵引入约束条件求应力等解方程组建立总刚度矩阵及外载§1.1有限元程序的基本内容有限元法的解题步骤结构离散化4有限元程序的基本内容：数据输入阶段有限元矩阵的计算、组集和求解数据输出阶段前处理处理器后处理有限元程序的基本内容：数据输入阶段有限元矩阵的计算、组集和求5一、数据输入阶段—前处理

读入和生成数据，形成有限元网格，为有限元矩阵计算作准备。（数据或图形输入）标题和控制信息计算存贮分配节点坐标和约束信息单元信息材料信息载荷信息一、数据输入阶段—前处理读入和生成数据，形成有限6二、有限元矩阵的计算、组集和求解—处理器计算插值函数矩阵[N

]、几何矩阵[B

]、Jacob矩阵[J

]在高斯点进行数值积分，求得单刚、单元载荷组集成总刚、总载求解二、有限元矩阵的计算、组集和求解—处理器计算插值函数矩阵[7三、数据输出阶段—后处理输出结果（场变量和场变量导数等），打印结果。场变量包括：位移、温度、流场的速度势等；场变量导数包括：应力、应变、热流、流场速度势等输出形式：数据输出和图形输出三、数据输出阶段—后处理输出结果（场变量和场变量导数等），打8§1-2有限元求解器带宽法：只存储总刚矩阵的半带宽内的元素。等带宽法：每行的带宽相等，常采用二维数组存储变带宽法：每行的带宽不等，常采用一维数组存储一、带宽法变带宽分块求解思想：一个变量的消元只影响刚度矩阵的一部分元素，较大的节点分量方程组可以分成较小的局部系数矩阵按节点逐步求解。§1-2有限元求解器带宽法：只存储总刚矩阵的半带宽内的元9

波前法是另一种分块储存的变带宽法，其消元次序是按单元编号进行，边组集边消元。调入内存的单元所保留的节点称为波前节点，所消去的节点称为消元节点，消元节点是与以后调入的单元没有联系的节点，即该节点的方程已组集完。波前法的回代按消元节点的反序进行。对内存的需求取决与最大波前宽。二、波前法波前法是另一种分块储存的变带宽法，其消10V-FortranFortran语言的基本格式变量基本语句子例子程序函数子程序其它功能、模块V-FortranFortran语言的基本格式11有限元程序设计课件12有限元程序设计课件13有限元程序设计课件14第二章FEP2程序的总体设计与输入数据FEP2的设计任务FEP2的结构设计FEP2的主程序FEP2的主控程序FEP2的数据输入格式与程序实现第二章FEP2程序的总体设计与输入数据FEP215§2-1FEP2的设计任务1.FEP2的简要题目说明FEP2是一个具有通用性的教学程序，可用于计算一般的线性静力学问题。已设计了平面梁单元与平面3-9节点元两种单元，但留下接口。2.支持软件与硬件FORTRAN77以上编译器、各种微机3.FEP2的功能§2-1FEP2的设计任务1.FEP2的简要题目说163.FEP2的功能1）输入文件名由用户自由定义，但限制为4个字符，输入文件扩展名为“DAT”，输出文件扩展名为“OUT”。2）节点坐标、单元信息等具有线性自动生成功能。3）可以处理多种工况、多种类型单元组合结构问题。4）可以处理固定约束和指定位移。5）采用变带宽存储、三角分解法求解刚度方程。6）为多种单元留下接口。3.FEP2的功能17§2-2FEP2的结构设计FEP2：由形式主程序、主控程序、功能模块组成。主程序的主要功能：定义输入、输出文件名，调用主控程序PCONTR主控程序PCONTR的主要功能模块：1）内存分配2）网格生成3）变带宽存储设置4）单刚的形成与组集5）刚度矩阵的分解6）形成节点载荷7）形成与组集单元载荷8）求解位移9）求单元应力10）求结构反力§2-2FEP2的结构设计FEP2：由形式主程序、主控18FEP2程序框图FEP2（主程序）PCONTR(主控程序)—SETMEM检查内存—PROFIL形成变带宽刚度矩阵地址

—PFORM(3)形成单刚并组集—LDLT总刚的三角分解

—GENVEC形成节点载荷—PFOM(3)构造并组集单元载荷

—FORBACK前消回代求出位移—PFORM(4)计算单元应力

—PFORM(6)计算结构反力FEP2程序框图FEP2（主程序）PCONTR—S19§2-3FEP2的主程序一、文件名FEP2的文件名可由用户自由定义（限制为4个字符），通过人机交互方式确定。设计中引入以下技巧：1)规定输入文件名FI与输出文件名FO为8个字符2)规定两个字符数组NAMINP(8)，NAMOUT(8)3)用IQUIVALENCE语句使FI与NAMINP、FO与NAMOUT等价4)用DATA语句给FI和FO赋初值:“.DAT”、“.OUT”5)输入NAMINP的前四个字符作为输入输出文件名§2-3FEP2的主程序一、文件名1)规定输入文件20COMMON/PSIZE/MAXM,MAXACHARACTER*8FI,FOCHARACTERNAMINP(8),NAMOUT(8),HEAD*50COMMON/HEAD/HEAD1EQUIVALENCE(FI,NAMINP(1)),(FO,NAMOUT(1))DATAFI,FO/'.DAT','.OUT'/WRITE(*,'(A\)')'INPUTFILENAME(4LETTERSONLY):'READ(*,'(4A1)')(NAMINP(I),I=1,4)DO10I=1,410NAMOUT(I)=NAMINP(I)

COMMON/PSIZE/MAXM,MAXA21

OPEN(6,FILE=FO)OPEN(5,FILE=FI)MAXM=16000MAXA=16380CALLPCONTR(FO)END

OPEN(6,FILE=FO)22二、定义数组FEP2中设置了两个数组M(1600)和A(16380)，数组M是作为存放坐标、单元信息、载荷、位移等，数组A则是存放刚度矩阵用的，其上界与FORTRAN（编译器）版本和计算机内存有关，确定了程序的计算规模。这两个数组也可合并成一个。二、定义数组23§2-4FEP2的主控程序

PCONTR：实质上的主程序。分以下几大段：一、程序头（前12语句）SUBROUTINEPCONTR(FO)CHARACTER*8FOLOGICALERR,ASEMCOMMON/M/M(6000)COMMON/A/A(16380)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/PSIZE/MAXM,MAXACOMMON/MDATA/NUMMATCHARACTERFD*25,FORC*6,PCOMMON/PRINT/PCOMMON/ASEM/ASEMDATAFORC/'FORC'/FD/'NODALFORCE/DISPL'/§2-4FEP2的主控程序PCONTR：24FO：输出文件名ERR：逻辑变量，作出错信息控制（.TRUE.为错）ASEM：逻辑变量，作为判断是否要组集总刚数组A：存放总刚的元素数组M：存放单刚、单元信息、节点坐标等NDF—最大自由度数NDM—空间维数NEN—单元的最多节点数NJ—节点总数NE—总单元数NEQ—总方程数NUMMAT—材料类型数FD,FORC：文字变量，作为输出的标题用P：文字变量，判断是否要输出单刚FO：输出文件名25二、输入控制信息、确定主要数组的地址（14—34语句）NEN1=NEN+1每个单元信息的存储量（节点号+材料号）NNEQ=NJ*NDF总节点节点自由度数NST=NDF*NEN单刚阶数NXC=1节点坐标数组XC的首址NIX=NXC+NJ*NDM单元信息数组IX的首址NID=NIX+NEN1*NE结构方程号编码数组ID的首址ND=NID+NNEQ材料常数数组的首址NXL=ND+20*NUMMAT单元节点坐标数组XL的首址NLD=NXL+NEN*NDM单元节点自由度数组LD的首址NUL=NLD+NST单元节点位移数组UL的首址NS=NUL+NST*2单刚数组S的首址NP=NS+NST*NST总刚对角线地址向量JP的首址

……二、输入控制信息、确定主要数组的地址（14—34语句）NE26READ(5,*)NJ,NE,NUMMAT,NDF,NDM,NENWRITE(*,2000)NJ,NE,NUMMAT,NDF,NDM,NENWRITE(6,2000)NJ,NE,NUMMAT,NDF,NDM,NENNEN1=NEN+1NNEQ=NJ*NDFNST=NDF*NENNXC=1NIX=NXC+NJ*NDM

……NS=NUL+NST*2NP=NS+NST*NSTNF=NP+NSTNU=NF+NNEQNJP=NU+NNEQNEND=NJP+NNEQREAD(5,*)NJ,NE,NUMM27§2-5FEP2的数据输入格式与程序实现FEP2的数据信息（自由格式）：1）控制信息：NJ，NE，NUMMAT，NDM，NEN2）坐标信息（一组）N，NG，（XL(I)，I=1，NDM）

节点号节点号增量坐标或载荷由GENVEC生成节点坐标或载荷向量，通过NDM,CD，FC进行识别生成的是节点坐标，还是载荷向量。CD：’NODALCOORDINATES’FC：’COOR’CD：’NODALFORCE/DISPL’FC：’FORC’§2-5FEP2的数据输入格式与程序实现FEP2的数据信28

SUBROUTINEGENVEC(NDM,X,CD,FC,NJ,ERR)CHARACTERCD*18,FC*6LOGICALERRREALX(NDM,*),XL(6)ERR=.FALSE.N=0NG=0102L=NLG=NGREAD(5,*,ERR=999)N,NG,(XL(I),I=1,NDM)101IF(N.LE.0.OR.N.GT.NJ)GOTO109DO103I=1,NDM103X(I,N)=XL(I)IF(LG)104,102,104104LG=ISIGN(LG,N-L)SUBROUTINEGENVEC(NDM,X,CD29LI=(IABS(N-L+LG)-1)/IABS(LG)DO105I=1,NDM105XL(I)=(X(I,N)-X(I,L))/LI106L=L+LGIF((N-L)*LG.LE.0)GOTO102IF(L.LE.0.OR.L.GT.NJ)GOTO108DO107I=1,NDMLMLG=L-LG107X(I,L)=X(I,LMLG)+XL(I)GOTO106WRITE(6,3000)L,CDWRITE(*,3000)L,CDERR=.TRUE.GOTO102CONTINUE……LI=(IABS(N-L+LG)-1)303）单元信息（一组）L，LX，LK，（IDL(K)，K=1，NEN）单元号单元号增量材料号节点号0，0，0，……结束由ELDA子程序生成单元信息，IX(NEN1,NE)二维数组。NEN=NEN+1NE—单元总数3）单元信息（一组）314）材料信息（一组或多组）对于每一类材料，要求输入两个记录：a）材料类型号，单元类型号b）一批与单元类型相关的常数（不超过20个）对于梁单元，输入两个如下12个常数：NP,E,G,A或Ix,Iy或Iz，Rh0，As，N1，N2，Ni，W,a类型（NP=0：平面梁作用平面载荷；1：平面梁作用空间载荷）Rh0：单位长度的质量密度N1：梁左端铰为1，非铰为0N2：右端铰为1，非铰为0Ni：单元载荷类型（0，1，2，3，4）（0为无载荷）单元载荷由EQLOAD处理成等效节点载荷4）材料信息（一组或多组）对于梁单元，输入两个如下12个32对于平面元，输入6个常数：E,XNU,TH，L，K，ITH：单元厚度L：沿一个方向的高斯积分点数（2，3，4）K：沿一个方向应力输出的高斯点数I：=0平面应变0平面应力PMESH调用ELEMLIB，再调用BEAMorPLANE对于平面元，输入6个常数：PMESH调用ELEM33SUBROUTINEPMESH(XC,IX,ID,D,NEN1)CHARACTERCOOR*6,XX*18COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/MDATA/NUMMATCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONXC(NDM,*),IX(NEN1,*),ID(NDF,*),D(20,*)LOGICALERRDATAXX/'NODALCOORDINATES'/COOR/'COOR'/1CALLGENVEC(NDM,XC,XX,COOR,NJ,ERR)

SUBROUTINEPMESH(XC,IX,ID341CALLGENVEC(NDM,XC,XX,COOR,NJ,ERR)IF(ERR)WRITE(*,3000)IF(ERR)STOP2CALLELDAT(IX,NEN1)3DO304N=1,NUMMATWRITE(6,2000)MA,IELWRITE(*,2000)MA,IELIF(MA.EQ.0)GOTO4CALLPZERO(D(1,MA),20)D(20,MA)=IELCALLLEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,NST,1)304CONTINUE2000FORMAT(/'MATERIALTYPE',I4,'ELEMENTTYPE',I4/)4CALLBOUN(ID)END1CALLGENVEC(NDM,XC,XX,CO355）约束信息（包括指定位移）（一组）PMESH调用BOUN输入：N，NG，（IDL(I)，I=1，NDF）

节点号节点号增量坐标或载荷0，0，0，……（平面元4个0，梁元5个0）IDL(I)：约束信息，被约束的自由度给非零值（1或-1），其它为0。1：该自由度被约束，但不继续生成-1：该自由度被约束，要求继续生成5）约束信息（包括指定位移）（一组）366）载荷信息载荷输入仍由GENVEC完成，维数：NDF

对于指定位移，生成方式与节点载荷相同，5）输入约束信息，6）输入其值。6）载荷信息对于指定位移，生成方式与节点载荷相37第三章总刚度矩阵的变带宽存贮与求解总刚度矩阵的变带宽存贮与对角线地址LDLT三角分解自由度指示矩阵ID与对角线地址矩阵第三章总刚度矩阵的变带宽存贮与求解总刚度矩阵的变带宽存38§3-1总刚度矩阵的变带宽存贮与对角线地址总刚度矩阵的性质：对称、正定、稀疏只存贮半带宽内的元素（下三角、变带宽）稀疏对称矩阵A中的元素可用两个一维数组B和JP完全确定。B的上界为A中的下三角、变带宽内元素总数，存贮A内的元素。JP的上界为A的阶数，存贮A的各行对角线元素的序号，称为对角线元素地址数组。对应关系：A(I,J)=B(JP(I)–I+J)每一行第一个非零元素的列号Mi:M1=1,Mi=I–JP(I)+JP(I-1)+1§3-1总刚度矩阵的变带宽存贮与对角线地址总刚度矩阵的性质39SUBROUTINEPROFIL(JP,ID,IX,NEN1,NK)

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

DIMENSIONJP(1),ID(NDF,1),IX(NEN1,*),IDL(18)

NEQ=0

NAD=0

DO50N=1,NJ

DO40I=1,NDF

J=ID(I,N)

IF(J)30,20,30

20NEQ=NEQ+1

ID(I,N)=NEQ

JP(NEQ)=0

GOTO40

30NAD=NAD-1

ID(I,N)=NAD

40CONTINUE

50CONTINUESUBROUTINEPROFIL(JP,ID40C.....COMPUTEROWBANDWIDTHDO500N=1,NEDO400I=1,NENII=IX(I,N)IF(II.EQ.0)GOTO400DO300K=1,NDFKK=ID(K,II)IF(KK.LE.0)GOTO300DO200J=1,NENJJ=IX(J,N)IF(JJ.EQ.0)GOTO200DO100L=1,NDFLL=ID(L,JJ)IF(LL.LE.0)GOTO100M=MAX0(KK,LL)JP(M)=MAX0(JP(M),IABS(KK-LL))100CONTINUE200CONTINUE300CONTINUECONTINUE500CONTINUEC.....COMPUTEROWBANDWIDTH41C.....COMPUTEDIAGONALPOINTERSFORPROFILNK=1JP(1)=1IF(NEQ.EQ.1)RETURNDO600N=2,NEQ600JP(N)=JP(N)+JP(N-1)+1NK=JP(NEQ)WRITE(*,2000)NK2000FORMAT('THEMAXIMUMSTORAGEOFGLOBALSTIFFNESS=',I6/)RETURNENDC.....COMPUTEDIAGONALPOINTER42

SUBROUTINEADDSTF(LD,JP,S,NST)CHARACTERYCOMMON/PRINT/YCOMMON/A/A(16380)COMMON/CDATA/NDF,NDM,NEN,NEJ,NE,NEQCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONLD(*),JP(*),S(NST,NST)IF(Y.EQ.'Y'.OR.Y.EQ.'y')WRITE(*,2000)N,SFORMAT('ELEMENTNUMBER',I4,'ELEMENTSTIFFNESS'/(6E12.5))SUBROUTINEADDSTF(LD,J43

DO200J=1,NSTK=LD(J)IF(K.LE.0)GOTO200L=JP(K)-KDO100I=1,NSTM=LD(I)IF(M.GT.K.OR.M.LE.0)GOTO100M=L+MA(M)=A(M)+S(J,I)100CONTINUE200CONTINUERETURNENDDO200J=1,NST44§3-2LDLT三角分解求解方程：Ax=b

A—对称、正定矩阵=>Ux=yU—单位上三角矩阵对A进行三角分解=>A=LLTL—下三角矩阵，通过比较系数法确定，需要开方运算LDLT三角分解=>A=LDLTL—下三角矩阵，D—对角矩阵i=1,…,n;j=2,…,iLDLT(A,JP,N)DIMENSIONA(*),JP(*)§3-2LDLT三角分解求解方程：Ax=b45Ax=b=>LDLTx=b令LTx=y=>LDy=b

前消公式：i=2,…,n回代公式：i=n-1,…,1

SUBROUTINEFORBACK(A,B,JP,N)DIMENSIONA(*),B(*),JP(*)Ax=b=>46§3-3自由度指示矩阵ID与对角线地址矩阵ID(NDF,NJ)：非零值（1或-1）为被约束的自由度，0为活化自由度。在PROFIL中，改造数组ID，原来的零元素（活化自由度）用正数计，原来的非零元素（非活化自由度）用负数计，全部零元素（活化）的个数记为NEQ（总刚矩阵的阶数）。§3-3自由度指示矩阵ID与对角线地址矩阵ID(NDF,47计算刚度矩阵的带宽的方案比较：1）以各节点自由度为研究对象，利用单元信息IX，找出与该点联系的节点，然后比较节点自由度之最大者。运算量：NJ*NDF*NE*NEN*NDF2）逐一研究各单元，分别取出各节点的各个自由度与该单元中其它节点的自由度比较，差值最大者再与所研究的自由度之带宽值（初值为零）比较，取其中最大者代替原有的带宽值。运算量：NE*NEN*NDF*NEN*NDF两者运算量比：NJ:NE选方案2带宽=自由度之最大差值+1暂存放于JP中计算刚度矩阵的带宽的方案比较：1）以各节点自由度为研究对象，48若已知对角线地址数组JP，则第I行的带宽：JP(1)=1,JP(I)–JP(I–1)

若JP存放的是带宽,则对角线地址数组：JP(1)=JP(1)=1JP(2)=JP(2)+JP(1)……JP(I+1)=JP(I+1)+JP(I)总刚矩阵(按活化自由度存贮)全部存贮量:NK=JP(NEQ)若已知对角线地址数组JP，则第I行的带宽：49第四章单刚与载荷的组集，指定约束的处理，单元分析的预处理子程序PFORM的功能子程序MODIFY完成将指定位移转化成等价内力将单刚组集成总刚子程序ADDSTF子程序BASBLY组集载荷与节点反力向量PRTREAC输出节点反力PRTDIS输出节点位移ELEMLIB单元库管理程序第四章单刚与载荷的组集，指定约束的处理，单元分析的预处50§4-1子程序PFORM的功能子程序PFORM起到承上启下的作用，为调用与单元运算有关的各种子程序做必要的前后处理。功能参数ISW、ASEM：ISW=3，ASEM=.TRUE.构造单刚，并组集ISW=3，ASEM=.FALSE.构造单元载荷，组集，将指定位移转化为成等效载荷ISW=4计算单元内力或应力，并输出ISW=6构造单元内力，并组集成结构反力ISW=1读入并构造有关的材料常数§4-1子程序PFORM的功能子程序P51SUBROUTINEPFORM(F,LD,UL,XL,S,P,U,X,IX,D,ID,JP,NST,NEN1,ISW)

LOGICALASEMCOMMON/ASEM/ASEM

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL

DIMENSIONUL(NDF,*),U(*),F(NDF,*),S(NST,*),P(*),JP(*),+XL(NDM,*),X(NDM,*),D(20,*),IX(NEN1,*),ID(NDF,*),+LD(NDF,*)

IEL=0

MCT=0！设置打印输出DO110N=1,NE

CALLPZERO(UL,NST*(NST+3))

MA=IX(NEN1,N)！材料号

DO108I=1,NEN

II=IX(I,N)

IF(II.NE.0)GOTO105SUBROUTINEPFORM(F,LD,UL,52DO102J=1,NDF

102LD(J,I)=0

GOTO108

105IID=II*NDF-NDF

NEL=I！确定单元的实际节点数DO106J1=1,NDM

106XL(J1,I)=X(J1,II)！确定单元的节点坐标DO107J=1,NDF

K=ID(J,II)！确定自由度标志数组IF(K.GT.0)UL(J,I)=U(K)！确定单元的活化自由度数INEN=I+NEN

IF(K.LT.0)UL(J,I)=U(NEQ-K)！确定单元的约束位移IF(K.LT.0)UL(J,INEN)=F(J,II)！确定单元的指定位移

IF(ISW.EQ.6)K=IID+J

107LD(J,I)=K！确定自由度标志数组108CONTINUEDO102J=1,NDF

102LD(J,I53IE20=D(20,MA)

IF(IE20.NE.IEL)MCT=0！设置打印输出

IEL=IE20！材料号CALLELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)！确定联系信息IX

IF(ASEM.AND.ISW.EQ.3)CALLADDSTF(LD,JP,S,NST)

130IF(ISW.EQ.6)CALLBASBLY(F,P,LD,NST)

120CONTINUE

IF(.NOT.ASEM.AND.ISW.EQ.3)CALLMODIFY(U,LD,S,UL(1,NEN+1),NST)

IF(.NOT.ASEM.AND.ISW.EQ.3)CALLBASBLY(U,P,LD,NST)

109CONTINUE

110CONTINUE

END

UL：单元位移S：单刚P：单元内力IE20=D(20,MA)

IF(I54§4-2子程序MODIFY完成将指定位移转化成等价内力按活化自由度和被约束自由度，将总刚方程分块：u1

：活化自由度的位移u2

：被约束自由度的位移（包括指定位移）F1

已知F2

：未知§4-2子程序MODIFY完成将指定位移转化成等价内力55

SUBROUTINEMODIFY(U,LD,S,DUL,NS)REALU,DUL,SDIMENSIONU(1),LD(1),S(NS,1),DUL(1)DO110I=1,NSK=LD(I)IF(K.LE.0)GOTO110DO100J=1,NS100U(K)=U(K)-S(I,J)*DUL(J)110CONTINUERETURNENDSUBROUTINEMODIFY(U,LD56§4-3将单刚组集成总刚子程序ADDSTFCALLELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)PFOR35构造了单刚，存放于S中，ADDSTF完成组集PFORM中，LD(NDF,NEN)ADDSTF中,LD(NST)一维数组与二维满存数组的元素对应关系：IJ=JP(I)–I+JIJ§4-3将单刚组集成总刚子程序ADDSTFCALLEL57SUBROUTINEADDSTF(LD,JP,S,NST)CHARACTERY……DIMENSIONLD(*),JP(*),S(NST,NST)IF(Y.EQ.'Y'.OR.Y.EQ.'y')WRITE(*,2000)N,S2000FORMAT('ELEMENTNUMBER',I4,'ELEMENTSTIFFNESS'/(6E12.5))DO200J=1,NSTK=LD(J)IF(K.LE.0)GOTO200L=JP(K)-KDO100I=1,NSTM=LD(I)IF(M.GT.K.OR.M.LE.0)GOTO100M=L+MA(M)=A(M)+S(J,I)100CONTINUE200CONTINUEENDSUBROUTINEADDSTF(LD,JP,58§4-4子程序BASBLY组集载荷与节点反力向量单元自由度指示数组LDISW=3组集载荷向量，按节点活化自由度次序排列ISW=6节点反力向量，按节点排列次序§4-4子程序BASBLY组集载荷与节点反力向量单元自由59130IF(ISW.EQ.6)CALLBASBLY(F,P,LD,NST)PFOR37IF(.NOT.ASEM.AND.ISW.EQ.3)CALLBASBLY(U,P,LD,NST)PFOR40SUBROUTINEBASBLY(B,P,LD,NS)BASB1DIMENSIONB(*),P(*),LD(*)BASB2DO100I=1,NSBASB4II=LD(I)BASB5IF(II.GT.0)B(II)=B(II)+P(I)BASB6100CONTINUEBASB7RETURNBASB8ENDBASB9130IF(ISW.EQ.6)CALLBASBLY60§4-5PRTREAC输出节点反力节点反力两种计算方式：1）由总刚方程求得2）分别计算每个单元内力（单刚乘以单元节点位移），然后组集成总体内力向量因总刚分解后，已不存在了，故选2）§4-5PRTREAC输出节点反力节点反力两种计算方式61SUBROUTINEPRTREAC(R)

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

REALR(NDF,*),RSUM(6)

NNEQ=NDF*NJ

DO50K=1,NDF

RSUM(K)=0.0

50CONTINUE

DO100N=1,NJ,50

J=MIN0(NJ,N+49)

WRITE(6,2000)(K,K=1,NDF)

DO100I=N,J

DO75K=1,NDF

R(K,I)=-R(K,I)

RSUM(K)=RSUM(K)+R(K,I)

75CONTINUESUBROUTINEPRTREAC(R)

62DO76K=1,NDF76IF(ABS(R(K,I)).GT.1.0E-3)GOTO77GOTO10077WRITE(6,2001)I,(R(K,I),K=1,NDF)100CONTINUEWRITE(6,2002)(RSUM(K),K=1,NDF)RETURNFORMAT(/5X,'NODALREACTIONS'/2X,'NODE',6(I8,'DOF'))2001FORMAT(I6,6F12.4)2002FORMAT(3X,'SUM',6F12.4)ENDDO76K=1,NDF63§4-6PRTDIS输出节点位移节点位移按节点次序重新排列注意：1）解得得位移向量为活化自由度，编号小于NEQ2）约束(包括指定位移)得值存放在二维数组F（载荷）中§4-6PRTDIS输出节点位移节点位移按节点次序重新64SUBROUTINEPRTDIS(ID,U,F)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQDIMENSIONID(NDF,*),U(*),F(NDF,*),UL(6)DO102II=1,NJ,50WRITE(6,2000)(K,K=1,NDF)JJ=MIN0(NJ,II+49)DO102N=II,JJDO100I=1,NDFK=ID(I,N)IF(K.LT.0)U(NEQ-K)=F(I,N)IF(K.LT.0)UL(I)=U(NEQ-K)100IF(K.GT.0)UL(I)=U(K)WRITE(6,2001)N,(UL(I),I=1,NDF)102CONTINUE2000FORMAT(//1X,'NODALDISPLACEMENTS'/2X,'NODE',6(I8,'DISP'))2001FORMAT(I6,6E12.4)ENDSUBROUTINEPRTDIS(ID,U,F65§4-7ELEMLIB单元库管理程序FEP2中已有两种单元：BEAM、PLANE若要增加单元，可修改第七句。例如，已编好一个板单元，名为：PLATE（形参），则程序修改如下：……GOTO(1,2,3)IEL……CALLPLATE（D,U,X,IX,S,P,I,J,K,ISW)RETURN……§4-7ELEMLIB单元库管理程序FEP2中已有两种66SUBROUTINEELEMLIB(D,U,X,IX,S,P,I,J,K,ISW)REALP(K),S(K,K),U(1)DIMENSIONIX(1),D(1),X(1)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELIF(IEL.LE.0.OR.IEL.GT.2)GOTO400GOTO(1,2)IEL1CALLBEAM(D,U,X,IX,S,P,I,J,K,ISW)GOTO1002CALLPLANE(D,U,X,IX,S,P,I,J,K,ISW)GOTO100100RETURN400WRITE(*,4000)IELSTOPFORMAT(5X,'**FATALERROR**ELEMENTCLASSNUMBER',I5,'INPUT')ENDSUBROUTINEELEMLIB(D67第五章平面单元设计平面四节点单元的基本公式数值积分形函数及其导数四节点单刚的构造3–9节点平面单元平面单元PLANE子程序第五章平面单元设计平面四节点单元的基本公式68§5-1平面四节点单元的基本公式在FEP2中，提供了3—9节点平面单元，其基本形式是四节点等参元（四边形单元）。子程序：PLANE(D,UL,XL,IX,S,P,NDF,NDM,NST,ISW)

PLANE中包含的子程序PGAUSS(L,LINT,SG,TG,WG)SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.)SHAP2(SS,TT,SHP,IX,NEL)PSTRES(SIG,SIG(4),SIG(5),SIG(6))§5-1平面四节点单元的基本公式在FEP269§5-2数值积分

SUBROUTINEPGAUSS(L,LINT,R,Z,W)COMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONLR(9),LZ(9),LW(9),R(3),Z(3),W(3),G4(4),H4(4)DATALR/-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/LINT=L*LGOTO(1,2,3,4),L1R(1)=0.0E0Z(1)=0.0E0W(1)=4.0E0IF(NEL.EQ.3)THENZ(1)=-1.0E0/3.0E0W(1)=3.0E0ENDIFRETURN§5-2数值积分SUBROUTINEPG702G=1.0E0/SQRT(3.0E0)DO21I=1,4R(I)=G*LR(I)Z(I)=G*LZ(I)21W(I)=1.0E0RETURN3G=SQRT(0.6E0)H=1.0E0/81.0E0DO31I=1,9R(I)=G*LR(I)Z(I)=G*LZ(I)31W(I)=H*LW(I)RETURN2G=1.0E0/SQRT(3.0E0)714G=SQRT(4.8E0)H=SQRT(30.0E0)/36.0E0G4(1)=SQRT((3.0E0+G)/7.0E0)G4(4)=-G4(1)G4(2)=SQRT((3.0E0-G)/7.0E0)G4(3)=-G4(2)H4(1)=0.50E0–H；H4(2)=0.50E0+HH4(3)=0.50E0+H；H4(4)=0.50E0-HI=0DO41J=1,4DO41K=1,4I=I+1R(I)=G4(K)Z(I)=G4(J)W(I)=H4(J)*H4(K)41CONTINUEEND4G=SQRT(4.8E0)72§5-3形函数及其导数SUBROUTINESHAPE(SS,TT,X,SHP,XSJ,NDM,NEL,IX,FLG)DIMENSIONSHP(3,*),X(NDM,1),S(4),T(4),XS(2,2),SX(2,2),IX(*)DATAS/-0.5E0,0.5E0,0.5E0,-0.5E0/,T/-0.5E0,-0.5E0,0.5E0,0.5E0/DO100I=1,4SHP(3,I)=(0.5e0+S(I)*SS)*(0.5e0+T(I)*TT)SHP(1,I)=S(I)*(0.5e0+T(I)*TT)100SHP(2,I)=T(I)*(0.5e0+S(I)*SS)IF(NEL.GE.4)GOTO120C...FORMTRIANGLEBYADDINGTHIRDANDFOURTHTOGETHERDO110I=1,3110SHP(I,3)=SHP(I,3)+SHP(I,4)！3节点平面单元

C....ADDQUADRATICTERMSIFNECESSARY120IF(NEL.GT.4)CALLSHAP2(SS,TT,SHP,IX,NEL)§5-3形函数及其导数SUBROUTINESHAPE73C....CONSTRUCTJACOBIANANDITSINVERSEDO130I=1,NDMDO130J=1,2XS(I,J)=0.0DO130K=1,NEL130XS(I,J)=XS(I,J)+X(I,K)*SHP(J,K)XSJ=XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1)IF(FLG)RETURNSX(1,1)=XS(2,2)/XSJSX(2,2)=XS(1,1)/XSJSX(1,2)=-XS(1,2)/XSJSX(2,1)=-XS(2,1)/XSJC....FORMGLOBALDERIVATIVESDO140I=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)140SHP(1,I)=TPENDC....CONSTRUCTJACOBIANANDI74§5-4四节点单刚的构造比较三种不同的算法，选运算量最少的一种。1）按矩阵直接运算2）预先算出DBj

，再与BiT相乘，展开3）按矩阵运算全部展开计算：BiT

DBj§5-4四节点单刚的构造比较三种不同的算法，选运算量最少753L=D(5)IF(L*L.NE.LINT)CALLPGAUSS(L,LINT,SG,TG,WG)DO320L=1,LINTCALLSHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.)XSJ=XSJ*WG(L)*D(7)J1=1DO320J=1,NELW11=SHP(1,J)*XSJW12=SHP(2,J)*XSJK1=J1DO310K=J,NELS(J1,K1)=S(J1,K1)+W11*SHP(1,K)S(J1,K1+1)=S(J1,K1+1)+W11*SHP(2,K)S(J1+1,K1)=S(J1+1,K1)+W12*SHP(1,K)S(J1+1,K1+1)=S(J1+1,K1+1)+W12*SHP(2,K)310K1=K1+NDF320J1=J1+NDF3L=D(5)76§5-53–9节点平面单元

3节点平面单元SHAPE12~135–9节点平面单元§5-53–9节点平面单元3节点平面单元SHA77SUBROUTINESHAP2(S,T,SHP,IX,NEL)DIMENSIONIX(*),SHP(3,*)S2=(1.0E0-S*S)/2.0E0T2=(1.0E0-T*T)/2.0E0DO100I=5,9DO100J=1,3100SHP(J,I)=0.0IF(IX(5).EQ.0)GOTO101SHP(1,5)=-S*(1.0E0-T)SHP(2,5)=-S2SHP(3,5)=S2*(1.0E0-T)101IF(NEL.LT.6)GOTO107IF(IX(6).EQ.0)GOTO102SHP(1,6)=T2SHP(2,6)=-T*(1.0E0+S)SHP(3,6)=T2*(1.0E0+S)SUBROUTINESHAP2(S,T,SHP,IX,N78102IF(NEL.LT.7)GOTO107SHAP19IF(IX(7).EQ.0)GOTO103SHAP20SHP(1,7)=-S*(1.0E0+T)SHAP21SHP(2,7)=S2SHAP22SHP(3,7)=S2*(1.0E0+T)SHAP23103IF(NEL.LT.8)GOTO107SHAP24IF(IX(8).EQ.0)GOTO104SHAP25SHP(1,8)=-T2SHAP26SHP(2,8)=-T*(1.0E0-S)SHAP27SHP(3,8)=T2*(1.0E0-S)SHAP28C....INTERIORNODE(LAGRANGIAN)104IF(NEL.LT.9)GOTO107SHAP30IF(IX(9).EQ.0)GOTO107SHAP31SHP(1,9)=-4.0E0*S*T2SHAP32SHP(2,9)=-4.0E0*T*S2SHAP33SHP(3,9)=4.0E0*S2*T2102IF(NEL.LT.7)GOTO10779C....CORRECTEDGENODESFORINTERIORNODEDO106J=1,3SHAP36DO105I=1,4SHAP37105SHP(J,I)=SHP(J,I)-0.25E0*SHP(J,9)SHAP38DO106I=5,8