《[计算机软件及应用]材料成形有限元法_有限元程序设计.doc》由会员分享，可在线阅读，更多相关《[计算机软件及应用]材料成形有限元法_有限元程序设计.doc（30页珍藏版）》请在三一文库上搜索。
1、华 中 科 技 大 学 教 案章节第一章 第一节周次第三周课次第一次教案内容一、 有限元程序设计的基本思想：1. 一般采用FORTRAN77语言，FORTRAN90语言也开始逐渐被采用。2. 程序流程标准化（通用性、可读性）。3. 程序模块化（可移植性、可读性、可调性）。4. 程序变量名统一化（可读性、可调性）。二、 有限元程序的一般过程： 最小势能原理（位移法的数学基础）： （1）其中：e为构形的应变；s为构形的应力；u为构形的虚位移；P为构形上所受的面力；G为构形的体力；PP为构形的势能。对于线弹性问题而言，构形就是线弹性体。某构形处于平衡状态时，其势能应为最小，即 （2） （3）1. 单
2、元离散过程：1) 几何关系（单元内任意点的位移u与应变e之间的关系） （4）其中L为单元几何关系算子。2) 本构关系（单元内任意点的应变e与应力s之间的关系） （5）其中为弹性矩阵。3) 插值关系（单元内任意点的位移u与单元结点位移ue之间的关系） （6）其中N为单元插值形函数矩阵对于某一有限元单元来说，由式（3）式（6）可得：华 中 科 技 大 学 教 案章节第一章 第一节周次第三周课次第一次教案内容 （7）其中B=LN。由于的任意性，得 （8）最终得到单元刚度方程 （9）其中：单元刚度矩阵 单元外载荷向量2. 单元集成过程 （10） （11）3. 刚度方程约束过程给出所要求解问题的边界条件
3、。4. 刚度方程求解过程由式（11）解出所有结点的位移U。5. 有限元回代过程1) 求单元的应变 （12）2) 求单元的应力 （13）华 中 科 技 大 学 教 案章节第一章 第二节周次第四周课次第二次教案内容一、 单元几何关系1. 三维情况，单元的应变与位移的关系为： （1）即 （2）2. 二维情况单元的应变与位移的关系可由三维情况退化而来，1) 平面应力和平面应变情况： （3）华 中 科 技 大 学 教 案章节第一章 第二节周次第四周课次第二次教案内容2) 轴对称情况 （4）二、 单元本构关系1. 三维情况单元的应力应变关系为 （5）其中 （6）2. 二维情况可由三维情况退化而来，1) 平
4、面应力（） （7）2) 平面应变（） （8）华 中 科 技 大 学 教 案章节第一章 第二节周次第四周课次第二次教案内容3) 轴对称情况（） （9）华 中 科 技 大 学 教 案章节第一章 第三节周次第五周课次第三次教案内容一、 有限元法中常采用的单元类型及其形状函数1. 线性单元（一次单元）：1) 一维线单元 （1）其中，l为单元长度；x1和x2为单元结点坐标；x为单元内任意点的坐标。2) 二维面单元a) 三角形单元 （2）其中，A为单元面积 （3） （4）xi、 yi、zi（i =1,2,3）为单元结点坐标。b) 四边形单元 （5）四边形四个结点的自然坐标为（r1，s1）=（-1，-1）；
5、（r2，s2）=（1，-1）（r3，s3）=（1，1）；（r4，s4）=（-1，1）华 中 科 技 大 学 教 案章节第一章 第三节周次第五周课次第三次教案内容3) 三维体单元 （5）其中，V为单元体积 （6） （7） （8） （9） （10）循环轮换脚标1、2、3、4，可得其它常数。华 中 科 技 大 学 教 案章节第一章 第三节周次第五周课次第三次教案内容2. 高阶单元1) 一维单元a) 一维二次单元： （11）b) 一维三次单元： （12）2) 二维单元a) 二维三角形单元（六结点三角形单元）： （13）b) 二维四边形单元（八结点四边形单元）： （14）华 中 科 技 大 学 教 案章
6、节第二章 第一节周次第六、七周课次第四、五次教案内容一、平面三结点三角形单元的有限元一般过程1. 有限元平衡方程1) 几何关系 （1）2) 本构关系a) 平面应力（） （2）b) 平面应变（） （3）c) 轴对称情况（） （4）3) 插值关系 （5）华 中 科 技 大 学 教 案章节第二章 第一节周次第六、七周课次第四、五次教案内容其中，ue假设为 则形函数矩阵可表示为 （6）将公式（1）（6）代入最小势能原理可得有限元平衡方程 （7）其中：单元刚度矩阵 单元外载荷向量1) 平面应力和平面应变情况的B矩阵为 （8）2) 轴对称情况的 B矩阵为 （9）华 中 科 技 大 学 教 案章节第二章 第
7、一节周次第六、七周课次第四、五次教案内容1) 对于平面应力情况，假设板厚为t，单元面积为a，则 （9）2) 对于平面应变情况，任取厚向长度为t，单元面积为a，则 （10）3) 对于轴对称情况，单元面积为a，则 （11）2. 有限元平衡方程组的边界条件约束处理由所有单元刚度方程组装所形成的总体刚度方程在未进行边界条件处理之前是奇异的，方程组有无数组解。只有在正确实施所要分析问题的边界条件之后，总体刚度方程才有唯一一组正确的解。有限元平衡方程组的边界条件约束处理方法一般有两种，1) 乘大数法：在所要约束的自由度位置的总刚对角线元素上乘以一个相对其是很大的数。2) 赋0赋1法：通过边界条件（即约束方
8、程）的引入，并根据初等变换使平衡方程仍然保持原有的对称性。具体实施如下：假设边界条件（即约束方程）为，原方程组为 （12）则将约束方程代入以后得， （13）华 中 科 技 大 学 教 案章节第二章 第一节周次第六、七周课次第四、五次教案内容 上式写成矩阵形式 （14）将方程第3行乘以a加到第2行，得 （15）将方阵的第三行和第三列按如下方式赋0赋1，就得到了约束后仍然是对称的方程组 （16）从式（16）可以解出x1，x2， x3，x4的值，但其中x3不是方程组的真实解，它必须通过约束方程来求得。两种边界条件约束处理方法的比较：1) 乘大数法是一种近似方法，其近似程度取决于所乘大数的大小；而且这
9、种方法只能处理约束位移是零的情况，即xi=0。2) 赋0赋1法是一种完全精确的约束处理方法；而且这种方法不仅能处理约束位移是零的情况，也可以处理约束位移是非零的情况，即类似于的情况，甚至 更一般的约束情况。3) 这两种约束处理方法都可以在组装后的总刚上直接进行约束处理，也可以在单元刚度方程上先进行约束处理，然后再进行单元刚度方程的组装。对于三角形单元而言，单刚组装总刚、方程组求解以及应力应变的回代过程与标准有限元法完全相同。华 中 科 技 大 学 教 案章节第二章 第二节周次第八周课次第六次教案内容一、平面四结点四边形单元的有限元一般过程1. 有限元平衡方程四结点四边形单元的几何关系、本构关系
10、与平面三结点三角形单元的完全相同，这里不再给出。而插值关系为 （1）其中， （2） （3）单元平衡方程为 （4）其中：单元刚度矩阵 单元外载荷向量1. 平面应力/应变问题B矩阵为 （5）华 中 科 技 大 学 教 案章节第二章 第二节周次第八周课次第六次教案内容2. 轴对称问题的B矩阵 （6）形函数Ni为 （7）B矩阵的具体求法如下：把Ni看成是x、y的函数对r和s求导， （8）其中为雅可比矩阵，表示xy坐标与自然坐标系rs之间的变换矩阵，这里等于单元的面积。由于此单元是等参单元，所以 （9）华 中 科 技 大 学 教 案章节第二章 第二节周次第八周课次第六次教案内容即 （10） （11）所以
11、， （12）由于 是关于自然坐标r、s的函数，所以单元刚度矩阵（平面应力/应变情况）也是关于 r、s的函数。一般要采用数值积分的方法（高斯积分方法）。令，并将积分限取为-1到+1，即 （13）高斯积分：一维情况， （14）二维情况， （15）一般来说，如果被积函数是一个小于或等于（2n-1）阶的多项式，n点高斯积分是精确的。对于四结点四边形单元来说，常规面内取4个高斯积分点（n=2） 是精确积分模式，面内只取1个积分点（n=1）是减缩积分模式。而实际上，多项式矩阵的阶次是3，理论上积分点数n=1.5，使此单元的数值积分有一定的困难。华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一
12、周课次第七、八、九次教案内容三结点三角形单元的有限元程序讲解：*=* ANALYSIS OF PLATES UNDER IPLANE LOADS *=* COMMON ALL(10000) !定义无名共用区 CALL OPFILE() !打开文件子程序 CALL INITDAT(NN,NE,ND,NFIX,E,ANU,T) !输入主控信息子程序 CALL ALLOC(NN,NE,ND,NFIX,MLOC,MCX,MCY,MIFIX,MP,MGS,MSTRES, & MSTRAIN) !地址(无名共用区)分配子程序 CALL CST(NN,NE,ND,ALL(MLOC),ALL(MCX),ALL
13、(MCY),E,ANU,T,NFIX, & ALL(MIFIX),ALL(MP),ALL(MGS),ALL(MSTRES),ALL(MSTRAIN)!子程序 CALL CLFILE() !关闭文件子程序 STOP !停止运行 END*=* SUBROUTINE OPEN FILE *=* SUBROUTINE OPFILE() OPEN(10,FILE=INP.DAT) OPEN(30,FILE=OUTP.DAT) RETURN END*=* SUBROUTINE CLOSE FILE *=* SUBROUTINE CLFILE() CLOSE(10) CLOSE(30) RETURN END
14、华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容*=* SUBROUTINE READ INITIAL DATA *=* SUBROUTINE INITDAT(NN,NE,ND,NFIX,E,ANU,T)* NN:结点数 NE:单元数 ND:自由度数 NFIX:约束数* E :弹性模量 ANU:泊松比 T:板厚 READ(10,*) NN,NE,ND,NFIX READ(10,*) E,ANU,T RETURN END*=* SUBROUTINE ALLOCATE MEMORY *=* SUBROUTINE ALLOC(NN,NE,ND,NFIX,ML
15、OC,MCX,MCY,MIFIX,MP,MGS,MSTRES,& MSTRAIN) MLOC =1 MCX =MLOC +NE*3 MCY =MCX +NN MIFIX =MCY +NN MP =MIFIX+NFIX MGS =MP +ND MSTRES =MGS +ND*ND MSTRAIN=MSTRES+NE*3 RETURN END*=* SUBROUTINE CST *=* SUBROUTINE CST(NN,NE,ND,LOC,CX,CY,E,ANU,T,NFIX,IFIX,P,GS,STRES,& STRAIN) DIMENSION LOC(NE,3),CX(NN),CY(NN),
16、IFIX(NFIX),P(ND),GS(ND,ND), & STRES(NE,3),STRAIN(NE,3),AA(6,6),N(6),CXL(3),CYL(3), & QG(6),QL(6),DD(3,3),BB(3,6),DB(3,6),SLOC(3),AT(18) REAL KL(6,6),KG(6,6),LAMDA(6,6),LAMDAT(6,6)华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容 CALL READDAT(LOC,CX,CY,IFIX,NE,NN,NFIX) !数据输入子程序 DO 95 I =1,ND DO 95 J =1,N
17、D95 GS(I,J)=0.0 CALL RDFORCE(P,ND) !右端项结点载荷输入子程序 DO 100 I=1,NE J =LOC(I,1) K =LOC(I,2) L =LOC(I,3)* 形成转换矩阵 LAMDA 子程序 CALL LAMBDA(J,K,L,NN,CX,CY,C1,C2,CXL,CYL,E,T,ANU,LAMDA,AREA) CALL CALCKL(KL,C1,C2,CXL,CYL,ANU) !形成局部坐标系的单元刚度 KL(半) CALL SYMTRA(KL,LAMDA,LAMDAT) !形成转置 LAMDAT 及所有 KL CALL MATMUL(KL,LAMD
18、A,AA,6,6,6) CALL MATMUL(LAMDAT,AA,KG,6,6,6) CALL ASSEMB(J,K,L,N,KG,GS,ND) !单刚组装总刚子程序100 CONTINUE CALL PROCEGS(GS,IFIX,P,NFIX,ND) !约束处理子程序 CALL DECOMP(ND,GS,AT) !方程求解子程序 CALL SOLVE(ND,GS,P) DO 200 I=1,NE J=LOC(I,1) K=LOC(I,2) L=LOC(I,3) CALL LAMBDA(J,K,L,NN,CX,CY,C1,C2,CXL,CYL,E,T,ANU,LAMDA,AREA) CAL
19、L TRANSF(J,K,L,ND,N,QG,QL,P,LAMDA) !形成局部坐标系的位移向量 QL CALL MATRIXB(CXL,CYL,AREA,BB) !形成 B 矩阵(BB) CALL MATRIXD(DD,ANU,E) !形成弹性矩阵 DD CALL MATMUL(DD,BB,DB,3,3,6)* 求单元应力、应变子程序子程序 CALL BKBRING(I,NE,SLOC,DB,BB,QL,LAMDA,STRES,STRAIN)200 CONTINUE CALL OUTPUT(P,STRES,STRAIN,NE,ND) !结果输出子程序 RETURN !主子程序返回 END华
20、中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容*=* SUBROUTINE LAMBDA *=* SUBROUTINE LAMBDA(J,K,L,NN,CX,CY,C1,C2,CXL,CYL,E,T,ANU, & LAMDA,AREA) REAL LAMDA(6,6) DIMENSION CX(NN),CY(NN),CXL(3),CYL(3) D12=SQRT(CX(K)-CX(J)*2+(CY(K)-CY(J)*2) DCL12=(CX(K)-CX(J)/D12 DCM12=(CY(K)-CY(J)/D12 D14=DCL12*(CX(L)-CX(J)
21、+DCM12*(CY(L)-CY(J) D43=SQRT(CX(L)-CX(J)*2+(CY(L)-CY(J)*2-D14*2) CXL(1)=0.0 CYL(1)=0.0 CXL(2)=0.0 CYL(2)=D12 CXL(3)=D43 CYL(3)=D14 DCL43=(CX(L)-CX(J)-DCL12*D14)/D43 DCM43=(CY(L)-CY(J)-DCM12*D14)/D43 AREA=(CXL(3)-CXL(2)*(CYL(2)-CYL(1)- & (CXL(2)-CXL(1)*(CYL(3)-CYL(2)/2.0 C1=E*T/(4.0*AREA*(1.0-ANU*2)
22、C2=E*T/(8.0*AREA*(1.0+ANU) DO 10 II=1,6 DO 10 JJ=1,6 10 LAMDA(II,JJ)=0.0 LAMDA(1,1)=DCL43 LAMDA(1,2)=DCM43 LAMDA(2,1)=DCL12 LAMDA(2,2)=DCM12 LAMDA(3,3)=DCL43 LAMDA(3,4)=DCM43 LAMDA(4,3)=DCL12华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容 LAMDA(4,4)=DCM12 LAMDA(5,5)=DCL43 LAMDA(5,6)=DCM43 LAMDA(6,5)=D
23、CL12 LAMDA(6,6)=DCM12 RETURN END*=* SUBROUTINE READ DATA *=* SUBROUTINE READDAT(LOC,CX,CY,IFIX,NE,NN,NFIX)DIMENSION LOC(NE,3),CX(NN),CY(NN),IFIX(NFIX)DO I=1, NEDO J=1,3READ(10,*) LOC(I,J)CONTINUE DO I=1,9READ(10,*) CX(I),CY(I)CONTINUEDO J=1,NFIXREAD(10,*) IFIX(J)CONTINUERETURN END*=* SUBROUTINE READ
24、 EXTERNAL FORCE *=* SUBROUTINE RDFORCE(P,ND) DIMENSION P(ND) READ(10,*)(P(I),I=1,ND) RETURN END华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容*=* SUBROUTINE CALCULATE KL *=* SUBROUTINE CALCKL(KL,C1,C2,CXL,CYL,ANU) DIMENSION CXL(3),CYL(3) REAL KL(6,6) Y32 =CYL(3) -CYL(2) X32 =CXL(3) -CXL(2) Y31 =CYL(3)
25、 -CYL(1) X31 =CXL(3) -CXL(1) Y21 =CYL(2) -CYL(1) X21 =CXL(2) -CXL(1) KL(1,1) =C1*(Y32*2) +C2*(X32*2) KL(2,1) =-C1*ANU*Y32*X32 -C2*X32*Y32 KL(2,2) =C1 *(X32*2) +C2 *(Y32*2) KL(3,1) =-C1*Y32*Y31 -C2*X32*X31 KL(3,2) =C1*ANU*X32*Y31 +C2*Y32*X31 KL(3,3) =C1*(Y31*2) +C2*(X31*2) KL(4,1) =C1*ANU*Y32*X31 +C2
26、*X32*Y31 KL(4,2) =-C1*X32*X31 -C2*Y32*Y31 KL(4,3) =-C1*ANU*Y31*X31 -C2*X31*Y31 KL(4,4) =C1*(X31*2) +C2*(Y31*2) KL(5,1) =C1*Y32*Y21 +C2*X32*X21 KL(5,2)=-C1*ANU*X32*Y21-C2*X21*Y32 KL(5,3)=-C1*Y31*Y21-C2*X31*X21 KL(5,4)=C1*ANU*X31*Y21+C2*Y31*X21 KL(5,5)=C1*(Y21*2)+C2*(X21*2) KL(6,1)=-C1*ANU*Y32*X21-C2*
27、X32*Y21 KL(6,2)=C1*X32*X21+C2*Y32*Y21 KL(6,3)=C1*ANU*Y31*X21+C2*X31*Y21 KL(6,4)=-C1*X31*X21-C2*Y31*Y21 KL(6,5)=-C1*ANU*Y21*X21-C2*X21*Y21 KL(6,6)=C1*(X21*2)+C2*(Y21*2) RETURN END华 中 科 技 大 学 教 案章节第三章 第一节周次第九 十一周课次第七、八、九次教案内容*=* SUBROUTINE SYMETRY TRANSFORM *=* SUBROUTINE SYMTRA(KL,LAMDA,LAMDAT) REAL KL(6,6),LAMDA(6,6),LAMDAT(6