有限元程序设计课件.ppt

《有限元程序设计课件.ppt》由会员分享，可在线阅读，更多相关《有限元程序设计课件.ppt（84页珍藏版）》请在三一办公上搜索。

1、教学程序：FEP1 李明瑞编 FEP2 李明瑞编 ZFEP,第一章 绪论,有限元程序的基本内容 有限元求解器,1.1 有限元程序的基本内容,有限元法的解题步骤,求应力等,有限元程序的基本内容：,数据输入阶段,前处理,一、数据输入阶段前处理,读入和生成数据，形成有限元网格，为有限元矩阵计算作准备。（数据或图形输入）,标题和控制信息 计算存贮分配 节点坐标和约束信息 单元信息 材料信息 载荷信息,二、有限元矩阵的计算、组集和求解处理器,计算插值函数矩阵 N、几何矩阵B、Jacob矩阵 J 在高斯点进行数值积分，求得单刚、单元载荷组集成总刚、总载 求解,三、数据输出阶段后处理,输出结果（场变量和场变

2、量导数等），打印结果。场变量包括：位移、温度、流场的速度势等；场变量导数包括：应力、应变、热流、流场速度势等,输出形式：数据输出和图形输出,1-2 有限元求解器,带宽法：只存储总刚矩阵的半带宽内的元素。等带宽法：每行的带宽相等，常采用二维数组存储变带宽法：每行的带宽不等，常采用一维数组存储,一、带宽法,变带宽分块求解思想：一个变量的消元只影响刚度矩阵的一部分元素，较大的节点分量方程组可以分成较小的局部系数矩阵按节点逐步求解。,波前法是另一种分块储存的变带宽法，其消元次序是按单元编号进行，边组集边消元。调入内存的单元所保留的节点称为波前节点，所消去的节点称为消元节点，消元节点是与以后调入的单元没

3、有联系的节点，即该节点的方程已组集完。波前法的回代按消元节点的反序进行。对内存的需求取决与最大波前宽。,二、波前法,V-Fortran,Fortran语言的基本格式变量基本语句子例子程序 函数子程序 其它功能、模块,第二章 FEP2 程序的总体设计与输入数据,FEP2 的设计任务 FEP2 的结构设计 FEP2 的主程序 FEP2 的主控程序 FEP2 的数据输入格式与程序实现,2-1 FEP2 的设计任务,1.FEP2的简要题目说明 FEP2是一个具有通用性的教学程序，可用于计算一般的线性静力学问题。已设计了平面梁单元与平面3-9节点元两种单元，但留下接口。2.支持软件与硬件 FORTRAN

4、77以上编译器、各种微机3.FEP2的功能,3.FEP2的功能1）输入文件名由用户自由定义，但限制为4个字符，输入文件扩展名为“DAT”，输出文件扩展名为“OUT”。2）节点坐标、单元信息等具有线性自动生成功能。3）可以处理多种工况、多种类型单元组合结构问题。4）可以处理固定约束和指定位移。5）采用变带宽存储、三角分解法求解刚度方程。6）为多种单元留下接口。,2-2 FEP2 的结构设计,FEP2：由形式主程序、主控程序、功能模块组成。主程序的主要功能：定义输入、输出文件名，调用主控程序PCONTR,主控程序PCONTR的主要功能模块：1）内存分配 2）网格生成3）变带宽存储设置 4）单刚的形

5、成与组集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)规

6、定输入文件名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,.

7、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（编译器）版本和计算机内存有关，确定了程序的计算规模。这两个

8、数组也可合并成一个。,2-4 FEP2 的主控程序,PCONTR：实质上的主程序。分以下几大段：,一、程序头（前12语句）,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

9、/NODAL FORCE/DISPL/,FO：输出文件名 ERR：逻辑变量，作出错信息控制（.TRUE.为错）ASEM：逻辑变量，作为判断是否要组集总刚 数组A：存放总刚的元素 数组M：存放单刚、单元信息、节点坐标等 NDF 最大自由度数 NDM 空间维数 NEN 单元的最多节点数 NJ 节点总数 NE 总单元数 NEQ 总方程数NUMMAT 材料类型数 FD,FORC：文字变量，作为输出的标题用 P：文字变量，判断是否要输出单刚,二、输入控制信息、确定主要数组的地址（1434语句）,NEN1=NEN+1 每个单元信息的存储量（节点号+材料号）NNEQ=NJ*NDF 总节点 节点自由度数 NS

10、T=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)N

11、J,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，NEN,2）坐标信息（一组）N，NG，（XL(I)，I=1，NDM）节点号 节点号增量 坐标

12、或载荷 由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

13、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 106 WRITE(6,3000)L,CD WRITE(*,3

14、000)L,CD ERR=.TRUE.GO TO 102 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：平面梁作用空间载荷

