《平面应力应变有限元分析程序设计课程设计.doc》由会员分享，可在线阅读，更多相关《平面应力应变有限元分析程序设计课程设计.doc（20页珍藏版）》请在三一文库上搜索。
1、平面应力/应变有限元分析程序设计 引言有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。摘要有限元法的基本思想是将物体（即连续的求解域）离散成有限个且按一定方式相互联结在一起的单元的组合，来模拟或逼近原来的物体，从而将一个连续的无限自由度问题简化为离散的有限自由度问题求解的一种数值分析方法。物体被离散后，通过对其中各个单元进行单元分析，最终得到对整个物体的分析。网格划分中每一个小的块体称为单元。确定单元形状、单元之间相互联结的点称为节点。单元上节点处的
2、结构内力为节点力，外力（有集中力、分布力等）为节点载荷。三角形常应变单元求解结构主程序l 功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。l 基本思想：单元结点按右手法则顺序编号。l 荷载类型：可计算结点荷载。l 说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。1.1 程序流程图计算单刚并将单刚加入总纲引入边界条件生成载荷向量计算节点位移及支反力计算单元应力及主应力输出结果开始导入数据根据题意将得到的程序运行所需的初始数据的input.txt导入。由MATLAB中的LinearTriangleElementStiffness
3、计算单刚矩阵和LinearTriangleAssemble 依次加入总刚矩阵。根据问题将约束信息引入 计算得到载荷向量，由F=KU计算支反力，再次调用MATLAB中的LinearTriangleElementStresses 和 LinearTriangleElementPStresses计算得到到单元应力和主应力。详细说明见附录。1.2 程序应用举例【例1】如图所示，受均匀分布载荷作用的薄平板结构。将平板离散化成两个线性三角单元，如图所示。假定E=210Gpa，t=0.025m，w=3000kn/m。解答如图所示 单元连通性如下表所示单元编号节点i节点j节点m11342123由题意可得输入数
4、据文件input.txt为：4 2 2 4 210e6 0.3 0.0251 3 41 2 30.0 0.00.5 0.00.5 0.250.0 0.252 9.375 03 9.375 01 1 12 0 03 0 04 1 1说明：第一行：读入程序控制信息NPION=fscanf(FP1,%d,1) %结点个数（结点编码总数）NELEM=fscanf(FP1,%d,1) %单元个数（单元编码总数）NFORCE=fscanf(FP1,%d,1) 结点荷载个数NVFIX=fscanf(FP1,%d,1) 受约束边界点数YOUNG=fscanf(FP1,%e,1) 弹性模量POISS=fscan
5、f(FP1,%f,1) 泊松比 THICK=fscanf(FP1,%d,1) %厚度第二、三行：读入单元连接信息：LNODS=fscanf(FP1,%d,3,NELEM); %单元定义数组，单元结点号，逆时针输入第四、五、六、七行：读入结点坐标COORD=fscanf(FP1,%f,2，NPOIN); %结点坐标值，第1列为x坐标值，第2列为y坐标值第八、九行：读入结点荷载信息FORCE=fscanf(FP1,%f,3,NFORCE); %结点号，X向结点荷载数值，Y向结点荷载数值（与坐标轴方向一致为正）第十、十一、十二、十三行：读入零位移信息FIXED=fscanf(FP1,%d,3,inf
6、); %结点号，X向约束，Y向约束下述两个例子思路相同 不再赘述%-运行该程序得到输出文件output.txt如下所示节点号 x位移 y位移 1 0.00000000e+000 0.00000000e+000 2 7.11111747e-006 1.11517786e-006 3 6.53122498e-006 4.46071143e-008 4 0.00000000e+000 0.00000000e+000节点号 Fx Fy 1 -9.37500000e+000 -5.62950360e+000 2 9.37500000e+000 6.66133815e-016 3 9.37500000e+
7、000 2.22044605e-016 4 -9.37500000e+000 5.62950360e+000单元号 单元应力x 单元应力y 单元应力xy 1 3.01441153e+003 9.04323459e+002 7.20576461e+000 2 2.98558847e+003 -3.60288231e+000 -7.20576461e+000单元号 单元主应力 单元主应力 单元主应力角 1 3.01443614e+003 9.04298852e+002 1.95656990e-001 2 2.98560584e+003 -3.62025247e+000 -1.38116518e-0
8、01【例2】考虑如图所示的薄平板结构。平板中间有一个洞，并且该平板的角上受到集中载荷作用 。假定E=70GPa，v=0.25，t=0.02，P=20KN。用如图所示的有限元网格结构求解 薄平板结构 单元离散化解答由题意知单元连通性如下所示单元号节点i节点j节点m11262156323742765348638777812871211911121610111615111011151210151413910141491413155610165109由以上可得输入数据文件input.txt为：16 16 1 8 70e6 0.25 0.021 2 61 6 52 3 72 7 63 4 83 8 77
9、 8 127 12 1111 12 1611 16 1510 11 1510 15 149 10 149 14 135 6 105 10 90.0 0.00.3 0.00.6 0.00.9 0.00.0 0.30.3 0.30.6 0.30.9 0.30.0 0.60.3 0.60.6 0.60.9 0.60.0 0.90.3 0.90.6 0.90.9 0.916 0 -201 1 15 1 19 1 113 1 1运行该程序得到输出文件output.txt如下节点号 x位移 y位移 1 0.00000000e+000 0.00000000e+000 2 -2.00086423e-005 -
10、2.25205457e-005 3 -2.90761491e-005 -5.80543773e-005 4 -3.04536942e-005 -8.53667727e-005 5 0.00000000e+000 0.00000000e+000 6 -8.04661454e-007 -1.72917194e-005 7 -7.16598120e-006 -5.85435184e-005 8 -7.73311578e-006 -8.67443178e-005 9 0.00000000e+000 0.00000000e+000 10 1.15903318e-007 -1.75810380e-005
11、11 1.01419534e-006 -6.38955125e-005 12 6.37605900e-006 -9.59678117e-005 13 0.00000000e+000 0.00000000e+000 14 2.07317472e-005 -1.99029003e-005 15 3.45609853e-005 -6.35081494e-005 16 3.55507483e-005 -1.16719705e-004节点号 Fx Fy 1 1.88054201e+001 1.07884164e+000 2 -3.55271368e-015 -9.32587341e-015 3 -1.4
12、2108547e-014 -1.95399252e-014 4 3.55271368e-015 0.00000000e+000 5 1.33664311e+000 9.25376599e+000 6 4.44089210e-016 1.24344979e-014 7 0.00000000e+000 -1.24344979e-014 8 -7.10542736e-015 2.84217094e-014 9 9.10453465e-001 2.24654136e-001 10 -3.55271368e-015 -1.06581410e-014 11 0.00000000e+000 1.776356
13、84e-014 12 1.77635684e-014 -2.84217094e-014 13 -2.10525167e+001 9.44273824e+000 14 0.00000000e+000 -7.10542736e-015 15 -7.10542736e-015 -1.42108547e-014 16 -7.10542736e-015 -2.00000000e+001单元号 单元应力x 单元应力y 单元应力xy 1 -4.65457955e+003 5.64145758e+001 -3.09546055e+002 2 -2.00271295e+002 -5.00678238e+001
14、-1.61389381e+003 3 -2.28723714e+003 -6.85942203e+002 -1.27154194e+003 4 -1.25791261e+003 9.05581311e+002 -2.05779636e+003 5 -4.28569588e+002 -4.28569588e+002 -4.28569588e+002 6 -1.71588939e+002 -1.57030152e+002 -5.87125613e+002 7 -7.15059778e+002 -2.33091351e+003 -1.31521831e+003 8 1.00149532e+003 -
15、9.98424792e+002 -2.22993145e+003 9 4.32793837e+001 -4.83128855e+003 -2.70443597e+002 10 2.70443597e+002 1.57995614e+002 -1.83537811e+003 11 2.47677493e+002 1.52304088e+002 -1.19165055e+003 12 3.29747226e+003 2.82600190e+002 -2.14567781e+003 13 -1.15624386e+002 -5.70673973e+002 2.83248549e+002 14 5.1
16、5990153e+003 1.28997538e+003 -1.85760403e+003 15 -2.18273340e+002 -1.22076005e+002 -1.52797443e+003 16 2.88470479e+001 7.21176198e+000 -1.64089688e+003单元号 单元主应力 单元主应力 单元主应力角 1 7.66669049e+001 -4.67483188e+003 3.74329288e+000 2 1.49047072e+003 -1.74080984e+003 4.36678459e+001 3 1.60264731e+001 -2.989
17、20582e+003 2.89013749e+001 4 2.14863592e+003 -2.50096722e+003 3.11349466e+001 5 -5.68434189e-013 -8.57139176e+002 -4.50000000e+001 6 4.22861192e+002 -7.51480283e+002 4.46448314e+001 7 2.05628414e+001 -3.06653613e+003 -2.92189826e+001 8 2.44540809e+003 -2.44233756e+003 -3.29236390e+001 9 5.82378344e+
18、001 -4.84624701e+003 -3.16584940e+000 10 2.05045868e+003 -1.62201947e+003 -4.41226902e+001 11 1.39259511e+003 -9.92613529e+002 -4.38541986e+001 12 4.41230560e+003 -8.32233152e+002 -2.74551090e+001 13 2.01651093e+001 -7.06463468e+002 2.56130823e+001 14 5.90724620e+003 5.42630711e+002 -2.19157345e+001
19、 15 1.35855662e+003 -1.69890596e+003 4.40984990e+00116 1.65896194e+003 -1.62290313e+003 -4.48111410e+001 %-【例3】考虑如图3.1所示的由均匀分布载荷和集中载荷作用的薄平板结构。将平板离散化成12个线性三角单元，如图4所示。假定E=210GPa，v=0.3，t=0.025m，w=100kN/m和P=12.5kN。 薄平板结构 单元离散化解答如图所示 单元连通性如下表所示单元编号123456789101112节点i113344565589节点j34546767891111节点m2322568
20、8910910由题意可得输入数据文件input.txt为：11 12 1 6 210e6 0.3 0.0251 3 21 4 33 5 23 4 54 6 5 4 7 65 6 86 7 85 8 95 9 108 11 99 11 100.0 0.50.25 0.50.125 0.3750.0 0.250.25 0.250.125 0.1250.0 0.00.25 0.00.375 0.1250.5 0.250.5 0.05 0 -12.51 1 14 1 17 1 1运行该程序得到输出文件output.txt如下所示节点号 x位移 y位移 1 0.00000000e+000 0.00000
21、000e+000 2 1.15783194e-006 -3.01027541e-006 3 4.32399366e-007 -1.65732863e-006 4 0.00000000e+000 0.00000000e+000 5 1.33346561e-008 -4.17009336e-006 6 -4.08900882e-007 -1.76168681e-006 7 0.00000000e+000 0.00000000e+000 8 -1.33987898e-006 -3.52820118e-006 9 -5.29884406e-007 -4.52575409e-006 10 8.32321
22、713e-008 -5.11440654e-006 11 -1.26998147e-006 -5.29031528e-006节点号 Fx Fy 1 -6.09390871e+000 4.35913827e+000 2 -7.21644966e-016 -4.44089210e-016 3 1.77635684e-015 0.00000000e+000 4 -3.12182583e-001 3.91807210e+000 5 -4.44089210e-016 -1.25000000e+001 6 -1.77635684e-015 -3.55271368e-015 7 6.40609129e+00
23、0 4.22278963e+000 8 2.72004641e-015 1.77635684e-015 9 -1.77635684e-015 -1.42108547e-014 10 -2.72004641e-015 1.77635684e-015 11 -1.77635684e-015 3.55271368e-015单元号 单元应力x 单元应力y 单元应力xy 1 1.15305830e+003 6.01598247e+002 -8.77878249e+002 2 7.98275752e+002 2.39482726e+002 -1.07088927e+003 3 6.03981463e+00
24、2 1.15544152e+003 -8.79161518e+002 4 2.49198910e+002 7.93325998e+002 -1.07217254e+003 5 -1.66782708e+002 -5.93279398e+002 -1.07873993e+003 6 -7.54893935e+002 -2.26468181e+002 -1.13832071e+003 7 -6.47363279e+002 -7.33398416e+002 -9.11628507e+002 8 -1.23547451e+003 -3.66587199e+002 -9.71209279e+002 9
25、6.84995649e+001 -5.18639563e+002 -4.54747351e-013 10 2.27373675e-013 -1.95713043e+002 6.84995649e+001 11 -9.09494702e-013 -1.95713043e+002 -6.84995649e+001 12 -6.84995649e+001 1.27213478e+002 -2.27373675e-013单元号 单元主应力 单元主应力 单元主应力角 1 1.79748982e+003 -4.28332686e+001 -3.62815905e+001 2 1.62561599e+003
26、 -5.87857513e+002 -3.76887090e+001 3 1.80109742e+003 -4.16744346e+001 3.62935466e+001 4 1.62741458e+003 -5.84889671e+002 3.78808636e+001 5 7.19584654e+002 -1.47964676e+003 -3.94088915e+001 6 6.77900336e+002 -1.65926245e+003 3.84662997e+001 7 2.22262043e+002 -1.60302374e+003 -4.36491766e+001 8 2.6291
27、8747e+002 -1.86498045e+003 3.29499816e+001 9 6.84995649e+001 -5.18639563e+002 -4.43763713e-014 10 2.15925857e+001 -2.17305628e+002 1.74960101e+001 11 2.15925857e+001 -2.17305628e+002 -1.74960101e+001 12 1.27213478e+002 -6.84995649e+001 6.65645570e-0141.3 程序的改进要点上述三角形三节点程序反映了有限元的基本思路，可以计算简单的平面应力或平面应变
28、问题。在熟练掌握了程序的编制与使用后，可在以下几方面对程序进行改进，以加深对矩阵位移法及MATLAB语言编程的理解：1、本程序的弹性模量仅能输入一个数值，意味着程序仅能计算由同种材料构成的结构。考虑如何改进使程序可计算由不同材料构成的组合结构。2、本程序仅能计算结点集中荷载类型，考虑如何编制体积力、表面力、跨中集中力等类型的程序。3、考虑如何编制有支座位移的程序。4、本程序最后的结果没有生成输出文件，可以编制输出文件。5、本程序计算的应力没有进行结果处理，可以编制最后结果处理的程序。6、可以在此程序基础上编制三角形六节点、四边形四节点等程序。综上所述，本章的三角形三节点常应变程序体现了如何将有
29、限元法的计算方法和过程用MATLAB程序语言表达出来，重点放在程序架构和流程的建立以及算法实现方面，主要依赖手工操作手工输入初始数据（前处理）、手工绘制计算结果（后处理），目的是使学生能够清晰、明确地掌握有限元法的基本理论和概念，熟练掌握应用MATLAB语言编制、修改和调试简单程序的能力。1.4 总结通过这两周的课程设计，对有限元结合MATLAB软件解决平面应力/应变问题的方法和过程逐渐清晰，特别是在应用MATLAB软件的过程中学到了很多东西。同时，通过此次课程设计，我对以前学过的有限元理论课程知识也有了进一步了解。一开始还没有头绪，耗了半天也没正处东西来，后来把课本拿来仔细看模仿学习，逐渐对
30、软件有所熟悉，对解题思路更加清晰，我深刻的认识到遇到什么困难自己动手自己找方法的重要性。以后要加强练习以增强解决各种问题的能力。附录程序说明1 程序准备format short e %设定输出类型clear all%清除所有已定义变量clc%清屏l 说明：format short e 设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法； clear all 清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况； clc 清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD
31、LNODS YOUNG POISS THICKglobal FORCE FIXED global ASTIF ASLOD ASDISP ESTIF ELEDISP SIGMA ELPST Force astifglobal FP1l 说明：NNODE单元结点数，NPION总结点数， NELEM单元数，NVFIX受约束边界点数，NFORCE结点力数，COORD结构结点坐标数组，LNODS单元定义数组，YOUNG弹性模量，POISS泊松比，THICK厚度 FORCE 节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED 约束信息数组(n,
32、3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0ASTIF总体刚度矩阵，ASLOD总体荷载向量，ASDISP结点位移向量SIGMA单元应力，ELPST单元主应力 ， Force支反力FP1数据文件指针3 打开文件FP1=fopen(input.txt,rt); %打开输入数据文件 存放初始数据l 说明：FP1=fopen(input.txt,rt); 打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen(output.t
33、xt,wt);打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。4 读入程序控制信息NPION=fscanf(FP1,%d,1) %结点个数（结点编码总数）NELEM=fscanf(FP1,%d,1) %单元个数（单元编码总数）NFORCE=fscanf(FP1,%d,1) 结点荷载个数NVFIX=fscanf(FP1,%d,1) 受约束边界点数YOUNG=fscanf(FP1,%e,1) 弹性模量POISS=fscanf(FP1,%f,1) 泊松比 THICK=fscanf(FP1,%d,1)
34、 %厚度LNODS=fscanf(FP1,%d,3,NELEM) %单元定义数组（单元结点号）l 说明：建立LNODS矩阵，该矩阵指出了每一单元的连接信息。矩阵的每一行针对每一单元，共计 NELEM；每一列相应为单元结点号（编码）、按逆时针顺序输入。命令中，3,NELEM表示读取NELEM行3列数据赋值给LNODS矩阵。显然，LNODS(i,1:3)依次表示i单元的i，j，k结点号。COORD=fscanf(FP1,%f,2,NPION) %结点坐标数组l 说明：建立COORD矩阵，该矩阵用来存储各结点x，y方向的坐标值。从FP1文件中读取全部结点个数NPOIN的坐标值，从1开始按顺序读取。C
35、OORD(i，1：2)表示第i个结点的x，y坐标。FORCE=fscanf(FP1,%f,3,NFORCE) %结点力数组l 说明： (n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向FIXED=fscanf(FP1,%d,3,inf) %约束信息数组l 说明： (n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0l 总体说明：从输入文件FP1中读入结点个数，单元个数，结点荷载个数，受约束边界点数，弹性模量，泊松比，厚度，单元定义数组，结点坐标数组，结点力数组，约束信息数组；程
36、序中弹性模量仅输入了一个值，表明本程序仅能求解一种材料构成的结构，如：钢筋混凝土结构、钢结构，不能求解钢筋混凝土钢组合结构。采用了命令fscanf，其中：%d表示读入整数格式，%f表示读入浮点数；1表示读取1个数，A,B形式表示读A行B列数组，A,B表示将A,B转置，inf表示正无穷。5 计算单刚并生成总刚及加入约束信息 %单刚加入总刚-ASTIF(1:2*NPION,1:2*NPION)=0; %张成特定大小总刚矩阵并置0l 说明：建立单元刚度矩阵ASTIF，该矩阵的行列数均为2*NPION ，NPION表示结点数，每个结点有两个方向的力和位移。%-for i=1:NELEM %第i个单刚矩
37、阵 ESTIF=LinearTriangleElementStiffness(YOUNG,POISS,THICK,COORD(LNODS(i,1),1), COORD(LNODS(i,1),2),COORD(LNODS(i,2),1), COORD(LNODS(i,2),2),COORD(LNODS(i,3),1), COORD(LNODS(i,3),2),1); % 第i个单刚矩阵加入总刚矩阵 ASTIF=LinearTriangleAssemble(ASTIF,ESTIF,LNODS(i,1),LNODS(i,2),LNODS(i,3); EndAstif= ASTIF%引入边界条件-%将
38、约束信息加入总刚（置0置1法）NUM=1; %计数器(当前已分析的节点数) i=0; %计数器(当前已处理的约束数) tmp(NVFIX)=0; %临时存被约束的序号while iNVFIX for j=-1:0 if FIXED(NUM,j+3)=1 %若发现约束 i=i+1; %计数器自增 tmp(i)=FIXED(NUM)*2+j; %求约束序号 end end NUM=NUM+1;endl 说明： tmp(NVFIX)=0，形成一个元素值均为0的一行NVFIX列的行向量，执行while语句，首先判断i是否小于控制数据NVFIX，若小于则往下进行，若不小于则退出。执行for语句，FIXED(NUM,j+3)表示约束信息数组中第NUM行第j+3列的元素，j从-1到0，即j+3表示2到3列，即约束信息数组中描述结点x和y方向受约束的情况，判断FIXED(NUM,j+3)若等于1，则约束数自增，若不等于1，跳出。FIXED(NUM)表示FIXED(NUM，1)，tmp(i)=FIXED(NUM)*2+j计算整体约束序号，将序号放入tmp行向量中的i列。%-