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

《平面三角形单元有限元程序设计》由会员分享，可在线阅读，更多相关《平面三角形单元有限元程序设计（9页珍藏版）》请在装配图网上搜索。

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

2、有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。2）单元分析，建立单元刚度矩阵。3）整体分析，建立总刚矩阵。4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。6）解方程，求出节点位移。7）求出各单元的单元应力。8）计算结果整理。一、程序设计网格划分如图，将薄板如图划分为 6行，并建立坐标系，则1PPXXYY节点编号单元编号刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。通过循环逐

3、个计算：（1）每个单元对应 2 种单元刚度矩阵中的哪一种；（ 2）该单元对应总刚度矩阵的那几行哪几列（ 3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为 3 层循环：（1）最外层：逐行计算（ 2）中间层：该行逐个计算（ 3）最里层：区分为第 奇/偶 数个计算k1e66k1e 5656ke66ke565622kZe6 6 kZe 56 56单元刚度的集成：Kk1e k2ekZe 边界约束的处理：划 0置 1法2适用：这种方法适用于边界节点位移分量为已知(含为 0)的各种约束。做法 ：（ 1） 将总刚矩阵K中相应于已知位移行主对角线元素置 1，其他元素改为零；同时将载荷列阵R中相应元素

4、用已知位移置换。 这样，由该方程求得的此位移值一定等于已知量。（ 2） 将K中已知位移相应的列的非主对角成元素也置 0，以保持K的对称性。 当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献 已知位移和相应总刚元素的乘积。 若约束为零位移约束时，此步则可省去。特点 ：（1） 经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。（ 2） 但这种方法不改变方程阶数，利于存贮。（ 3） 不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。程序如下：变量说明NNODE单元节点数NPION总结点数NELEM单

5、元数NVFIX受约束边界点数FIXED约束信息数组NFORCE节点力数FORCE节点力数组COORD结构节点坐标数组LNODS单元定义数组YOUNG弹性模量POISS泊 松比THICK厚 度B单元应变矩阵(3*6)D单元弹性矩阵(3*3)S单元应力矩阵(3*6)A单元面积ESTIF单元刚度矩阵ASTIF总体刚度矩阵ASLOD总体荷载向量3ASDISP节点位移向量ELEDISP单元节点位移向量STRESS单元应力%*%初始化clearformat short e%设定输出类型clear%清除内存变量NELEM=36%单元个数（单元编码总数）NPION=28%结点个数（结点编码总数）NVFIX=2

6、%受约束边界点数NFORCE=1%结点荷载个数YOUNG=2e11%弹性模量POISS=0.25%泊松比THICK=0.1%厚 度LNODS=1 2 3;2 4 5;2 5 3;3 5 6;478;485;589;596;6 9 10;7 11 12;7 12 8;8 12 13;8 13 9;9 13 14;9 14 10;10 14 15;11 16 17;11 17 12; 12 17 18; 12 18 13;13 18 19; 13 19 14;14 19 20;14 20 15;15 20 21;16 22 23;16 23 17;17 23 24;17 24 18;18 24 2

7、5;18 25 19;25 19 26;19 26 20;20 26 27;20 27 21;21 27 28 %单 元定 义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3;1.5 3;-2.25 4.5;-0.75 4.5;0.75 4.5;2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6;-3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5;2.25 7.5;3.75 7.5;-4.5 9;-3 9;-1.5 9;0 9;1.5 9;3 9;4.5 9

8、%结 点 坐标 数 组%坐 标： x， y坐标 （共NPOIN 组 ）FORCE=1 0 -15%结 点 力数 组（ 受力 结点 编号 , x方向 ,y方向 ）FIXED=221 1;281 1%约 束信 息 （约 束点 ， x约束 ， y约 束）%有约束为1，无约束为0%*%生成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION);%生成特定大小总体刚度矩阵并置0%*for i=1:NELEM4%生成弹性矩阵 DD= 1POISS0;POISS10;00(1-POISS)/2*YOUNG/(1-POISS2)%*%计算当前单元的面积A=-det(1 COOR

9、D(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)/2 %* %生成应变矩阵 Bfor j=0:2b(j+1)=COORD(LNODS(i,(rem(j+1),3)+1),2) -COORD(LNODS(i,(re m(j+2),3)+1),2);c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(r em(j+2),3)+1),1);endB=b(1)

10、0b(2)0b(3)0;0c(1)0c(2)0c(3);c(1)b(1) c(2)b(2) c(3)b(3)/(2*A);B1( :,:,i)=B;%*%求应力矩阵S=D*BS=D*B;ESTIF=B*S*THICK*A;%求解单元刚度矩阵a=LNODS(i,:);%临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠

11、加到总刚%度矩阵中endendend%*%将约束信息加入总体刚度矩阵（对角元素改一法）5for i=1:NVFIXif FIXED(i,2)=1ASTIF(:,(FIXED(i,1)*2-1)=0;%一 列为 零ASTIF(FIXED(i,1)*2-1),:)=0;%一 行为 零ASTIF(FIXED(i,1)*2-1),(FIXED(i,1)*2-1)=1;%对角元素为1end%* %生成单元刚度矩阵并组成总体刚度矩阵 %*if FIXED(i,3)=1ASTIF( :,FIXED(i,1)*2)=0;%一列 为零ASTIF(FIXED(i,1)*2,:)=0;%一行为零ASTIF(FIXE

12、D(i,1)*2 ,FIXED(i,1)*2)=1;%对 角 元素为1endend%*%生成荷载向量ASLOD(1:2*NPION)=0;%总体荷载向量置零for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*%求解内力ASDISP=ASTIFASLOD%计算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量fori=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP%求内 力end（程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些地方数据输入时出了问题，还在寻找具体原因）6二、有限元软件分析设置材料参数建模7网格划分边界约束8添加载荷9