15、）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,

16、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,N

17、UMMAT 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

18、，0，（平面元4个0，梁元5个0）IDL(I)：约束信息，被约束的自由度给非零值（1或-1），其它为0。1：该自由度被约束，但不继续生成-1：该自由度被约束，要求继续生成,6）载荷信息 载荷输入仍由 GENVEC 完成，维数：NDF,对于指定位移，生成方式与节点载荷相同，5）输入约束信息，6）输入其值。,第三章 总刚度矩阵的变带宽存贮与求解,总刚度矩阵的变带宽存贮与对角线地址LDLT 三角分解自由度指示矩阵ID与对角线地址矩阵,3-1 总刚度矩阵的变带宽存贮与对角线地址,总刚度矩阵的性质：对称、正定、稀疏只存贮半带宽内的元素（下三角、变带宽）,稀疏对称矩阵A中的元素可用两个一维数组B和JP完全

19、确定。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

20、(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 CONTINUE,C.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

21、,JJ)IF(LL.LE.0)GO TO 100 M=MAX0(KK,LL)JP(M)=MAX0(JP(M),IABS(KK-LL)100 CONTINUE200 CONTINUE300 CONTINUECONTINUE500 CONTINUE,C.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

22、 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,S FORMAT(ELEMENT NUMBER,I4,ELEMENT STIFFNESS/(6E12.5),DO 200 J

23、=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 END,3-2 LDLT 三角分解,求解方程：Ax=b A 对称、正定矩阵=Ux=y U 单位上三角矩阵,对A进行三角分解=A=LLTL 下三角矩阵，通过比较系数法确定，需要开方运算,LDLT三角分解=A=LDLTL 下三角矩阵，D 对角矩阵,i=1,n;j=2,i,LDLT(A,JP,N)DIME

24、NSION 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）以各节点自由度为研究对象

25、，利用单元信息IX，找出与该点联系的节点，然后比较节点自由度之最大者。运算量：NJ*NDF*NE*NEN*NDF,2）逐一研究各单元，分别取出各节点的各个自由度与该单元中其它节点的自由度比较，差值最大者再与所研究的自由度之带宽值（初值为零）比较，取其中最大者代替原有的带宽值。运算量：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

26、)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.构

27、造单刚，并组集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,*)

28、,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,

