文章编号:1003-1251(2008)03-0045-03基于MATLAB的有限元程序研究
李 永,吴玉斌,倪庆杰(沈阳理工大学装备工程学院,辽宁沈阳110168)摘 要:有限元方法和MATLAB软件都是求解偏微分方程的有力工具,有限元方法的每个过程都是对方程进行矩阵的转化,而MATLAB是进行矩阵计算的工具,针对FOR2TRAN语言在面向对象编程方面的缺陷,利用MATLAB强大的矩阵运算能力和面向对象的交互式编程语言实现了有限元的数值计算过程,用两个算例进行了验证,相同节点位移比在1.0222～0.9750之间,满足工程需要,且MATLAB所需的代码比较少,可读性高.关键词:有限元;MATLAB;数值计算;程序中图分类号:TP311.11 文献标识码:A
AStudyontheFiniteElementProgramBasedonMATLAB
LIYong,WUYu2bin,NIQing2jie(ShenyangLigongUniversity,Shenyang110168,China)Abstract:FiniteelementmethodandMATLABsoftwarearepowerfultoolstosolvepartialdifferentialequations.Eachprocessoffiniteelementmethodistotransformtheequationstothematrix,andMATLABisatooltocalculatematrix.ConsideringthedeficienciesofFOR2TRANlanguageintheobject2orientedprogrammingseries,Usingthepowerfulcomputingca2pacityandobject2orientedinteractiveprogramminglanguageofMATLAB,theprocessoffi2niteelementnumericalcalculationisachieved,andtwoexamplesareusedtoverifyit.Thedisplacementratioofthesamenodesisbetween1.0222to0.9750.Itcanmeettheengi2neeringrequirments,andMATLABhasfewercodesandhighreadability.Keywords:finiteelement;MATLAB;numericalcalculation;program
5 有限元的MATLAB解法
1.打开MATLAB。
2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写），需要的话可点击Options菜单下Grid命令，打开栅格。
3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进行编辑并（如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标）
用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。
4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。
5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击 6 “Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。
6.进入PDE模式：点击"PDE"菜单下“PDE Mode”命令，进入PDE模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。
7.对模型进行剖分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。
8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。
.
. 有限元的MATLAB解法
1.打开MATLAB。
2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写），需要的话可点击Options菜单下Grid命令，打开栅格。
3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set
formula栏进行编辑并（如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标）
用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。
4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。
5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点.
. 击“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。
6.进入PDE模式：点击"PDE"菜单下“PDE Mode”命令，进入PDE模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。
7.对模型进行剖分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。
8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。
MATLAB报告
Matlab程序求解简要过程如下：
（1）求取单元节点位移提取矩阵T
单元节点位移提取矩阵T本质上是置换矩阵群中的一个，结果可将任意杂乱的节点顺序置换成统一的顺序。另一方面其作用是对单元刚度矩阵进行“升维操作”，将单元刚度矩阵统筹到整体刚度矩阵上来，便于对总体节点位移矩阵和支座反力进行求取。
本程序分析过程中对单元1的节点提取是按顺序编号1-2-3，对单元2的节点提取是按顺序编号2-3-4。单元1的节点位移提取矩阵如下： 单元2的节点位移提取矩阵如下：
（2）求取单元几何矩阵B
单元1的节点按编号顺序1-2-3分别进行对几何函数矩阵或算子矩阵的bi逆时针操作，对ci顺时针操作；单元2的节点按编号顺序2-3-4分别进行对几何函数矩阵的bi顺时针操作，对ci逆时针操作.在MATLAB程序中通过mod()取模函数来达到对节点的顺时针或逆时针循环操作。
单元1的几何矩阵如下：
单元2的几何矩阵如下：
（3）求取应力矩阵S
单元应力矩阵满足S=D*B,其中D为弹性矩阵，B为单元几何矩阵 各单元的弹性矩阵如下：
单元1的应力矩阵如下：
单元2的应力矩阵如下：
（4）求取单元刚度矩阵K
单元刚度矩阵K满足公式K=B’*D*B*t*A,其中t为平面板的厚度，A为单元面积，且单元刚度矩阵为对称矩阵。
单元1的刚度矩阵如下：
单元2的刚度矩阵如下：
（5）求取总体刚度矩阵sumKK
由上述步骤求得的单元刚度矩阵K利用单元虚功原理和刚度方程可导出K’*δ=f，其中δ为单元节点位移列阵，f为单元等效节点载荷列阵，为了能将各个单元刚度方程统一到一个整体，便需要步骤（1）的单元节点提取矩阵对单元刚度方程进行变换，将两个变换结果联立便得到总体刚度方程，其中也可得到总体刚度矩阵sumKK，且总体刚度矩阵可由sumKK=Σ T’*K*T求得。
总体刚度矩阵如下：
（6）求取总体节点位移矩阵和支座反力
利用上述步骤提到的总体刚度方程sumKK*delta=F，其中delta为总体节点位移矩阵，F为总体等效节点载荷列阵。利用总体边界条件可将总体刚度方程分成两部分，一部分先求取总体节点位移矩阵delta，最后代入总体刚度方程求取总体等效节点载荷列阵F，由此便得到支座反力。
