有限元程序设计.ppt

有限元程序设计有限元程序设计教学程序：教学程序：FEP1 李明瑞编李明瑞编 FEP2 李明瑞编李明瑞编 ZFEP第一章第一章 绪论绪论 有限元程序的基本内容有限元程序的基本内容 有限元求解器有限元求解器1.1 有限元程序的基本内容有限元程序的基本内容有限元法的解题步骤有限元法的解题步骤结构离散化为结构离散化为有限元网格有限元网格计算单元刚度计算单元刚度矩阵矩阵引入约束条件引入约束条件求求应力等应力等解解方程组方程组建立总刚度矩建立总刚度矩阵及外载阵及外载有限元程序的基本内容：有限元程序的基本内容：数据输入阶数据输入阶段段有限元矩阵的计有限元矩阵的计算、组集和求解算、组集和求解数据输出数据输出阶段阶段前前处理处理处理器处理器后后处理处理一、一、数据输入阶段数据输入阶段前处理 读读入入和和生生成成数数据据，形形成成有有限限元元网网格格，为为有有限限元元矩矩阵计算作准备。（数据或图形输入）阵计算作准备。（数据或图形输入）1.1.标题和控制信息标题和控制信息2.2.计算存贮分配计算存贮分配3.3.节点坐标和约束信息节点坐标和约束信息4.4.单元信息单元信息5.5.材料信息材料信息6.6.载荷信息载荷信息二、二、有限元矩阵的计算、组集和求解有限元矩阵的计算、组集和求解处理器处理器计计算算插插值值函函数数矩矩阵阵 N N 、几几何何矩矩阵阵 B B 、JacobJacob矩矩阵阵 J J 在在高斯点进行数值积分，求得单刚、单元载荷高斯点进行数值积分，求得单刚、单元载荷组集成总刚、总载组集成总刚、总载 求解求解三、三、数据输出阶段数据输出阶段后后处理输出结果（场变量和场变量导数等），打印结果。输出结果（场变量和场变量导数等），打印结果。场变量包括：位移、温度、流场的速度势等；场变量包括：位移、温度、流场的速度势等；场变量导数包括：应力、应变、热流、流场速度势等场变量导数包括：应力、应变、热流、流场速度势等输出形式：输出形式：数据输出和图形输出数据输出和图形输出1-2 有限元求解器有限元求解器带宽法：只存储总刚矩阵的半带宽内的元素。带宽法：只存储总刚矩阵的半带宽内的元素。等带宽法：每行的带宽相等，常采用二维数组存储等带宽法：每行的带宽相等，常采用二维数组存储变带宽法：每行的带宽不等，常采用一维数组存储变带宽法：每行的带宽不等，常采用一维数组存储一、一、带宽法带宽法变带宽分块求解思想：变带宽分块求解思想：一个变量的消元只影响刚度矩阵的一部分元一个变量的消元只影响刚度矩阵的一部分元素，较大的节点分量方程组可以分成较小的局部素，较大的节点分量方程组可以分成较小的局部系数矩阵按节点逐步求解。系数矩阵按节点逐步求解。波波前前法法是是另另一一种种分分块块储储存存的的变变带带宽宽法法，其其消消元元次次序序是按单元编号进行，边组集边消元。是按单元编号进行，边组集边消元。调调入入内内存存的的单单元元所所保保留留的的节节点点称称为为波波前前节节点点，所所消消去去的的节节点点称称为为消消元元节节点点，消消元元节节点点是是与与以以后后调调入入的的单元没有联系的节点，即该节点的方程已组集完。单元没有联系的节点，即该节点的方程已组集完。波波前前法法的的回回代代按按消消元元节节点点的的反反序序进进行行。对对内内存存的的需求取决与最大波前宽。需求取决与最大波前宽。二、二、波前法波前法V-Fortran 1.Fortran语言的基本格式语言的基本格式2.变量变量3.基本语句基本语句4.子例子程序子例子程序5.函数子程序函数子程序 6.其它功能、模块其它功能、模块第二章第二章 FEP2 程序的总体设计与输入程序的总体设计与输入数据数据 FEP2 的设计任务的设计任务 FEP2 的结构设计的结构设计 FEP2 的主程序的主程序 FEP2 的主控程序的主控程序 FEP2 的数据输入格式与程序实现的数据输入格式与程序实现2-1 FEP2 的设计任务的设计任务1.FEP2的简要题目说明的简要题目说明 FEP2是一个具有通用性的教学程序，可用于计算是一个具有通用性的教学程序，可用于计算一般的线性静力学问题。已设计了平面梁单元与平面一般的线性静力学问题。已设计了平面梁单元与平面3-9节点元两种单元，但留下接口。节点元两种单元，但留下接口。2.支持软件与硬件支持软件与硬件 FORTRAN77以上编译器、各种微机以上编译器、各种微机3.FEP2的的功能功能3.FEP2的功能的功能1）输入文件名由用户自由定义，但限制为输入文件名由用户自由定义，但限制为4个字符，输个字符，输入文件扩展名为入文件扩展名为“DAT”，输出文件扩展名为输出文件扩展名为“OUT”。2）节点坐标、单元信息等具有线性自动生成功能。节点坐标、单元信息等具有线性自动生成功能。3）可以处理多种工况、多种类型单元组合结构问题。）可以处理多种工况、多种类型单元组合结构问题。4）可以处理固定约束和指定位移。）可以处理固定约束和指定位移。5）采用变带宽存储、三角分解法求解刚度方程。）采用变带宽存储、三角分解法求解刚度方程。6）为多种单元留下接口。）为多种单元留下接口。2-2 FEP2 的的结构设计结构设计FEP2：由形式主程序、主控程序、功能模块组成。由形式主程序、主控程序、功能模块组成。主程序的主要功能：定义输入、输出文件名，调用主控主程序的主要功能：定义输入、输出文件名，调用主控程序程序PCONTR主控程序主控程序PCONTR的主要功能模块：的主要功能模块：1）内存分配内存分配 2）网格生成）网格生成3）变带宽存储设置变带宽存储设置 4）单刚的形成与组集）单刚的形成与组集5）刚度矩阵的分解）刚度矩阵的分解 6）形成节点载荷）形成节点载荷7）形成与组集单元载荷）形成与组集单元载荷 8）求解位移）求解位移9）求单元应力）求单元应力 10）求结构反力）求结构反力FEP2 程序框图程序框图FEP2（主程序）主程序）PCONTR(主控程序主控程序)SETMEM 检查内存检查内存 PROFIL 形成变带宽刚度矩阵地址形成变带宽刚度矩阵地址 PFORM(3)形成单刚并组集形成单刚并组集 LDLT 总刚的三角分解总刚的三角分解 GENVEC 形成节点载荷形成节点载荷 PFOM(3)构造并组集单元载荷构造并组集单元载荷 FORBACK 前消回代求出位移前消回代求出位移 PFORM(4)计算单元应力计算单元应力 PFORM(6)计算结构反力计算结构反力2-3 FEP2 的主程序的主程序一、文件名一、文件名 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的前的前四个字符作为输入输出文件名四个字符作为输入输出文件名 COMMON/PSIZE/MAXM,MAXA CHARACTER*8 FI,FO CHARACTER NAMINP(8),NAMOUT(8),HEAD*50 COMMON/HEAD/HEAD1 EQUIVALENCE(FI,NAMINP(1),(FO,NAMOUT(1)DATA FI,FO/.DAT,.OUT/WRITE(*,(A)INPUT FILE NAME(4 LETTERS ONLY):READ(*,(4A1)(NAMINP(I),I=1,4)DO 10 I=1,4 10 NAMOUT(I)=NAMINP(I)OPEN(6,FILE=FO)OPEN(5,FILE=FI)MAXM=16000 MAXA=16380 CALL PCONTR(FO)END 二、定义数组二、定义数组 FEP2 中设置了两个数组中设置了两个数组M(1600)和和 A(16380)，数数组组M是作为存放坐标、单元信息、载荷、位移等，数组是作为存放坐标、单元信息、载荷、位移等，数组 A 则是存放刚度矩阵用的，其上界与则是存放刚度矩阵用的，其上界与FORTRAN（编译编译器）版本和计算机内存有关，确定了程序的计算规模。器）版本和计算机内存有关，确定了程序的计算规模。这两个数组也可合并成一个。这两个数组也可合并成一个。2-4 FEP2 的主控程序的主控程序 PCONTR PCONTR：实质上的主程序。实质上的主程序。分以下几大段：分以下几大段：一、程序头（前一、程序头（前1212语句）语句）SUBROUTINE PCONTR(FO)CHARACTER*8 FO LOGICAL ERR,ASEM COMMON/M/M(6000)COMMON/A/A(16380)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON/PSIZE/MAXM,MAXA COMMON/MDATA/NUMMAT CHARACTER FD*25,FORC*6,P COMMON/PRINT/P COMMON/ASEM/ASEM DATA FORC/FORC/FD/NODAL FORCE/DISPL/FO：输出文件名输出文件名 ERR：逻辑变量，作出错信息控制（逻辑变量，作出错信息控制（.TRUE.为错）为错）ASEM：逻辑变量，作为判断是否要组集总刚逻辑变量，作为判断是否要组集总刚 数组数组A：存放总刚的元素存放总刚的元素 数组数组M：存放单刚、单元信息、节点坐标等存放单刚、单元信息、节点坐标等 NDF 最大自由度数最大自由度数 NDM 空间维数空间维数 NEN 单元的最多节点数单元的最多节点数 NJ 节点总数节点总数 NE 总单元数总单元数 NEQ 总方程数总方程数NUMMAT 材料类型数材料类型数 FD,FORC：文字变量，作为输出的标题用文字变量，作为输出的标题用 P：文字变量，判断是否要输出单刚文字变量，判断是否要输出单刚二、输入控制信息、确定主要数组的地址（二、输入控制信息、确定主要数组的地址（14341434语句）语句）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的首址的首址 READ(5,*)NJ,NE,NUMMAT,NDF,NDM,NEN WRITE(*,2000)NJ,NE,NUMMAT,NDF,NDM,NEN WRITE(6,2000)NJ,NE,NUMMAT,NDF,NDM,NEN NEN1=NEN+1 NNEQ=NJ*NDF NST=NDF*NEN NXC=1 NIX=NXC+NJ*NDM NS=NUL+NST*2 NP=NS+NST*NST NF=NP+NST NU=NF+NNEQ NJP=NU+NNEQ NEND=NJP+NNEQ 2-5 FEP2的数据输入格式与程序实现的数据输入格式与程序实现FEP2的数据信息（自由格式）：的数据信息（自由格式）：1）控制信息：）控制信息：NJ，NE，NUMMAT，NDM，NEN2）坐标信息（一组）坐标信息（一组）N，NG，（XL(I)，I=1，NDM）节点号节点号 节点号增量节点号增量 坐标或载荷坐标或载荷 由由GENVEC生成节点坐标或载荷向量，通过生成节点坐标或载荷向量，通过NDM,CD，FC 进行识别生成的是节点坐标，还是载荷向量。进行识别生成的是节点坐标，还是载荷向量。CD：NODAL COORDINATES FC：COORCD：NODAL FORCE/DISPL FC：FORC SUBROUTINE GENVEC(NDM,X,CD,FC,NJ,ERR)CHARACTER CD*18,FC*6 LOGICAL ERR REAL X(NDM,*),XL(6)ERR=.FALSE.N=0 NG=0102 L=N LG=NG READ(5,*,ERR=999)N,NG,(XL(I),I=1,NDM)101 IF(N.LE.0.OR.N.GT.NJ)GO TO 109 DO 103 I=1,NDM103 X(I,N)=XL(I)IF(LG)104,102,104104 LG=ISIGN(LG,N-L)LI=(IABS(N-L+LG)-1)/IABS(LG)DO 105 I=1,NDM 105 XL(I)=(X(I,N)-X(I,L)/LI106 L=L+LG IF(N-L)*LG.LE.0)GO TO 102 IF(L.LE.0.OR.L.GT.NJ)GO TO 108 DO 107 I=1,NDM LMLG=L-LG107 X(I,L)=X(I,LMLG)+XL(I)GO TO 106108 WRITE(6,3000)L,CD WRITE(*,3000)L,CD ERR=.TRUE.GO TO 102109 CONTINUE 3）单元信息（一组）单元信息（一组）L，LX，LK，（IDL(K)，K=1，NEN）单元号单元号 单元号增量单元号增量 材料号材料号 节点号节点号 0，0，0，结束结束由由ELDA子程序生成单元信息，子程序生成单元信息，IX(NEN1,NE)二维数组。二维数组。NEN=NEN+1 NE 单元总数单元总数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，非铰为，非铰为0 N2：右端铰为右端铰为1，非铰为，非铰为0 Ni：单元载荷类型（单元载荷类型（0，1，2，3，4）（0为无载荷）为无载荷）单元载荷由单元载荷由 EQLOAD 处理成等效节点载荷处理成等效节点载荷 对于平面元，输入对于平面元，输入6个常数：个常数：E,XNU,TH，L，K，I TH：单元厚度单元厚度 L：沿一个方向的高斯积分点数沿一个方向的高斯积分点数（2，3，4）K：沿一个方向应力输出的高斯点数沿一个方向应力输出的高斯点数 I：=0 平面应变平面应变 0 平面应力平面应力 PMESH 调用调用 ELEMLIB，再调用再调用 BEAM or PLANE SUBROUTINE PMESH(XC,IX,ID,D,NEN1)CHARACTER COOR*6,XX*18 COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON/MDATA/NUMMAT COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL DIMENSION XC(NDM,*),IX(NEN1,*),ID(NDF,*),D(20,*)LOGICAL ERR DATA XX/NODAL COORDINATES/COOR/COOR/1 CALL GENVEC(NDM,XC,XX,COOR,NJ,ERR)1 CALL GENVEC(NDM,XC,XX,COOR,NJ,ERR)IF(ERR)WRITE(*,3000)IF(ERR)STOP 2 CALL ELDAT(IX,NEN1)3 DO 304 N=1,NUMMAT WRITE(6,2000)MA,IEL WRITE(*,2000)MA,IEL IF(MA.EQ.0)GO TO 4 CALL PZERO(D(1,MA),20)D(20,MA)=IEL CALL LEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,NST,1)304 CONTINUE 2000 FORMAT(/MATERIAL TYPE,I4,ELEMENT TYPE,I4/)4 CALL BOUN(ID)END 5）约束信息（包括指定位移）（一组）约束信息（包括指定位移）（一组）PMESH 调用调用 BOUN 输入：输入：N，NG，（IDL(I)，I=1，NDF）节点号节点号 节点号增量节点号增量 坐标或载荷坐标或载荷 0，0，0，（平面元（平面元4个个0，梁元，梁元5个个0）IDL(I)：约束信息，被约束的自由度给非零值（约束信息，被约束的自由度给非零值（1或或-1），），其它为其它为0。1：该自由度被约束，但不继续生成该自由度被约束，但不继续生成-1：该自由度被约束，要求继续生成：该自由度被约束，要求继续生成6）载荷信息）载荷信息 载荷输入仍由载荷输入仍由 GENVEC 完成，维数：完成，维数：NDF 对于指定位移，生成方式与节点载荷相同，对于指定位移，生成方式与节点载荷相同，5）输入）输入约束信息，约束信息，6）输入其值。）输入其值。第三章第三章 总刚度矩阵的变带宽存贮与求总刚度矩阵的变带宽存贮与求解解总刚度矩阵的变带宽存贮与对角线地址总刚度矩阵的变带宽存贮与对角线地址LDLT 三角分解三角分解自由度指示矩阵自由度指示矩阵ID与对角线地址矩阵与对角线地址矩阵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 SUBROUTINE PROFIL(JP,ID,IX,NEN1,NK)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ DIMENSION JP(1),ID(NDF,1),IX(NEN1,*),IDL(18)NEQ=0 NAD=0 DO 50 N=1,NJ DO 40 I=1,NDF J=ID(I,N)IF(J)30,20,3020 NEQ=NEQ+1 ID(I,N)=NEQ JP(NEQ)=0 GO TO 4030 NAD=NAD-1 ID(I,N)=NAD40 CONTINUE50 CONTINUEC.COMPUTE ROW BANDWIDTH DO 500 N=1,NE DO 400 I=1,NEN II=IX(I,N)IF(II.EQ.0)GO TO 400 DO 300 K=1,NDF KK=ID(K,II)IF(KK.LE.0)GO TO 300 DO 200 J=1,NEN JJ=IX(J,N)IF(JJ.EQ.0)GO TO 200 DO 100 L=1,NDF LL=ID(L,JJ)IF(LL.LE.0)GO TO 100 M=MAX0(KK,LL)JP(M)=MAX0(JP(M),IABS(KK-LL)100 CONTINUE200 CONTINUE300 CONTINUE400CONTINUE500 CONTINUEC.COMPUTE DIAGONAL POINTERS FOR PROFIL NK=1 JP(1)=1 IF(NEQ.EQ.1)RETURN DO 600 N=2,NEQ600 JP(N)=JP(N)+JP(N-1)+1 NK=JP(NEQ)WRITE(*,2000)NK2000 FORMAT(THE MAXIMUM STORAGE OF GLOBAL STIFFNESS=,I6/)RETURN END SUBROUTINE ADDSTF(LD,JP,S,NST)CHARACTER Y COMMON /PRINT/Y COMMON /A/A(16380)COMMON /CDATA/NDF,NDM,NEN,NEJ,NE,NEQ COMMON /ELDATA/LM,N,MA,MCT,IEL,NEL DIMENSION LD(*),JP(*),S(NST,NST)IF(Y.EQ.Y.OR.Y.EQ.y)WRITE(*,2000)N,S2000 FORMAT(ELEMENT NUMBER ,I4,ELEMENT STIFFNESS/(6E12.5)DO 200 J=1,NST K=LD(J)IF(K.LE.0)GO TO 200 L=JP(K)-K DO 100 I=1,NST M=LD(I)IF(M.GT.K.OR.M.LE.0)GO TO 100 M=L+M A(M)=A(M)+S(J,I)100 CONTINUE200 CONTINUE RETURN END3-2 LDLT 三角分解三角分解求解方程：求解方程：Ax=b A 对称、正定矩阵对称、正定矩阵 =Ux=y U 单位上三角矩阵单位上三角矩阵对对A进行三角分解进行三角分解 =A=LLTL 下三角矩阵，通过比较系数法确定，需要开方运算下三角矩阵，通过比较系数法确定，需要开方运算LDLT三角分解三角分解 =A=LDLTL 下三角矩阵，下三角矩阵，D 对角矩阵对角矩阵i=1,n;j=2,iLDLT(A,JP,N)DIMENSION A(*),JP(*)Ax=b =LDLT x=b 令令 LT x=y =LD y=b 前消前消公式：公式：i=2,n回代公式：i=n-1,1 SUBROUTINE FORBACK(A,B,JP,N)DIMENSION A(*),B(*),JP(*)3-3 自由度指示矩阵自由度指示矩阵ID与对角线地址矩阵与对角线地址矩阵ID(NDF,NJ)：非零值（非零值（1或或-1）为被约束的自由度，）为被约束的自由度，0为活化自由度。为活化自由度。在在PROFIL 中，改造数组中，改造数组ID，原来的零元素（活化原来的零元素（活化自由度）用正数计，原来的非零元素（非活化自由度）自由度）用正数计，原来的非零元素（非活化自由度）用负数计，全部零元素（活化）的个数记为用负数计，全部零元素（活化）的个数记为NEQ（总刚总刚矩阵的阶数）。矩阵的阶数）。计算刚度矩阵的带宽的方案比较：计算刚度矩阵的带宽的方案比较：1）以各节点自由度为研究对象，利用单元信息）以各节点自由度为研究对象，利用单元信息IX，找找出与该点联系的节点，然后比较节点自由度之最大者。出与该点联系的节点，然后比较节点自由度之最大者。运算量：运算量：NJ*NDF*NE*NEN*NDF2）逐一研究各单元，分别取出各节点的各个自由度与）逐一研究各单元，分别取出各节点的各个自由度与该单元中其它节点的自由度比较，差值最大者再与所研该单元中其它节点的自由度比较，差值最大者再与所研究的自由度之带宽值（初值为零）比较，取其中最大者究的自由度之带宽值（初值为零）比较，取其中最大者代替原有的带宽值。代替原有的带宽值。运算量：运算量：NE*NEN*NDF*NEN*NDF两者运算量比：两者运算量比：NJ:NE 选方案选方案 2 带宽带宽=自由度之最大差值自由度之最大差值+1 暂存放于暂存放于 JP 中中 若已知对角线地址数组若已知对角线地址数组 JP，则，则第第 I 行的带宽：行的带宽：JP(1)=1,JP(I)JP(I 1)若若 JP 存放的是带宽存放的是带宽,则对角线地址数组则对角线地址数组：JP(1)=JP(1)=1 JP(2)=JP(2)+JP(1)JP(I+1)=JP(I+1)+JP(I)总刚矩阵总刚矩阵(按活化自由度存贮按活化自由度存贮)全部存贮量全部存贮量:NK=JP(NEQ)第四章第四章 单刚与载荷的组集，指定约束的单刚与载荷的组集，指定约束的处理，单元分析的预处理处理，单元分析的预处理子程序 PFORM 的功能子程序 MODIFY 完成将指定位移转化成等价内力将单刚组集成总刚子程序 ADDSTF子程序 BASBLY 组集载荷与节点反力向量PRTREAC 输出节点反力PRTDIS 输出节点位移ELEMLIB 单元库管理程序4-1 子程序 PFORM 的功能 子程序子程序 PFORM 起到承上启下的作用，为调用与单起到承上启下的作用，为调用与单元运算有关的各种子程序做必要的前后处理。元运算有关的各种子程序做必要的前后处理。功能参数功能参数 ISW、ASEM：ISW=3，ASEM=.TRUE.构造单刚，并组集构造单刚，并组集ISW=3，ASEM=.FALSE.构造单元载荷，组集，将指构造单元载荷，组集，将指 定位移转化为成等效载荷定位移转化为成等效载荷ISW=4 计算单元内力或应力，并输出计算单元内力或应力，并输出ISW=6 构造单元内力，并组集成结构反力构造单元内力，并组集成结构反力ISW=1 读入并构造有关的材料常数读入并构造有关的材料常数 SUBROUTINE PFORM(F,LD,UL,XL,S,P,U,X,IX,D,ID,JP,NST,NEN1,ISW)LOGICAL ASEM COMMON/ASEM/ASEM COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL DIMENSION UL(NDF,*),U(*),F(NDF,*),S(NST,*),P(*),JP(*),+XL(NDM,*),X(NDM,*),D(20,*),IX(NEN1,*),ID(NDF,*),+LD(NDF,*)IEL=0 MCT=0 ！设置打印输出 DO 110 N=1,NE CALL PZERO(UL,NST*(NST+3)MA=IX(NEN1,N)！材料号 DO 108 I=1,NEN II=IX(I,N)IF(II.NE.0)GO TO 105 DO 102 J=1,NDF102 LD(J,I)=0 GO TO 108105 IID=II*NDF-NDF NEL=I ！确定单元的实际节点数确定单元的实际节点数 DO 106 J1=1,NDM106 XL(J1,I)=X(J1,II)！确定单元的节点坐标确定单元的节点坐标DO 107 J=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+J107 LD(J,I)=K ！确定自由度标志数组确定自由度标志数组108 CONTINUE IE20=D(20,MA)IF(IE20.NE.IEL)MCT=0 ！设置打印输出 IEL=IE20 ！材料号CALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)！确定联系信息 IX IF(ASEM.AND.ISW.EQ.3)CALL ADDSTF(LD,JP,S,NST)130 IF(ISW.EQ.6)CALL BASBLY(F,P,LD,NST)120 CONTINUE IF(.NOT.ASEM.AND.ISW.EQ.3)CALL MODIFY(U,LD,S,UL(1,NEN+1),NST)IF(.NOT.ASEM.AND.ISW.EQ.3)CALL BASBLY(U,P,LD,NST)109 CONTINUE110 CONTINUE ENDUL：单元位移单元位移 S：单刚单刚 P：单元内力单元内力4-2 子程序 MODIFY 完成将指定位移转化成等价内力按活化自由度和被约束自由度，将总刚方程分块：u1：活化自由度的位移 u2：被约束自由度的位移（包括指定位移）F1 已知 F2：未知 SUBROUTINE MODIFY(U,LD,S,DUL,NS)REAL U,DUL,S DIMENSION U(1),LD(1),S(NS,1),DUL(1)DO 110 I=1,NS K=LD(I)IF(K.LE.0)GO TO 110 DO 100 J=1,NS100 U(K)=U(K)-S(I,J)*DUL(J)110 CONTINUE RETURN END4-3 将单刚组集成总刚子程序将单刚组集成总刚子程序 ADDSTFCALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)PFOR 35构造了单刚，存放于构造了单刚，存放于 S 中，中，ADDSTF 完成组集完成组集PFORM 中，LD(NDF,NEN)ADDSTF 中,LD(NST)一维数组与二维满存数组的元素对应关系：IJ=JP(I)I+J I J SUBROUTINE ADDSTF(LD,JP,S,NST)CHARACTER Y DIMENSION LD(*),JP(*),S(NST,NST)IF(Y.EQ.Y.OR.Y.EQ.y)WRITE(*,2000)N,S2000 FORMAT(ELEMENT NUMBER ,I4,ELEMENT STIFFNESS/(6E12.5)DO 200 J=1,NST K=LD(J)IF(K.LE.0)GO TO 200 L=JP(K)-K DO 100 I=1,NST M=LD(I)IF(M.GT.K.OR.M.LE.0)GO TO 100 M=L+M A(M)=A(M)+S(J,I)100 CONTINUE200 CONTINUE END4-4子程序 BASBLY 组集载荷与节点反力向量单元自由度指示数组LDISW=3 组集载荷向量，按节点活化自由度次序排列ISW=6 节点反力向量，按节点排列次序130 IF(ISW.EQ.6)CALL BASBLY(F,P,LD,NST)PFOR 37IF(.NOT.ASEM.AND.ISW.EQ.3)CALL BASBLY(U,P,LD,NST)PFOR 40 SUBROUTINE BASBLY(B,P,LD,NS)BASB 1 DIMENSION B(*),P(*),LD(*)BASB 2 DO 100 I=1,NS BASB 4 II=LD(I)BASB 5 IF(II.GT.0)B(II)=B(II)+P(I)BASB 6100 CONTINUE BASB 7 RETURN BASB 8 END BASB 94-5 PRTREAC 输出节点反力节点反力两种计算方式：节点反力两种计算方式：1）由总刚方程求得）由总刚方程求得 2）分别计算每个单元内力（单刚乘以单元节点位移），）分别计算每个单元内力（单刚乘以单元节点位移），然后组集成总体内力向量然后组集成总体内力向量因总刚分解后，已不存在了，故选因总刚分解后，已不存在了，故选 2）SUBROUTINE PRTREAC(R)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ REAL R(NDF,*),RSUM(6)NNEQ=NDF*NJ DO 50 K=1,NDF RSUM(K)=0.0 50 CONTINUE DO 100 N=1,NJ,50 J=MIN0(NJ,N+49)WRITE(6,2000)(K,K=1,NDF)DO 100 I=N,J DO 75 K=1,NDF R(K,I)=-R(K,I)RSUM(K)=RSUM(K)+R(K,I)75 CONTINUE DO 76 K=1,NDF76 IF(ABS(R(K,I).GT.1.0E-3)GO TO 77 GO TO 10077 WRITE(6,2001)I,(R(K,I),K=1,NDF)100 CONTINUE WRITE(6,2002)(RSUM(K),K=1,NDF)RETURN2000 FORMAT(/5X,NODAL REACTIONS/2X,NODE,6(I8,DOF)2001 FORMAT(I6,6F12.4)2002 FORMAT(3X,SUM,6F12.4)END 4-6 PRTDIS 输出节点位移节点位移按节点次序重新排列注意：1）解得得位移向量为活化自由度，编号小于NEQ2）约束(包括指定位移)得值 存放在二维数组 F（载荷）中 SUBROUTINE PRTDIS(ID,U,F)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ DIMENSION ID(NDF,*),U(*),F(NDF,*),UL(6)DO 102 II=1,NJ,50 WRITE(6,2000)(K,K=1,NDF)JJ=MIN0(NJ,II+49)DO 102 N=II,JJ DO 100 I=1,NDF K=ID(I,N)IF(K.LT.0)U(NEQ-K)=F(I,N)IF(K.LT.0)UL(I)=U(NEQ-K)100 IF(K.GT.0)UL(I)=U(K)WRITE(6,2001)N,(UL(I),I=1,NDF)102 CONTINUE2000 FORMAT(/1X,NODAL DISPLACEMENTS/2X,NODE,6(I8,DISP)2001 FORMAT(I6,6E12.4)END4-7 ELEMLIB 单元库管理程序FEP2 中已有两种单元：BEAM、PLANE若要增加单元，可修改第七句。例如，已编好一个板单元，名为：PLATE（形参），则程序修改如下：GOTO(1,2,3)IEL 3CALL PLATE（D,U,X,IX,S,P,I,J,K,ISW)RETURN SUBROUTINE ELEMLIB(D,U,X,IX,S,P,I,J,K,ISW)REAL P(K),S(K,K),U(1)DIMENSION IX(1),D(1),X(1)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL IF(IEL.LE.0.OR.IEL.GT.2)GO TO 400 GO TO(1,2)IEL1 CALL BEAM(D,U,X,IX,S,P,I,J,K,ISW)GO TO 1002 CALL PLANE(D,U,X,IX,S,P,I,J,K,ISW)GO TO 100100 RETURN400 WRITE(*,4000)IEL STOP4000 FORMAT(5X,*FATAL ERROR*ELEMENT CLASS NUMBER,I5,INPUT)END第五章第五章 平面单元设计平面单元设计平面四节点单元的基本公式数值积分形函数及其导数四节点单刚的构造3 9 节点平面单元平面单元 PLANE 子程序5-1平面四节点单元的基本公式 在在 FEP2 中，提供了中，提供了3 9节点平面单元，其