有限元程序设计

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

读入和生成数据，形成有限元网格，为有限元矩阵计算作准备。（数据或图形输入）标题和控制信息计算存贮分配节点坐标和约束信息单元信息材料信息载荷信息第六页，共八十五页，2022年，8月28日二、有限元矩阵的计算、组集和求解—处理器计算插值函数矩阵[N

]、几何矩阵[B

]、Jacob矩阵[J

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

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

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

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

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

—PFORM(6)计算结构反力第十九页，共八十五页，2022年，8月28日§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的前四个字符作为输入输出文件名第二十页，共八十五页，2022年，8月28日COMMON/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)

第二十一页，共八十五页，2022年，8月28日

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

第二十二页，共八十五页，2022年，8月28日二、定义数组FEP2中设置了两个数组M(1600)和A(16380)，数组M是作为存放坐标、单元信息、载荷、位移等，数组A则是存放刚度矩阵用的，其上界与FORTRAN（编译器）版本和计算机内存有关，确定了程序的计算规模。这两个数组也可合并成一个。第二十三页，共八十五页，2022年，8月28日§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'/第二十四页，共八十五页，2022年，8月28日FO：输出文件名ERR：逻辑变量，作出错信息控制（.TRUE.为错）ASEM：逻辑变量，作为判断是否要组集总刚数组A：存放总刚的元素数组M：存放单刚、单元信息、节点坐标等NDF—最大自由度数NDM—空间维数NEN—单元的最多节点数NJ—节点总数NE—总单元数NEQ—总方程数NUMMAT—材料类型数FD,FORC：文字变量，作为输出的标题用P：文字变量，判断是否要输出单刚第二十五页，共八十五页，2022年，8月28日二、输入控制信息、确定主要数组的地址（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的首址

……第二十六页，共八十五页，2022年，8月28日READ(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+NNEQ第二十七页，共八十五页，2022年，8月28日§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’第二十八页，共八十五页，2022年，8月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)第二十九页，共八十五页，2022年，8月28日LI=(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……第三十页，共八十五页，2022年，8月28日3）单元信息（一组）L，LX，LK，（IDL(K)，K=1，NEN）单元号单元号增量材料号节点号0，0，0，……结束由ELDA子程序生成单元信息，IX(NEN1,NE)二维数组。NEN=NEN+1NE—单元总数第三十一页，共八十五页，2022年，8月28日4）材料信息（一组或多组）对于每一类材料，要求输入两个记录：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处理成等效节点载荷第三十二页，共八十五页，2022年，8月28日对于平面元，输入6个常数：E,XNU,TH，L，K，ITH：单元厚度L：沿一个方向的高斯积分点数（2，3，4）K：沿一个方向应力输出的高斯点数I：=0平面应变0平面应力PMESH调用ELEMLIB，再调用BEAMorPLANE第三十三页，共八十五页，2022年，8月28日SUBROUTINEPMESH(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)

第三十四页，共八十五页，2022年，8月28日1CALLGENVEC(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)END第三十五页，共八十五页，2022年，8月28日5）约束信息（包括指定位移）（一组）PMESH调用BOUN输入：N，NG，（IDL(I)，I=1，NDF）

节点号节点号增量坐标或载荷0，0，0，……（平面元4个0，梁元5个0）IDL(I)：约束信息，被约束的自由度给非零值（1或-1），其它为0。1：该自由度被约束，但不继续生成-1：该自由度被约束，要求继续生成第三十六页，共八十五页，2022年，8月28日6）载荷信息载荷输入仍由GENVEC完成，维数：NDF

对于指定位移，生成方式与节点载荷相同，5）输入约束信息，6）输入其值。第三十七页，共八十五页，2022年，8月28日第三章总刚度矩阵的变带宽存贮与求解总刚度矩阵的变带宽存贮与对角线地址LDLT三角分解自由度指示矩阵ID与对角线地址矩阵第三十八页，共八十五页，2022年，8月28日§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第三十九页，共八十五页，2022年，8月28日SUBROUTINEPROFIL(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

50CONTINUE第四十页，共八十五页，2022年，8月28日C.....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))100CONTINUE200CONTINUE300CONTINUECONTINUE500CONTINUE第四十一页，共八十五页，2022年，8月28日C.....COMPUTEDIAGONALPOINTERSFORPROFILNK=1JP(1)=1IF(NEQ.EQ.1)RETURNDO600N=2,NEQ600JP(N)=JP(N)+JP(N-1)+1NK=JP(NEQ)WRITE(*,2000)NK2000FORMAT('THEMAXIMUMSTORAGEOFGLOBALSTIFFNESS=',I6/)RETURNEND第四十二页，共八十五页，2022年，8月28日

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))第四十三页，共八十五页，2022年，8月28日

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)100CONTINUE200CONTINUERETURNEND第四十四页，共八十五页，2022年，8月28日§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(*)第四十五页，共八十五页，2022年，8月28日Ax=b=>LDLTx=b令LTx=y=>LDy=b

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

