基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板.doc

《基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板.doc》由会员分享，可在线阅读，更多相关《基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板.doc（24页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、基于Matlab语言按平面三角形单元划分结构有限元程序设计专 业： 建筑及土木工程班 级： 建工研12-2姓 名： 韩志强学 号： 471220580基于Matlab语言按平面三角形单元划分结构有限元程序设计一、 有限单元发及Matlab语言概述1. 有限单元法随着现代工业、生产技术发展，不断要求设计高质量、高水平大型、复杂与精密机械及工程结构。为此目，人们必须预先通过有效计算手段，确切预测即将诞生机械与工程结构，在未来工作时所发生应力、应变与位移因此，需要寻求一种简单而又精确数值分析方法。有限单元法正是适应这种要求而产生与发展起来一种十分有效数值计算方法。有限元法把一个复杂结构分解成相对简单

2、“单元”，各单元之间通过结点相互连接。单元内物理量由单元结点上物理量按一定假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成结构。我们只要分析每个单元力学特性，然后按照有限元法规则把这些单元“拼装”成整体，就能够得到整体结构力学特性。有限单元法基本步骤如下：(1)结构离散：结构离散就是建立结构有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成有限元模型。在该步骤中，需要根据结构几何特性、载荷情况等确定单元体内任意一点位移插值函数。 (2)单元分析：根据弹性力学几何方程以及物理方程确定单元刚度矩阵。(3)整体分析:把各个单元按原来结构重新连接起来，并在单元刚

3、度矩阵基础上确定结构总刚度矩阵，形成如下式所示整体有限元线性方程：式中，是载荷矩阵, 是整体结构刚度矩阵, 是节点位移矩阵。 (4)载荷移置:根据静力等效原理，将载荷移置到相应节点上，形成节点载荷矩阵。(5)边界条件处理:对式所示有限元线性方程进行边界条件处理。 (6)求解线性方程:求解式所示有限元线性方程，得到节点位移。在该步骤中，若有限元模型节点越多，则线性方程数量就越多，随之有限元分析计算量也将越大。(7)求解单元应力及应变 根据求出节点位移求解单元应力与应变。 (8)结果处理及显示 进入有限元分析后处理部分，对计算出来结果进行加工处理，并以各种形式将计算结果显示出。2. Matlab简

4、介在用有限元法进行结构分析时，将会遇到大量数值计算，因而在实用上是一定要借助于计算机与有限元程序，才能完成这些复杂而繁重数值计算工作。而Matlab是当今国际科学界最具影响力与活力软件。它起源于矩阵运算，并已经发展成一种高度集成计算机语言。它提供了强大科学计算，灵活程序设计流程，高质量图形可视化及界面设计，便捷及其他程序与语言接口功能。Matlab在各国高校及研究单位起着重大作用。“工欲善其事，必先利其器”。如果有一种十分有效工具能解决在教学及研究中遇到问题，那么Matlab语言正是这样一种工具。它可以将使用者从繁琐、无谓底层编程中解放出来，把有限宝贵时间更多地花在解决问题中，这样无疑会提高工

5、作效率。 目前，Matlab已经成为国际上最流行科学及工程计算软件工具，现在Matlab已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景全新计算机高级编程语言了，有人称它为“第四代”计算机语言，它在国内外高校与研究部门正扮演着重要角色。Matlab语言功能也越来越强大，不断适应新要求提出新解决方法。可以预见，在科学运算、自动控制及科学绘图领域Matlab语言将长期保持其独一无二地位。为此，本例采用Matlab语言编程，以利用其快捷强大矩阵数值计算功能。二、 问题描述一矩形薄板，一边固定，承受150kN集中荷载，结构简图及按平面三角形单元划分有限元模型图如下所示。材料参数：弹性

