韩志强-有限元程序设计

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

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

3、程确定单元的刚度矩阵。(3) 整体分析 : 把各个单元按原来的结构重新连接起来， 并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程：FK式中，F 是载荷矩阵 ,K 是整体结构的刚度矩阵,是节点位移矩阵。(4) 载荷移置 : 根据静力等效原理， 将载荷移置到相应的节点上， 形成节点载荷矩阵。(5) 边界条件处理 : 对式所示的有限元线性方程进行边界条件处理。(6) 求解线性方程 : 求解式所示的有限元线性方程， 得到节点的位移。 在该步骤中，若有限元模型的节点越多， 则线性方程的数量就越多， 随之有限元分析的计算量也将越大。(7) 求解单元应力及应变 ：根据求出的

4、节点位移求解单元的应力和应变。.(8) 结果处理与显示 ：在进行程序设计时，还要进入有限元分析的后处理部分，对计算出来的结果进行加工处理，并以各种形式将计算结果显示出。2. Matlab 简介在用有限元法进行结构分析时， 将会遇到大量的数值计算， 因而在实用上是一定要借助于计算机和有限元程序，才能完成这些复杂而繁重的数值计算工作。而 Matlab 是当今国际科学界最具影响力和活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。 它提供了强大的科学计算， 灵活的程序设计流程，高质量的图形可视化与界面设计， 便捷的与其他程序和语言接口的功能。 Matlab 在各国高校与研究单位起着

5、重大的作用。“工欲善其事，必先利其器” 。如果有一种十分有效的工具能解决在教学与研究中遇到的问题，那么 Matlab 语言正是这样的一种工具。它可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作效率。 目前， Matlab 已经成为国际上最流行的科学与工程计算的软件工具，现在的 Matlab 已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了， 有人称它为 “第四代”计算机语言，它在国内外高校和研究部门正扮演着重要的角色。 Matlab 语言的功能也越来越强大， 不断适应新的要求提出新的解决方法。

6、可以预见，在科学运算、自动控制与科学绘图领域 Matlab 语言将长期保持其独一无二的地位。为此，本例采用 Matlab 语言编程，以利用其快捷强大的矩阵数值计算功能。二、 问题描述一矩形薄板， 一边固定，承受 150kN 集中荷载，结构简图及按平面三角形单元划分的有限元模型图如下所示。材料参数：弹性模量E2108 KN / m 2 ；泊松比：0.2 ；薄板厚度 2mm。在本例中，所取结构模型及数据主要用于程序设计理论分析，与工程实际无关。.F=150kN1m2myx654431m121231m1m三、 参数输入：单元个数 NELEM=4节点个数 NNODE=6受约束边界点数NVFIX=2节点

7、荷载个数 NFORCE=1弹性模量 YOUNG=2e8泊松比 POISS=0.2厚度 THICK=0.0021 2 6单元节点编码数组LNODS =23 424 525 6.0 01 0节点坐标数组 COORD =202 11 10 1节点力数组 FORCE =40-1501 1 1约束信息数组 FIXED =6 1 1以上数值数据为程序运行前输入的初始数据，存为“471220580.txt ”文本格式，同时必须放在 Matlab 工作目录下，路径不对程序不能自动读取指定初始文件，运行出错。初始数据文本文件输入格式如下图：.四、 Matlab语言程序源代码：1. 程序中变量说明NNODE单元节

8、点数B单元应变矩阵 (3*6)NPION总结点数D单元弹性矩阵 (3*3)NELEM单元数S单元应力矩阵 (3*6)NVFIX受约束边界点数A单元面积FIXED约束信息数组ESTIF单元刚度矩阵NFORCE节点力数ASTIF总体刚度矩阵FORCE节点力数组ASLOD总体荷载向量COORD结构节点坐标数组ASDISP节点位移向量LNODS单元定义数组ELEDISP单元节点位移向量YOUNG弹性模量STRESS单元应力POISS泊松比FP1数据文件指针THICK厚度2. Matlab语言程序代码%*%初始化及数据调用format short e%设定输出类型clear%清除内存变量FP1=fope

9、n( 471220580.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,3,NELEM)%单元定义数组（单元结点号）%相应为单元结点

10、号（编码） 、按逆时针顺序输入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);%生成特定大小总体刚度矩阵并置0%*.for i=1:NELEM%生成弹性矩阵 DD= 1POISS0;POISS10

11、;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%*%生成应变矩阵 Bfor 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

12、+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1);endB=b(1)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,(

13、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;%对角元素为1end.if FIXED(i,3)=1ASTIF( :,FIXED(i,1)*2)=0;%一列为零ASTI

14、F(FIXED(i,1)*2,:)=0;%一行为零ASTIF(FIXED(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:LN

15、ODS(i,j)*2);%取出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP%求内力endfclose(FP1);%关闭数据文件五、 程序运行结果NELEM =4NPION =6NVFIX =2NFORCE =1YOUNG =200000000POISS =2.0000e-001THICK =2.0000e-003LNODS =126234.245256COORD =0 01 02 02 11 10 1FORCE =4 0 -150FIXED =111611D =2.0833e+0084.1667e+00704.1667e+0072.0833e+0080

16、008.3333e+007A =5.0000e-001D =2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A =5.0000e-001D =2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A =5.0000e-001D =2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A =5.0000e-001ASDISP =00-8.0607e-004-1.5848e-003-1.0281e-00

17、3.-4.4727e-0031.1937e-003-4.6947e-0038.6670e-004-1.8880e-00300（说明：以上 12 个 ASDISP输出结果数据表示节点位移向量， 即各节点分别在 X， Y 方向上产生的位移。单位： m）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.6793e+005-2.7040e+0

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

19、M LISTING *THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEMNODEUXUYUZUSUM10.00000.00000.00000.00002 -0.80607E-03-0.15848E-020.00000.17780E-023 -0.10281E-02-0.44727E-020.00000.45893E-0240.11937E-02-0.46947E-020.00000.48441E-0250.86670E-03-0.18880E-020.00000.20774E-0260.00000.00000.00000.0000MAXIMUM ABSOLUTE VALUESNODE4404VALUE0.11937E-02-0.46947E-02 0.0000 0.48441E-02.ANSYS 求解单元应力结果列表显示如下： （单位： KN/m 2 ）PRINT SELEMENT SOLUTION PER ELEMENT* POST1 ELEMENT NODAL STRESS LISTING *THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATESELEMENT=1PLANE182NOD