SUBROUTINEFORBACK(A,B,JP,N)DIMENSIONA(*),B(*),JP(*)第四十六页，共八十五页，2022年，8月28日§3-3自由度指示矩阵ID与对角线地址矩阵ID(NDF,NJ)：非零值（1或-1）为被约束的自由度，0为活化自由度。在PROFIL中，改造数组ID，原来的零元素（活化自由度）用正数计，原来的非零元素（非活化自由度）用负数计，全部零元素（活化）的个数记为NEQ（总刚矩阵的阶数）。第四十七页，共八十五页，2022年，8月28日计算刚度矩阵的带宽的方案比较：1）以各节点自由度为研究对象，利用单元信息IX，找出与该点联系的节点，然后比较节点自由度之最大者。运算量：NJ*NDF*NE*NEN*NDF2）逐一研究各单元，分别取出各节点的各个自由度与该单元中其它节点的自由度比较，差值最大者再与所研究的自由度之带宽值（初值为零）比较，取其中最大者代替原有的带宽值。运算量：NE*NEN*NDF*NEN*NDF两者运算量比：NJ:NE选方案2带宽=自由度之最大差值+1暂存放于JP中第四十八页，共八十五页，2022年，8月28日若已知对角线地址数组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)第四十九页，共八十五页，2022年，8月28日第四章单刚与载荷的组集，指定约束的处理，单元分析的预处理子程序PFORM的功能子程序MODIFY完成将指定位移转化成等价内力将单刚组集成总刚子程序ADDSTF子程序BASBLY组集载荷与节点反力向量PRTREAC输出节点反力PRTDIS输出节点位移ELEMLIB单元库管理程序第五十页，共八十五页，2022年，8月28日§4-1子程序PFORM的功能子程序PFORM起到承上启下的作用，为调用与单元运算有关的各种子程序做必要的前后处理。功能参数ISW、ASEM：ISW=3，ASEM=.TRUE.构造单刚，并组集ISW=3，ASEM=.FALSE.构造单元载荷，组集，将指定位移转化为成等效载荷ISW=4计算单元内力或应力，并输出ISW=6构造单元内力，并组集成结构反力ISW=1读入并构造有关的材料常数第五十一页，共八十五页，2022年，8月28日SUBROUTINEPFORM(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)GOTO105第五十二页，共八十五页，2022年，8月28日DO102J=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！确定自由度标志数组108CONTINUE第五十三页，共八十五页，2022年，8月28日IE20=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：单元内力第五十四页，共八十五页，2022年，8月28日§4-2子程序MODIFY完成将指定位移转化成等价内力按活化自由度和被约束自由度，将总刚方程分块：u1

：活化自由度的位移u2

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

已知F2

：未知第五十五页，共八十五页，2022年，8月28日

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)110CONTINUERETURNEND第五十六页，共八十五页，2022年，8月28日§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第五十七页，共八十五页，2022年，8月28日SUBROUTINEADDSTF(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)100CONTINUE200CONTINUEEND第五十八页，共八十五页，2022年，8月28日§4-4子程序BASBLY组集载荷与节点反力向量单元自由度指示数组LDISW=3组集载荷向量，按节点活化自由度次序排列ISW=6节点反力向量，按节点排列次序第五十九页，共八十五页，2022年，8月28日130IF(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)BASB6100CONTINUEBASB7RETURNBASB8ENDBASB9第六十页，共八十五页，2022年，8月28日§4-5PRTREAC输出节点反力节点反力两种计算方式：1）由总刚方程求得2）分别计算每个单元内力（单刚乘以单元节点位移），然后组集成总体内力向量因总刚分解后，已不存在了，故选2）第六十一页，共八十五页，2022年，8月28日SUBROUTINEPRTREAC(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)

75CONTINUE第六十二页，共八十五页，2022年，8月28日DO76K=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)END第六十三页，共八十五页，2022年，8月28日§4-6PRTDIS输出节点位移节点位移按节点次序重新排列注意：1）解得得位移向量为活化自由度，编号小于NEQ2）约束(包括指定位移)得值存放在二维数组F（载荷）中第六十四页，共八十五页，2022年，8月28日SUBROUTINEPRTDIS(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)END第六十五页，共八十五页，2022年，8月28日§4-7ELEMLIB单元库管理程序FEP2中已有两种单元：BEAM、PLANE若要增加单元，可修改第七句。例如，已编好一个板单元，名为：PLATE（形参），则程序修改如下：……GOTO(1,2,3)IEL……CALLPLATE（D,U,X,IX,S,P,I,J,K,ISW)RETURN……第六十六页，共八十五页，2022年，8月28日SUBROUTINEELEMLIB(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')END第六十七页，共八十五页，2022年，8月28日第五章平面单元设计平面四节点单元的基本公式数值积分形函数及其导数四节点单刚的构造3–9节点平面单元平面单元PLANE子程序第六十八页，共八十五页，2022年，8月28日§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))第六十九页，共八十五页，2022年，8月28日§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第七十页，共八十五页，2022年，8月28日2G=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)RETURN第七十一页，共八十五页，2022年，8月28日4G=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)41CONTINUEEND第七十二页，共八十五页，2022年，8月28日§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)第七十三页，共八十五页，2022年，8月28日C....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)=TPEND第七十四页，共八十五页，2022年，8月28日§5-4四节点单刚的构造比较三种不同的算法，选运算量最少的一种。1）按矩阵直接运算2）预先算出DBj

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

DBj第七十五页，共八十五页，2022年，8月28日3L=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+NDF第七十六页，共八十五页，2022年，8月28日§5-53–9节点平面单元

3节点平面单元SHAPE12~135–9节点平面单元第七十七页，共八十五页，2022年，8月28日SUBROUTINESHAP2(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)第七十八页，共八十五页，2022年，8月28日102IF(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*T2第七十九页，共八十五页，2022年，8月28日C....CORRECTEDGENODESFORINTERIORNODEDO106J=1,3SHAP36DO105I=1,4SHAP37105SHP(J,I)=SHP(J,I)-0.25E0*SHP(J,9)SHAP38DO106I=5,8