Matlab中有限元脚本程序的编程

Matlab¹(华中科技大学武汉430074)介绍在MATLAB中进行有限元计算时脚本编程的方法。脚本编程方法比使用图形界面作图方法求解更加灵活,对复杂的边界条件处理更加容易控制,它扩展了MATLAB在有限元计算方面的应用范围。:MATLAB脚本有限元:TP301ProgrammingoftheFiniteElementMethodwithScriptinMATLABWangShuchengZhanQionghua(HuazhongUniversityofScience&Technology,Wuhan430074)Abstract:ThispaperintroducestheprogrammingofthefiniteelementmethodwithMATLAB.ThecomputationmethodbywritingscriptismoreagilethanoneofGraphicalUserInterfaceandthiscandisposethecomplicatedboundaryproblemeasily.TheapplicationofMATLABonthefiniteelementmethodisextended.Keywords:MATLAB,script,finiteelementmethodClassnumber:TP3011引言MATLAB是一种以矩阵为基本编程单位的数值计算工具[1-2],既有高效的科学计算功能,又有强大的图形处理功能,是工程应用中广泛使用的软件。MATLAB中自带了有限元的工具)))偏微分方程工具箱(PartialDifferentialEquationTOOL2BOX)。对于有限元的求解,MATLAB的帮助只是涉及了一些较简单的例子,对于复杂的边界的处理和一些总体的编程思路很少提及,并且对于数据结构(特别是几何文件)介绍的也不够。通过对电磁场和温度场的计算表明,如果要进行复杂的场的计算,必需使用脚本程序编程。本文针对这一情况,对编程思路和数据结构作了一些具体的介绍。2MATLAB中的有限元求解方法和方程类型一般来说,在MATLAB中有限元的计算方法分为两种。第一种:PDE图形用户界面(GraphicalUserInterface)的方法。在命令行窗口中输入pdetool,就会出现这个界面,然后在界面中作图求解。这种方法采用纯粹的作图求解,不用编写程序,只能求解一些简单的问题。第二种:调用命令行函数(Command-LineFunctions)的方法,对于复杂的有限元的求解可以通过编写命令行函数的脚本格式来实现,本文介绍的就是这种方法。MATLAB中能够求解的偏微分方程(PDE)类型(式1-3)及边界条件(式4、5)如下:eigenvablueproblem:-ý#(cýu)+au=Kdu(1)parabolicpde:d9u9t-ý#(cýu)+au=f(2)hyperbolicpde:d92u9t2-ý#(cýu)+au=f(3)Neumann:ny.(cªýu)+qu=g(4)dirichlet:hu=r(5)24计算机与数字工程第32卷¹收到本文时间:2004年2月5日3MATLAB中的求解过程Matlab帮助中给出了场的计算流程,如图(1),图中椭圆代表数据,方形代表函数,其过程说明如下:(1)准备所有的数据:两个文件(几何文件和边界文件)以及参数(c、a、f,如果是时变场的则还要d);(2)先对几何文件进行剖分,产生p(点)、e(线条)和t(三角形);(3)输入剖分的信息,外加边界条件(边界条件文件),参数(caf等)后进行求解,得到求解结果;(4)绘图,对于稳态场采用绘图函数就可以了,但是对于时变场,最好的表达方式还是用动画。图1有限元程序的流程4MATLAB中的数据格式对于编程来说,最复杂的莫过于其数据格式。只有了解了数据格式,才能对数据进行处理。对于MATLAB中有限元的运算中的数据如下:4.1几何文件包含了确定求解区域的几何分布的回调函数(被系统调用的函数),几何文件的作用就是把求解区域的几何图形通过三个步骤传给MATLAB系统。在MATLAB中的函数有这样的特点:传给被调用的函数的参数的个数可以和被调用的函数的参数的个数不同,被调用的函数可以自己判断传入的参数的个数。设回调函数的名为SRMG(可以任意命名),对于用户编写的几何回调函数来说,三个步骤为:(1)返回边界总数,此时MATLAB系统输入的输入参数个数为0,函数的形式为:NE=SRMG()其中:NE给出几何边界的数目。(2)返回系统指定的边界信息,此时MATLAB系统输入参数为1个,函数的形式为:D=SRMG(BS)其中:BS系统输入的边界编号;D为对应边界上面的参数每条边界占一列,在每一列中格式如下:第1行包含起始参数;第2行包含中止参数;第3行包含边界的左手区域编号;第4行包含边界的右手区域编号;(3)返回边界上面具体的点,此时MATLAB系统的输入参数为2个,函数的形式为:[X,Y]=SRMG(BS,S)其中:BS为边界的编号,S为系统给定的中间参数,其值介于起始参数和中止参数之间,X,Y为对应边界对应参数上的点。举例说明:(1)对于BS对应的一条边界y=2x,xI[0,1],如果输入参数仅为BS,则返回起始参数0,1和左右手区域的编号。如果输入参数为BS和S,则返回x=Sy=2S(2)对于BS对应的一条边界圆r=2,HI[0,2pi],如果输入参数仅为BS,且对应其边界则返回起始参数0,2pi和左右手区域的编号。如果输入参数为BS和S,则返回x=2cos(S)y=2sin(S)通过上述步骤,完成了对图(1)中的边界文件的准备。4.2边界文件包含了边界的约束条件的矩阵b1。b1的列数等于边界总数,每一列包含了每一条边界的边界条件的信息。每一列的格式具体说来如下:第1行为Neumann边界条件的维数N,即式(4)的维数;第2行是Dirichlet边界条件的维数M,即式(5)的维数;第3行是3+N2-1行是对应的式(4)中q字符串的长度;第3+N2行到3+N2+N-1行对应式(4)中g字符串的长度;第3+N2+N行到3+N2+N+MN-1行对应式(5)中h字符串的长度;25第32卷(2004)第6期计算机与数字工程第3+N2+N+NM行到3+N2+N+MN+M-1行对应式(5)中r字符串的长度;上面只是字符串的长度,接下来就是对应的式(4)和式(5)中的q、g、h、r的字符串的表达式。举例说明:对于边界条件ny.(cªýu)=-x2,在MAT2LAB中N=1,M=0,q的表达式为-0.,长度q为1,g的表达式为--x.C2.,长度为5,所以其对应的边界条件的矢量写法为[1015-0.--x.C2.].。4.3c、a、fc、a、f是偏微分方程的参数,见式(1)、(2)和(3),其表达式是以字符串的形式给系统。对于不同的区域的偏微分方程的参数c、a、f不相同,采用/!0分开。例如:c对于不同的三个区域的参数分别为50,60和70则表达式为[-50!60!70.]。4.4p、e、tp为剖分点,p矩阵为2行,分别包含了剖分点的x和y坐标;e为剖分线段,e矩阵为7行,分别包含了起点终点的编号等信息;t为剖分三角形,t矩阵为4行,包含了三角形三个顶点的编号和所在区域的编号。4.5uu仅有一行,为对应的剖分点p的求解结果。5采用脚本的高级编程5.1上面介绍了有限元求解的数据格式,我们可以按照这些数据格式编写一些比较复杂场的求解:(1)关于几何边界,如采用第一种求解方式)))图形界面,只能求解由圆、椭圆、直线边界组成的简单的几何区域,但是对于其它的复杂的几何边界则无法实现。如采用第二种方式)))编写命令行函数的脚本程序,可以通过回调函数的形式,写出边界曲线的函数,被系统回调。这样只要知道曲线的函数就可以求解。(2)边界的控制。对于运动物体的求解(比如运动的铁心的磁场),对于一些边界的变换(比如不同的电机参数,比如电机的外径不同)。这些如果采用图形界面求解方式,则当不同情况、不同状态时,必须重新绘制图形。如果采用脚本编写程序,在几何文件中,只要改写边界参数,就可达到改变几何边界的目的。(3)耦合边界条件。当两个边界条件发生耦合关系,例如两条边界的热流量的和为零,可以通过假设边界上面的温度,来迭代求解,直到满足耦合的约束条件。5.2作者编写了开关磁阻电动机的温度分布场的计算程序,功能为在用户界面中输入电机参数(包含结构尺寸、材料导热系数、通风情况等),就能算出并显示温度场,而且找出最高温度点的值及其坐标。其中计算过程中有边界的耦合过程,需要迭代来求解。求解的几何区域和结果见图(2),根据求解区域的对称性,图中的求解区域为1/8电机的定子截面。图2求解区域与求解结果程序的流程图如图(3),图中实线表示数据的流向,虚线表示循环的控制,流程的大致方向为从上往下。Workspace是MATLAB中特有的,用来存储或者读取变量。主程序是对界面消息响应函数,从而获取输入的电机参数。由于几何文件和边界文件的格式是特定的,不能从主程序中获得数据,这里就采用了Workspace来间接获取主程序数据,这也是脚本编程的关键。26Matlab中有限元脚本程序的编程第32卷图3求解温度场的流程[1]蒙以正.MATLAB5.X应用与技巧.北京:科学出版社,1999.8[2]刘卫国.科学计算与MATLAB语言.中国铁道出版社,2000.4(4)[10]W.Landi,etal.Pointer-inducedAliasing:AproblemClassification.Conf.Rec.18thAnn.ACMSymp.Prin2ciplesofProgrammingLanguages.1991,93-103[11]W.Landi,etal.Undecidabilityofstaticanalysis.ACMLett.ProgrammingLanguagesSyst.,vol.1,Dec.1992,323-337[12]AckosKiss,etal.InterproceduralStaticSlicingofBinaryExecutables.ThirdIEEEInternationalWorkshoponSourceCodeAnalysisandManipulation.September26-27,2003,118-127[13]M.Harman,etal.AnInterproceduralAmorphousSlicerforWSL.SecondIEEEInternationalWorkshoponSourceCodeAnalysisandManipulation(SCAM.02).October01-01,2002,105-114[14]黄文伟.C程序文件间依赖性分析[硕士学位论文].东南大学,2004[15]张挺.基于依赖性分析的Ada程序切片[硕士学位论文].东南大学,1997[16]陈振强.基于依赖性分析的程序切片技术研究[博士学位论文].东南大学,2003[17]周毓明.软件度量中的若干问题研究[博士学位论文].东南大学,200227第32卷(2004)第6期计算机与数字工程