29、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(

30、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：单刚

31、 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 END,4-3 将单刚组集成

32、总刚子程序 ADDSTF,CALL 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(ELEMEN

33、T 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 END,4-4子程序 BASBLY 组集载荷与节点反力向量,单元自由度指示数组LDISW=3 组集载荷向量，按节点活化自由度次序排列ISW=6 节点反力向量，按节点排列次序,130 IF(ISW.EQ.6)CALL BAS

34、BLY(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 9,4-5 PRTREAC 输出节点反力,节点反力两种计算方式：1）由总刚方程求得

35、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

36、)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)RETURN 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）解得得位移向量为活化自由度

37、，编号小于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

38、(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)END,4-7 ELEMLIB 单元库管理程序,FEP2 中已有两种单元：BEAM、PLANE若要增加单元，可修改第七句。例如，已编好一个板单元，名为：PLATE（形参），则程序修改如下：,GOTO(1,2,3)IEL CALL PLATE（D,U,X,IX,S,P,I,J,K,ISW)RETURN,SUBROUTINE ELEMLIB(D,U,

39、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

40、 STOP FORMAT(5X,*FATAL ERROR*ELEMENT CLASS NUMBER,I5,INPUT)END,第五章 平面单元设计,平面四节点单元的基本公式数值积分形函数及其导数四节点单刚的构造3 9 节点平面单元平面单元 PLANE 子程序,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,ND

41、M,NEL,IX,.FALSE.)SHAP2(SS,TT,SHP,IX,NEL)PSTRES(SIG,SIG(4),SIG(5),SIG(6),5-2 数值积分,SUBROUTINE PGAUSS(L,LINT,R,Z,W)COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL DIMENSION LR(9),LZ(9),LW(9),R(3),Z(3),W(3),G4(4),H4(4)DATA LR/-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/LINT=L*L GO TO(1,2,3,4),L

42、1 R(1)=0.0E0 Z(1)=0.0E0 W(1)=4.0E0 IF(NEL.EQ.3)THEN Z(1)=-1.0E0/3.0E0 W(1)=3.0E0 END IF RETURN,2 G=1.0E0/SQRT(3.0E0)DO 21 I=1,4 R(I)=G*LR(I)Z(I)=G*LZ(I)21 W(I)=1.0E0 RETURN3 G=SQRT(0.6E0)H=1.0E0/81.0E0 DO 31 I=1,9 R(I)=G*LR(I)Z(I)=G*LZ(I)31 W(I)=H*LW(I)RETURN,4 G=SQRT(4.8E0)H=SQRT(30.0E0)/36.0E0 G4(

43、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+H H4(3)=0.50E0+H；H4(4)=0.50E0-H I=0 DO 41 J=1,4 DO 41 K=1,4 I=I+1 R(I)=G4(K)Z(I)=G4(J)W(I)=H4(J)*H4(K)41 CONTINUE END,5-3 形函数及其导数,SUBROUTINE SHAPE(SS,TT,X,SHP,XSJ,NDM,NEL,IX,FLG)DIMENSION SHP(3,*),X(

44、NDM,1),S(4),T(4),XS(2,2),SX(2,2),IX(*)DATA S/-0.5E0,0.5E0,0.5E0,-0.5E0/,T/-0.5E0,-0.5E0,0.5E0,0.5E0/DO 100 I=1,4 SHP(3,I)=(0.5e0+S(I)*SS)*(0.5e0+T(I)*TT)SHP(1,I)=S(I)*(0.5e0+T(I)*TT)100 SHP(2,I)=T(I)*(0.5e0+S(I)*SS)IF(NEL.GE.4)GO TO 120 C.FORM TRIANGLE BY ADDING THIRD AND FOURTH TOGETHER DO 110 I=1,

45、3110 SHP(I,3)=SHP(I,3)+SHP(I,4)！3 节点平面单元 C.ADD QUADRATIC TERMS IF NECESSARY120 IF(NEL.GT.4)CALL SHAP2(SS,TT,SHP,IX,NEL),C.CONSTRUCT JACOBIAN AND ITS INVERSE DO 130 I=1,NDM DO 130 J=1,2 XS(I,J)=0.0 DO 130 K=1,NEL130 XS(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)RETURN SX(1,

46、1)=XS(2,2)/XSJ SX(2,2)=XS(1,1)/XSJ SX(1,2)=-XS(1,2)/XSJ SX(2,1)=-XS(2,1)/XSJC.FORM GLOBAL DERIVATIVES DO 140 I=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)140 SHP(1,I)=TP END,5-4 四节点单刚的构造,比较三种不同的算法，选运算量最少的一种。1）按矩阵直接运算2）预先算出D Bj，再与BiT相乘，展开3）按矩阵运算全部展开,计算：BiT D B

47、j,3 L=D(5)IF(L*L.NE.LINT)CALL PGAUSS(L,LINT,SG,TG,WG)DO 320 L=1,LINT CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.)XSJ=XSJ*WG(L)*D(7)J1=1 DO 320 J=1,NEL W11=SHP(1,J)*XSJ W12=SHP(2,J)*XSJ K1=J1 DO 310 K=J,NEL S(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)

48、+W12*SHP(1,K)S(J1+1,K1+1)=S(J1+1,K1+1)+W12*SHP(2,K)310 K1=K1+NDF320 J1=J1+NDF,5-5 39 节点平面单元,3 节点平面单元 SHAPE 1213,59 节点平面单元,SUBROUTINE SHAP2(S,T,SHP,IX,NEL)DIMENSION IX(*),SHP(3,*)S2=(1.0E0-S*S)/2.0E0 T2=(1.0E0-T*T)/2.0E0 DO 100 I=5,9 DO 100 J=1,3 100 SHP(J,I)=0.0 IF(IX(5).EQ.0)GO TO 101 SHP(1,5)=-S*(

49、1.0E0-T)SHP(2,5)=-S2 SHP(3,5)=S2*(1.0E0-T)101 IF(NEL.LT.6)GO TO 107 IF(IX(6).EQ.0)GO TO 102 SHP(1,6)=T2 SHP(2,6)=-T*(1.0E0+S)SHP(3,6)=T2*(1.0E0+S),102 IF(NEL.LT.7)GO TO 107 SHAP 19 IF(IX(7).EQ.0)GO TO 103 SHAP 20 SHP(1,7)=-S*(1.0E0+T)SHAP 21 SHP(2,7)=S2 SHAP 22 SHP(3,7)=S2*(1.0E0+T)SHAP 23103 IF(NEL

50、.LT.8)GO TO 107 SHAP 24 IF(IX(8).EQ.0)GO TO 104 SHAP 25 SHP(1,8)=-T2 SHAP 26 SHP(2,8)=-T*(1.0E0-S)SHAP 27 SHP(3,8)=T2*(1.0E0-S)SHAP 28C.INTERIOR NODE(LAGRANGIAN)104 IF(NEL.LT.9)GO TO 107 SHAP 30 IF(IX(9).EQ.0)GO TO 107 SHAP 31 SHP(1,9)=-4.0E0*S*T2 SHAP 32 SHP(2,9)=-4.0E0*T*S2 SHAP 33 SHP(3,9)=4.0E0*