平面三角形单元有限元程序设计

《平面三角形单元有限元程序设计》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、一、题目如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。试计算薄板的位移及应力分布。要求：1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形单元，单元数不得少于30个）；2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）；3.提交程序编写过程的详细报告及计算机程序；

2、4.所有同学参加答辩，并演示有限元计算程序。有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。2）单元分析，建立单元刚度矩阵。123）整体分析，建立总刚矩阵。4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。6）解方程，求出节点位移。7）求出各单元的单元应力。8）计算结果整理。一、程序设计网格划分如图，将薄板如图划分为6行，并建立坐标系，则PPXXYY单元编号节点编号12刚度矩阵的集

3、成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种；（2）该单元对应总刚度矩阵的那几行哪几列（3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环：（1）最外层：逐行计算（2）中间层：该行逐个计算（3）最里层：区分为第奇/偶数个计算单元刚度的集成：边界约束的处理：划0置1法适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约

4、束。  做法：  （1） 将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎ 这样，由该方程求得的此位移值一定等于已知量。 12（2） 将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎ 当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。 ◎ 若约束为零位移约束时，此步则可省去。 特点： （1） 经以

5、上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2） 但这种方法不改变方程阶数，利于存贮。 （3） 不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。程序如下：变量说明NNODE单元节点数NPION总结点数NELEM单元数NVFIX受约束边界点数FIXED约束信息数组NFORCE节点力数FORCE节点力数组COORD结构节点坐标数组LNODS单元定义数组YOUNG弹性模量POISS泊松比THICK厚度B单元应变矩阵(3*6)D单元弹性矩阵(3*3)S单元应力矩阵(3

6、*6)12A单元面积ESTIF单元刚度矩阵ASTIF总体刚度矩阵ASLOD总体荷载向量ASDISP节点位移向量ELEDISP单元节点位移向量STRESS单元应力%**********************************************************%初始化clearformatshorte%设定输出类型clear%清除内存变量NELEM=36%单元个数（单元编码总数）NPION=28%结点个数（结点编码总数）NVFIX=2%受约束边界点数NFORCE=1%结点荷载个

7、数YOUNG=2e11%弹性模量POISS=0.25%泊松比THICK=0.1%厚度LNODS=[123;245;253;356;478;485;589;596;6910;71112;7128;81213;8139;91314;91410;101415;111617;111712;121718;121813;131819;131914;141920;142015;152021;162223;162317;172324;172418;182425;182519;251926;12192620;20

8、2627;202721;212728]%单元定义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=[00;-0.751.5;0.751.5;-1.53;03;1.53;-2.254.5;-0.754.5;0.754.5;2.254.5;-36;-1.56;06;1.56;36;-3.757.5;-2.257.5;-0.757.5;0.757.5;2.257.5;3.757.5;-4.59;-39;-1.59;09;1.59;39;4.59]%结点坐标数组%坐标：x，y坐标