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

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

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

2、料性质、荷载条件、约束条件等，离散构造并进展单元编码、结点编码、结点位移编码、选取坐标系。2单元分析，建立单元刚度矩阵。3整体分析，建立总刚矩阵。4建立整体构造的等效节点荷载和总荷载矩阵5边界条件处理。6解方程，求出节点位移。7求出各单元的单元应力。8计算结果整理。一、程序设计网格划分如图，将薄板如图划分为6行，并建立坐标系，则PP*YY单元编号节点编号刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。通过循环逐个计算：1每个单元对应2种单元刚度矩阵中的哪一种；2该单元对应总刚度矩阵的那几行哪几列3

3、将该单元的单元刚度矩阵参加总刚度矩阵的对应行列循环又分为3层循环：1最外层：逐行计算2中间层：该行逐个计算3最里层：区分为第奇/偶数个计算单元刚度的集成：边界约束的处理：划0置1法适用：这种方法适用于边界节点位移分量为(含为0)的各种约束。做法：1将总刚矩阵K中相应于位移行主对角线元素置1，其他元素改为零；同时将载荷列阵R中相应元素用位移置换。这样，由该方程求得的此位移值一定等于量。2将K中位移相应的列的非主对角成元素也置0，以保持K的对称性。当然，在位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去位移对该方程的奉献位移和相应总刚元素的乘积。假设约束为零位

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

5、ASTIF 总体刚度矩阵ASLOD 总体荷载向量ASDISP 节点位移向量ELEDISP 单元节点位移向量STRESS 单元应力%*%初始化clearformat short e %设定输出类型clear %去除内存变量NELEM=36 %单元个数单元编码总数NPION=28 %结点个数结点编码总数NVFI*=2 %受约束边界点数NFORCE=1 %结点荷载个数YOUNG=2e11 %弹性模量POISS=0.25 %泊松比THICK=0.1 %厚度LNODS=1 2 3;2 4 5;2 5 3;3 5 6; 4 7 8;4 8 5;5 8 9;5 9 6; 6 9 10;7 11 12;7 1

6、2 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 25;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

7、;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 %结点坐标数组%坐标：*，y 坐标共 NPOIN 组FORCE=1 0 -15 %结点力数组受力结点编号, * 方向,y 方向FI*ED=22 1 1;28 1 1 %约束信息约束点，* 约束，y 约束 %有约束为 1，无约束为 0 %*%生成单

8、元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置 0 %*for i=1:NELEM %生成弹性矩阵 D D= 1 POISS 0; POISS 1 0; 0 0 (1-POISS)/2*YOUNG/(1-POISS2) %*%计算当前单元的面积 A=-det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2); 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2); 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)/2 %*%生

9、成应变矩阵 B for j=0:2 b(j+1)=COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,(rem(j+2),3)+1),2); c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1); end B=b(1) 0 b(2) 0 b(3) 0; 0 c(1) 0 c(2) 0 c(3); c(1) b(1) c(2) b(2) c(3) b(3)/(2*A); B1( :,:,i)=B;%* %求应力矩阵 S=D*B S=D*B; ESTIF=B*S*

10、THICK*A; %求解单元刚度矩阵 a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号 for j=1:3 for k=1:3 ASTIF(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); %根据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中 end end end %*%将约束信息参加总体刚度矩阵对角元素改一法 for i=1:NVFI* if FI*ED(i,2)=1 ASTIF(:,(FI*ED(i,1)*

11、2-1)=0; %一列为零 ASTIF(FI*ED(i,1)*2-1),:)=0; %一行为零 ASTIF(FI*ED(i,1)*2-1),(FI*ED(i,1)*2-1)=1; %对角元素为 1 end %*%生成单元刚度矩阵并组成总体刚度矩阵%* if FI*ED(i,3)=1 ASTIF( :,FI*ED(i,1)*2)=0; %一列为零 ASTIF(FI*ED(i,1)*2,:)=0; %一行为零 ASTIF(FI*ED(i,1)*2 ,FI*ED(i,1)*2)=1; %对角元素为 1 end end %*%生成荷载向量 ASLOD(1:2*NPION)=0; %总体荷载向量置零 f

12、or i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end %* %求解内力 ASDISP=ASTIFASLOD %计算节点位移向量 ELEDISP(1:6)=0; %当前单元节点位移向量 for i=1:NELEM for j=1:3 ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); %取出当前单元的节点位移向量 end i STRESS=D*B1(:, :, i)*ELEDISP %求内力 end 程序计算结果和有限元软件得出的结果稍有偏差，可能是程序*些地方数据输入时出了问题，还在寻找具体原因二、 有限元软件分析建模设置材料参数添加载荷边界约束网格划分. z.