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

Matlab在线性立体有限元分析中的应用 摘要：Matlab具有强大的运算功能，本文以线性四面体元为例，详细介绍MATLAB在刚度矩阵推导，静力结构等有限元分析中的具体应用，编写了刚度矩阵，引用边界条件以及后处理各步骤的程序，该方法可以进一步推广到其他单元甚至更复杂的结构分析中。 关键词：Matlab 有限元 刚度矩阵 0 引言 Matlab是美国MathWork公司开发的用于数值计算，算法研究，建模仿真，实时实现的理想集成环境，因其完整的专业体系和强大的运算功能已广泛应用于工业、电子、信号处理、控制、建筑、教学等各个领域。 有限元是近代数值计算最有效方法之一．有限元法的基础是单元划分以及刚度矩阵的推导，目前，有限元分析已有一个相对固定的模式，而烦琐、复杂的矩阵运算、微分、积分是分析过程中的主要内容．通常，这种矩阵运算是由手工来完成的，工作量大，而且极易出错．利用MatLab丰富的符号运算功能，构建有限单元模型，完成刚度矩阵推导及后处理过程中的运算，不但速度快，而且准确性高。 利用Matlab编写函数M文件并在运算过程中调用，能够依据具体问题对模型进行分析运算，并能在类似问题中得到推广应用。 1 线性四面体有限元分析中的基本方程 线性四面体（立体）元（liner tetrahedral(solid)element）是既有局部坐标又有总体坐标的三维有限元，用线性函数描述。线性四面体元的系数有弹性模量E和泊松比，每个线性四面体与元有四个节点并且每个节点有三个自由度，如图1所示。这四个节点的总体坐标用、、、表示。单元刚度矩阵给定如下： （1.1） 式中V是单元的体积，由下式给出： （1.2） 图1 线性四面体（立体）元 矩阵由下式（1.3）确定： 在方程（1.3）中，形函数由下式给出： （1.4） 在方程（1.1）中，矩阵由下式（1.5）确定： 2 建立的Matlab函数 TetrahedronElementVolume(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4),该函数根据给出的第一个节点坐标，第二个节点坐标，第三个节点坐标和第四个节点坐标返回单元的体积。 TetrahedronElementStiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4),该函数用于计算弹性模量E、泊松比NU、第一个节点坐标，第二个节点坐标，第三个节点坐标和第四个节点坐标时的线性四面体（立体）单元的刚度矩阵。该函数返回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、第三个节点坐标、第四个节点坐标以及单元节点位移矢量为u时的单元应力矢量。该函数返回单元应力矢量。 QuadraticQuadElementStresses(sigma),该函数根据平面压力矢量sigma计算单元主应力。其形式为。其中sigma1、sigma2和sigma3为单元的主应力，该函数并不返回主应力方向角。 3 推广到线性多面体结构的有限元分析 显然，线性四面体元有12个自由度，每个节点有三个自由度。因此，对一个有n个节点的结构而言，其整体刚度矩阵K应该是3n×3n的（每个节点有三个自由度）。调用Matlab的TetrahedronAssenble函数可以得到整体刚度矩阵K，并列写以下方程组： （3.1） 式中，U是结构节点位移矢量，F是结构节点载荷矢量。在这一步中，边界条件被手动赋值给U和F。然后用分解法和高斯消去法求解方程组，得到位移和支反力，就可以用下式求出每个单元的应力矢量。 （3.2） 以线性六面体平板为例，平板受到均匀分布的载荷作用，将平板分为5个线性四面体单元，8个节点。 通过调用Matlab的TetrahedronElementStiffness函数，得到5个单元刚度矩阵，每个矩阵都是12