6、模量；泊松比：；薄板厚度。在本例中，所取结构模型及数据主要用于程序设计理论分析，及工程实际无关。三、 参数输入：单元个数NELEM=4 节点个数NNODE=6受约束边界点数NVFIX=2 节点荷载个数NFORCE=1弹性模量YOUNG=2e8 泊松比POISS=0.2 厚度THICK=0.002 单元节点编码数组LNODS =节点坐标数组COORD =节点力数组FORCE =4 0 -150约束信息数组FIXED =以上数值数据为程序运行前输入初始数据，存为“”文本格式，同时必须放在Matlab工作目录下，路径不对程序不能自动读取指定初始文件，运行出错。初始数据文本文件输入格式如下图：四、 M

7、atlab语言程序源代码：1. 程序中变量说明第 24 页NNODE 单元节点数 NPION 总结点数NELEM 单元数 NVFIX 受约束边界点数FIXED 约束信息数组 NFORCE 节点力数FORCE 节点力数组 COORD 结构节点坐标数组 LNODS 单元定义数组YOUNG 弹性模量 POISS 泊松比THICK 厚度 B 单元应变矩阵(3*6) D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6)A 单元面积 ESTIF 单元刚度矩阵ASTIF 总体刚度矩阵 ASLOD 总体荷载向量 ASDISP 节点位移向量 ELEDISP 单元节点位移向量STRESS 单元应力 FP1 数据

8、文件指针2. Matlab语言程序代码%初始化及数据调用format short e %设定输出类型clear %清除内存变量FP1=fopen(,rt); %打开输入数据文件，读入控制数据NELEM=fscanf(FP1,%d,1), %单元个数（单元编码总数）NPION=fscanf(FP1,%d,1), %结点个数（结点编码总数）NVFIX=fscanf(FP1,%d,1) %受约束边界点数NFORCE=fscanf(FP1,%d,1), %结点荷载个数YOUNG=fscanf(FP1,%e,1), %弹性模量POISS=fscanf(FP1,%f,1), %泊松比 THICK=fsca

9、nf(FP1,%f,1) %厚度LNODS=fscanf(FP1,%d,3,NELEM) %单元定义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=fscanf(FP1,%f,2,NPION) %结点坐标数组 %坐标：x，y坐标（共NPOIN组）FORCE=fscanf(FP1,%f,3,NFORCE) %结点力数组（受力结点编号, x方向,y方向）FIXED=fscanf(FP1,%d,3,NVFIX) %约束信息（约束点，x约束，y约束） %有约束为1，无约束为0%生成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成

10、特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵D D= 1 POISS 0; POISS 1 0; 0 0 (1-POISS)/2*YOUNG/(1-POISS2)%计算当前单元面积A 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%生成应变矩阵B for j=0:2 b(j+1)= COORD(LNODS(i,(rem(j+1),3)+1),2)-CO

