MATLAB在线性四面体有限元分析中的应用

Matlab在线性立体有限元分析中的应用摘要：Matlab具有强大的运算功能，本文以线性四面体元为例，详细介绍MATLAB在刚度矩阵推导，静力结构等有限元分析中的具体应用，编写了刚度矩阵，引用边界条件以及后处理各步骤的程序，该方法可以进一步推广到其他单元甚至更复杂的结构分析中。关键词：Matlab有限元刚度矩阵0引言Matlab是美国MathWork公司开发的用于数值计算，算法研究，建模仿真，实时实现的理想集成环境，因其完整的专业体系和强大的运算功能已广泛应用于工业、电子、信号处理、控制、建筑、教学等各个领域。有限元是近代数值计算最有效方法之一．有限元法的基础是单元划分以及刚度矩阵的推导，目前，有限元分析已有一个相对固定的模式，而烦琐、复杂的矩阵运算、微分、积分是分析过程中的主要内容．通常，这种矩阵运算是由手工来完成的，工作量大，而且极易出错．利用MatLab丰富的符号运算功能，构建有限单元模型，完成刚度矩阵推导及后处理过程中的运算，不但速度快，而且准确性高。利用Matlab编写函数M文件并在运算过程中调用，能够依据具体问题对模型进行分析运算，并能在类似问题中得到推广应用。1线性四面体有限元分析中的基本方程线性四面体（立体）元（linertetrahedral(solid)element）是既有局部坐标又有总体坐标的三维有限元，用线性函数描述。线性四面体元的系数有弹性模量E和泊松比，每个线性四面体与元有四个节点并且每个节点有三个自由度，如图1所示。这四个节点的总体坐标用111x（,y,z）、222x（,y,z）、333x（,y,z）、444x（,y,z）表示。单元刚度矩阵给定如下：[][][][]TkVBDB（1.1）式中V是单元的体积，由下式给出：11122233344411611xyzxyzVxyzxyz（1.2）图1线性四面体（立体）元矩阵B由下式（1.3）确定：3124312431243311224433112244311220000000000000000000000000000000000xxxxyyyyzzzzyxyxyxyxzyzyzyzyzxzxNNNNNNNNNNNNBNNNNNNNNNNNNNNNNNNNNN34400zxzxNNN在方程（1.3）中，形函数由下式给出：111111()6NxyzV222221()6NxyzV333331()6NxyzV444441()6NxyzV（1.4）在方程（1.1）中，矩阵D由下式（1.5）确定：10001000100012000002(1)(12)1200000212000002ED2建立的Matlab函数TetrahedronElementVolume(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4),该函数根据给出的第一个节点坐标111(,,)xyz，第二个节点坐标222(,,)xyz，第三个节点坐标333(,,)xyz和第四个节点坐标444(,,)xyz返回单元的体积。TetrahedronElementStiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4),该函数用于计算弹性模量E、泊松比NU、第一个节点坐标111(,,)xyz，第二个节点坐标222(,,)xyz，第三个节点坐标333(,,)xyz和第四个节点坐标444(,,)xyz时的线性四面体（立体）单元的刚度矩阵。该函数返回12×12的单元刚度矩阵k.TetrahedronAssenble(K,k,i,j,m,n),该函数将连接节点i、节点j、节点m和节点n的线性四面体元的单元刚度矩阵k集成到整体刚度矩阵K。每集成一个单元，该函数都会返回3n×3n的整体刚度矩阵K。TetrahedronElementStresses(E,NU,x1,y1,x2,y2,x3,y3,x4,y4,u),该函数计算在弹性模量E、泊松比为NU、第三个节点坐标333(,,)xyz、第四个节点坐标444(,,)xyz以及单元节点位移矢量为u时的单元应力矢量。该函数返回单元应力矢量。QuadraticQuadElementStresses(sigma),该函数根据平面压力矢量sigma计算单元主应力。其形式为sigma1,sigma2,sigma3T。其中sigma1、sigma2和sigma3为单元的主应力，该函数并不返回主应力方向角。3推广到线性多面体结构的有限元分析显然，线性四面体元有12个自由度，每个节点有三个自由度。因此，对一个有n个节点的结构而言，其整体刚度矩阵K应该是3n×3n的（每个节点有三个自由度）。调用Matlab的TetrahedronAssenble函数可以得到整体刚度矩阵K，并列写以下方程组：{}{}KUF（3.1）式中，U是结构节点位移矢量，F是结构节点载荷矢量。在这一步中，边界条件被手动赋值给U和F。然后用分解法和高斯消去法求解方程组，得到位移和支反力，就可以用下式求出每个单元的应力矢量。{}[][]{}DBu（3.2）以线性六面体平板为例，平板受到均匀分布的载荷作用，将平板分为5个线性四面体单元，8个节点。通过调用Matlab的TetrahedronElementStiffness函数，得到5个单元刚度矩阵，每个矩阵都是12×12的。由于结构有8个节点，所以整体刚度矩阵是24×24的。由于在该结构中有5个这样的单元，所以5次调用Matlab的TetrahedronAssenble函数，就可以得到整体刚度矩阵K。将边界条件带入到该结构的方程组中，就可以进行求解计算。提取边界条件对应的子矩阵，就可以求出节点的支反力；建立节点位移矢量U，计算结构节点力矢量F。接下来，调用Matlab的TetrahedronElementStresses函数就能够计算出每个单元的主应力。4结论以线性四面体作为有限元模型的基本单元，可以推广到更为复杂的立体结构进行有限元分析。Matlab在构建模型，推导刚度矩阵，完成后处理计算中发挥了重要作用，构建了分析过程函数，算法准确可靠，并能在类似工程问题中重复利用，简便有效。参考文献[1]黄以实.用Matlab分析静力结构有限元法[J].广西大学学报：自然科学版，2009,7:293-298.[2]李伟旭，刘艳珂，张永水.用Matlab进行结构的有限元法分析[J].重庆交通学院学报,2007,2：23-29.[3]周水兴，陈山林.Matlab在有限元刚度矩阵推导中的应用[J].重庆交通学院学报:2007,4:29-31.[4]P.I.Kattan著，韩来彬译.MATLAB有限元分析与应用[M].北京：清华大学出版社，2004.