11、ORD(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*BS=D*B; ESTIF=B*S*THICK*A; %求解单元刚度矩阵 a=LNODS(i,:); %临时向量,用来记录当前单元节点编号 for j=1:3 fo

12、r 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:NVFIX if FIXED(i,2)=1 ASTIF(:,(FIXED(i,1)*2-1)=0; %一列为零 ASTIF(FIXED(i,1)*2-1),:)=0; %一行为零 ASTIF(FIXED(i,1)*

13、2-1),(FIXED(i,1)*2-1)=1; %对角元素为1 end if FIXED(i,3)=1 ASTIF( :,FIXED(i,1)*2)=0; %一列为零 ASTIF(FIXED(i,1)*2,:)=0; %一行为零 ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1 end end%生成荷载向量 ASLOD(1:2*NPION)=0; %总体荷载向量置零 for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end%求解内力 ASDISP=ASTIFASLOD %计

14、算节点位移向量 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 fclose(FP1); %关闭数据文件五、 程序运行结果NELEM = 4NPION = 6NVFIX = 2NFORCE = 1YOUNG = 200000000POISS =THICK =LNODS = 1 2 6 2 3 4 2 4 5 2 5

15、 6COORD = 0 0 1 0 2 0 2 1 1 1 0 1FORCE = 4 0 -150FIXED = 1 1 1 6 1 1D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =D = 2.0833

16、e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =ASDISP = 0 0 0 0（说明：以上12个ASDISP输出结果数据表示节点位移向量，即各节点分别在X，Y方向上产生位移。）i = 1STRESS = -1.6793e+005 -3.3586e+004 -1.3207e+005i = 2STRESS = -5.5503e+004 -5.5503e+004 -5.5503e+004i = 3STRESS = 5.5503e+004 -4.9526e+004 -9.4497e+004i = 4STRESS = 1.

17、6793e+005 -2.7040e+004 -1.7932e+004（说明：以上四组STRESS输出结果数据表示单元应力SX,SY,SXY，i为单元号。）六、 ANSYS建模比较分析利用ANSYS有限元分析软件，完全按照前面Matlab程序设计各项条件参数以及单元划分方式建立ANSYS模型，其求解结果及以上程序计算结果比较，验证程序是否正确。ANSYS模型节点单元如下图所示：ANSYS模型求解变形图如下所示：ANSYS求解节点位移结果列表显示如下：（单位：m）PRINT U NODAL SOLUTION PER NODE* POST1 NODAL DEGREE OF FREEDOM LIST

18、ING * THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM NODE UX UY UZ USUM 1 0.0000 0.0000 0.0000 0.0000 6 0.0000 0.0000 0.0000 0.0000 MAXIMUM ABSOLUTE VALUES NODE 4 4 0 4ANSYS求解单元应力结果列表显示如下：（单位：KN/m2）PRINT S ELEMENT SOLUTION PER ELEMENT* POST1 ELEMENT NODAL STRESS LISTING

19、* THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATES ELEMENT= 1 PLANE182 NODE SX SY SZ SXY SYZ 1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 2 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182 NODE SX SY SZ SXY SYZ 2 -55503. -

20、55503. 0.0000 -55503. 0.0000 3 -55503. -55503. 0.0000 -55503. 0.0000 4 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182 NODE SX SY SZ SXY SYZ 2 55503. -49526. 0.0000 -94497. 0.0000 4 55503. -49526. 0.0000 -94497. 0.0000 5 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182 NODE SX SY S

21、Z SXY SYZ 2 0.16793E+06 -27040. 0.0000 -17932. 0.0000 5 0.16793E+06 -27040. 0.0000 -17932. 0.0000 6 0.16793E+06 -27040. 0.0000 -17932. 0.0000 结论通过比较Matlab语言设计程序运行结果与ANSYS建模分析结果可知，两种方式算出结果完全一致，说程序设计正确。所以本程序适用于按三角形单元划分平面结构有限元分析。format short eclearFP1=fopen(LinearTriangleElement of Thin plate.txt,rt);N

22、ELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUNG=fscanf(FP1,%e,1)POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%f,1)LNODS=fscanf(FP1,%d,3,NELEM)COORD=fscanf(FP1,%f,2,NPION)FORCE=fscanf(FP1,%f,3,NFORCE)FIXED=fscanf(FP1,%d,3,NVFIX)ASTIF=zeros(2*NPION,2*NPION); %生成

23、特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS2)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2%计算当前单元面积AA=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%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNOD

24、S(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);endB=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*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)*

25、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);endendend%将约束信息加入总体刚度矩阵for 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;endif FIXED(i,3)=1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1

26、)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIFASLODELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISPendfclose(FP1);format short eclearFP

27、1=fopen(LinearTriangleElement of Thin plate.txt,rt);NELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUNG=fscanf(FP1,%e,1)POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%f,1)LNODS=fscanf(FP1,%d, NELEM,3)COORD=fscanf(FP1,%f, NPION,2)FORCE=fscanf(FP1,%f, NFORCE,3)FIXED

28、=fscanf(FP1,%d, NVFIX,3)ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS2)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2%计算当前单元面积AA=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%生成应变矩阵Bfor

29、 j=0:2b(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);endB=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*BS=D*B;ESTIF=B*S*THICK*A;a=LNODS(i,:);

30、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);endendend%将约束信息加入总体刚度矩阵for 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;endif FIXED(i,3)=1ASTIF(

31、:,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIFASLODELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISPendfclose(FP1);