COMSOLMULTIPHYSICS和数值分析基础

第一章COMSOLMＵＬTＩPHＹＳICS及数值分析基础W。B．Ｊ．ZIMMERＭAN,B.Ｎ.HEWAKAＮDＡＭBＹDepａrtmentｏfCｈemicalandＰrｏｃessＥnｇinｅｅriｎg，UｎiｖersityoｆSｈｅfｆield，NｅwｃasｔｌeＳｔrｅet，ＳｈefｆｉeldＳ13JDUniｔeｄKiｎgｄomＥ—ｍail：ｗ。ｚｉｍｍeｒman@shｅf.aｃ．uk本章主要介绍ＣOＭＳOＬＭuｌtipｈｙsｉｃｓ在零维和一维模型数值分析方面的几个关键内容。这些内容包括求根、步进式数值积分、常微分方程数值积分和线性系统分析。这几乎是全部的化工过程数学分析方法。下面通过CＯＭSOLMｕltｉpｈｙsｉcｓ中的一些常见化工过程应用实例来介绍这些方法，包括：闪蒸、管式反应器设计、集中反应系统和固体中热传导。1．简介本章内容很多，可以分为几个不同的目标。首先介绍了CＯMSOＬMｕltｉpｈysiｃs的主要工作特性;其次介绍了如何使用这些特性来分析一些简洁的，位于零维空间、一维空间或“空间－时间"系统中的化工问题.本章还盼望通过展现COMＳOLMultiphｙｓiｃs和ＭＡＴＬAＢ工具在化工过程分析中的强大功能,激发读者对使用COMSOLＭultｉphｙsicｓ进行建模与仿真的爱好。由于CＯＭSOLMultｉphyｓｉｃｓ不是一个通用的问题求解工具，所以一些目标需要迂回实现。作者在使用FＯRＴRAN、Maｔｈｅｍａtｉca和MATＬＡB解决化工问题方面有着丰富的教学阅历，并用这些工具实现过这里全部的例子。而且,扩展化工问题的数值分析也已经在PＯLYMＡTH［1]中实现，这似乎只在化工委员会的CAＣHE项目中使用过。本书前一版已经介绍过在零维空间中求解非线性代数方程和与时间有关的常微分方程的内容.从概念上讲,零维域就是一个简洁的有限元。通过讨论某一特定有限元中的变化对理解有限元方法格外有用。但是,CＯMＳＯLMｕｌｔｉphyｓicｓ通过独立对话框设置，使得零维几何方程和与时间相关的常微分方程求解变得格外简洁。所以本章将同时采纳这两种方法求解这些例子。2．方法１:求根典型的数值分析课程会讲解多种求根方法，但是从实际阅历来看,只有两种算法格外有用——二分法和牛顿法。我们这里没有列出全部方法,而是重点考虑为什么求根是最有效的数值分析工具。在线性系统中求根格外简洁，但是对于非线性系统这就是一个挑战,而全部感爱好的动力学问题几乎都是非线性系统.对非线性系统的求根起源于对反函数的描述。为什么呢?由于对于大多数非线性函数,“正向"y=ｆ(ｕ）很好表示，但是它的反函数u＝f－-1(ｙ)可能不能显式表示、多值（无意义）或根本不存在。如果反函数存在的话,求解反函数其实就是求根的过程——求解满意Ｆ（u)＝0的ｕ等价于求解F(ｕ）=f(u)－y＝0。由于大多数数值分析的目标是在系统约束下计算求解，所以这也等价于对全部的约束取反.COＭSOLMｕｌtｉｐhysｉcｓ拥有求解非线性问题的核心函数—-ｆｅｍnｌiｎ，本节主要介绍用它求解零维非线性问题。feｍnlｉn函数使用牛顿方法求解,由于只有一个变量u,牛顿法通过对一阶倒数迭代来求根。该方法首先估量函数的斜率范围，然后再逼近根。该斜率可以通过理论分析(牛顿—拉夫逊方法）和数值（正割法)方法求得。如果能用任何一种方法求得斜率，就可以用泰勒定律来逼近根。其基本思想就是使用目前推测值u０的泰勒展开式：ﻩﻩﻩﻩﻩﻩ(1）该公式可以化简，忽视(ｕ—ｕ0）的高阶项,计算根如下:ﻩﻩ ﻩﻩﻩﻩﻩ(2)这个方法可以快速地扩展到多维求解空间,例如将ｕ看作未知矢量，“被除”看作“乘以f的雅克比矩阵的逆”。下一节介绍ＣＯMＳOLMultiphyｓiｃs中的求根过程。2。1求根:COMＳOＬMulｔiphysｉcｓ非线性求解器的应用实例如上节所述,求根本身是一个“零维”活动，至少对于“空间－时间”系统多维未知矢量ｕ来说是这样的。ＣOMSＯＬ多物理场没有零维模式,所以我们临时采纳一维模式。这在方面增加了我们不需要的冗余功能。但是由于问题规模较小,ＣOMＳＯLＭｕlｔiphyｓｉcｓ编码效率高，且现代微处理器的运算速度快,这点就不成为问题了。启动ＭＡTLＡB并在命令窗口键入COMＳOLMulｔｉｐｈysicｓ。屏幕闪烁几秒后，会消灭一个模型导航窗口。依据表1所示步骤,建立一个零维应用模式来解决非线性多项式方程：ﻩﻩﻩﻩﻩﻩﻩ (3）通过在“Phｙsiｃs”菜单、“Subｄomａinｓｅttingｓ"选项中的设定，使得表１中的每个子域都满意该方程.注意左上角,这是以矢量符号给出的方程。在一维模式下,这个方程可以简化为:ﻩ ﻩﻩ (4）显然，α、γ和β在一维简化里是多余的。既然我们是求零维的根,全部公式４左侧的系数都可以设为0。通过重新组合多项式,我们发现ａ=4和ｆ＝u＋ｕ＋２。注意我们是通过设定最大基元尺寸为1、将域离散化为只有一个基元的域,从而得到零维空间的！表1系数模式下的求根实例。文件名称:ｒｏotｆiｎdｅr。mｐh．MoｄeｌNaviｇator选择1D空间维数，COMＳOLMulｔiｐhｙsｉcs：PDEmｏdes：ＰＤE，cｏeffｉcienｔ形式Ｄraw菜单Sｐeｃiｆyoｂjeｃｔｓ:Ｌｉnｅ。Ｃoordｉnａtｅs弹出菜单。x:01名称:inｔeｒｖａl,完成Pｈysics菜单：Bouｎdａrｙsｅttｉｎgs选择域：1和2（按住Ctrl键同时选中)选择Nｅumanｎ边界条件采纳默认设置q=０ｇ＝0,完成Physiｃs菜单:Ｓｕbｄoｍａinsｅttings选择域：1设定:c＝0;a=4;ｆ＝ｕ^3+u^2+2；＝0应用，选择Init选项卡：设定ｕ（ｔ０)=－2,完成Meｓh菜单：Meｓhpａraｍeters设定最大基元尺寸1点击remｅsｈ，完成Sｏｌve菜单：Soｌｖeｒpaｒａmeｔeｒs求解器选择Stationaryｎoｎliｎeａr,求解,完成Pｏｓt－prｏcｅssｉｎgPointEｖalutｉon选择边界1，表达式:u，完成设定初始推测值为ｕ（ｔ）＝－2，下面我们来寻找最接近这个值的根.你可能对为什么设置ａ=４而不是把全部有关参数放进ｆ中感到奇怪，那是由于对公式(4)右侧有限元离散后得不到一个奇异刚度矩阵。后处理阶段在输出窗口中显示出结果:Ｖalue：-2.7320５1，Expreｓｓioｎ:u，Bounｄａry:1。解析解得出的根为-1—，数值解能够很好的与其吻合。依据代数中二次公式的解知道，另外一个根是—1＋，第三个根是1.回到子域设置中来，如果最初设置的推测值为u(t0）＝-0．5,则ＣOMSOLMultiｐhyｓics解出u=0．7６2051,又一个很好的近似数值解。如果u(t0)=１.2作为最初猜想时，得到ｕ＝１。该练习体现了非线性求解器和问题的两个特性——（１）非线性问题可以有多个解；(２）最初推测值对于求解很关键。牛顿法通常收敛到最接近的解,但是在高阶非线性问题中,这点可能不成立。这些特性在多维求解空间和“空间-时间”相关系统中同样适用。COMＳOＬＭultiｐhysiｃｓ模型ｍｐｈ文件（ｒｏotfindｅｒ.ｍph)中包含了全部的MATLAB／COMSOＬSｃrｉpt源代码和FEMLAB的扩展功能,以便再次恢复到ＦEMLAB当前图形用户界面。该文件在网页HYPERLIＮＫ”htｔp：/／ｅｙｒie。shef.ac.uｋ／femｌaｂ＂ｈｔtp：//eyriｅ。sheｆ.ａc.uｋ/ｆeｍlａb中可以得到.只需打开“Menu”菜单,选择“Opｅｎmｏｄeｌｍ-fiｌe”选项，在弹出的“Opｅnfile”对话框中选中它,你就可以在“Ｓubdｏｍainsettings"中快速设定非线性函数,给定初始推测值，然后用非线性求解器得到收敛解.但是如果函数中没有线性项可以放在公式（4）左侧怎么办？例如，FＥMLAB将taｎ（u）－u2+５=0放在公式(4）左侧时，会生成奇异刚度矩阵。建议在子域中设定一个u的二次派生系数，c＝１。通过与Neｕmann边界条件相结合，这个人为集中因素不能转变基元中解为常数的事实，但是它可以避开刚度矩阵奇异。在一般模式下求根tan(u）—ｕ2+５=０中奇异刚度矩阵的问题可以通过一般模式来避开，它可以求解下式：ﻩﻩﻩﻩﻩ ﻩ(5)其中Г(u，ｕｘ）和（4）式中的系数项功能类似,但是被求解器处理方式不同。在系数模式中，认为系数与ｕ无关,除非使用数值雅可比矩阵，这会产生一些非线性依靠关系-—迭代次数增加。一般模式下构建刚度矩阵时，精确雅可比矩阵会将Г和F对u进行微分。通常一般模式比使用数值雅可比矩阵的系数模式需要的迭代次数少。下面使用的雅可比矩阵，不需要任何特殊处理来避开刚度矩阵奇异.总的来说,一般模式在处理非线性问题时比系数模式更加有用。从我的个人观点看来,系数模式是ＣＯMSＯＬＭuｌtiｐhyｓiｃｓ的一个遗传特性--MAＴLAＢ的偏微分方程工具箱，它在很多方面是ＦＥＭＬAB和ＣOMSＯLＭｕltipｈyｓiｃs的先驱，它广泛使用了系数形式.此外，数值雅克比矩阵系数公式是一个长期存在的标准有限元方法，所以相对于其它代码，它是一种格外有用的公式.表2给出了应用一般模式的步骤——对我们刚才的步骤做了小小改动.表２一般模式下的求根实例。文件名：rｔfingen.ｍpｈMｏdeｌNaｖｉgaｔor1D，CＯMSOLMultiphyｓics:PＤEmｏdｅs,ｇｅｎerａl模式Oｐｔioｎs设定Ａxes/Ｇrid为［０，１]Dｒａw名称：Inｔeｒvａｌ；Ｓtaｒｔ＝0;Stop=1Ｐhyｓicｓ菜单:BouｎdａｒySettｉnｇs设定两个端点（域）为Neumａnn边界条件Phｙsics菜单:SubｄｏｍａｉｎSettinｇs设定Г=0;=0；F=ｕ３＋u2—4*u+２初值u（t０）＝－0。5Meｓh模式设定最大基元尺寸，generaｌ＝1;点击ＲeｍeｓhSolve使用默认设置(ｎonｌineaｒsolver，ｅｘacｔJacobiａｎ）Pｏsｔ－Pｒocessinｇ经过5次迭代，结果收敛。点击图片读出u=0.7３20５1。转变初始条件找到另外两个根虽然这里给出了单变量简洁函数求根的模板（ｒtｆｉｎdgeｎ．ｍｐh)，但是实际上ＭAＴＬAＢ有更简洁的计算子程序,通过内置函数ｆzeｒo和函数声明来求根，现在ＣＯMＳOLMｕｌtiphysics图形用户界面可以使用帮助方程中的帮助变量来求解几何约束,该变量称为状态变量。作为对一般模式求根模型的补充，使用状态变量v依据表3的步骤来求解一个非线性方程。表３在ＯDE设置中求根Physics菜单：ODＥｓetｔｉｎｇs名称：ｖ。方程：ｔａnh(v）—v^2＋5，完成Soｌve菜单:Solverｐaｒａmｅtｅrs选择Staｔionａrynｏnlｉneａｒ求解器，求解，完成Post－ｐroｃｅssiｎg选择Poiｎteｖａluａｔｉｏn,边界２，表达式：ｖReporｔ窗口值:—2．00８８1９，表达式:v，边界:2下一节我们将非线性求根方法应用到一个常见的化工实例中——闪蒸,它用到了COＭＳOＬMｕltiｐhｙsicｓ图形用户界面的更多新特性。２。2求根：闪蒸实例化学热力学中广泛存在求根法的应用实例，由于化学平衡和质量守衡的约束条件很充分，再加上状态方程,就产生了和问题中未知变量个数相同的约束条件。在本节中,我们以闪蒸为例来介绍简洁的求根方法，构建只有一个自由度的系统，这里以相分率作为未知变量。表4闪蒸单元的初始组分以及平衡时的安排系数K组分XI6５℃，3.4bar下的Kｉ乙烷0。00７91６。2丙烷0．１２８15.２i－丁烷０．０８４92。６n－丁烷0．26901.98i－正戊烷0．05890。91n-正戊烷0．１36１0．72己烷０．３15１0。２８液相碳氢化合物混合物通过闪蒸到65℃、3．4baｒ状态。液相组成和每种组分闪蒸“Ｋ"值见表4.我们盼望知道闪蒸过程中气相和液相产物的成分和液相的比例。表4给出了原始成分。组分i的质量平衡满意以下关系ﻩ ﻩﻩ （6)Ｘi是闪蒸前液体的莫尔分数,xi是闪蒸后液体莫尔分数,yi是蒸发产物莫尔分数,是液相产物与总供应摩尔流率的比值。平衡系数的定义为：Ｋｉ＝ｙi／xi。将该方程消去质量平衡中的xi,得到一个关于ｙi和Xｉ的方程：ﻩﻩﻩﻩ ﻩﻩ （７）由于各yi相加为１，所以我们有的非线性方程：ﻩﻩﻩﻩﻩﻩ(8）其中n是组分数量。该函数可用ｒｏot（s)求解，将计算结果回代就可以得到全部产物的莫尔分数。牛顿-拉夫逊方法需要知道当前导数来决定计算方向,在COMＳOLＭｕltiphｙｓiｃs中将解析计算该值。通过下式很容易得到牛顿-拉夫逊迭代方法:ﻩﻩﻩ ﻩﻩ(９)现在再回到ＣOＭＳOＬMｕlｔiphysics求根计算。作为一个练习，我们用一般ＰＤＥ模式来计算求解。我们可以只载入ｒootfindｅr。mph或ｒtfindgen.mｐｈ并加以修改来计算这个问题，但是熟识ＣＯMSOＬMulｔiphｙsics功能也是格外重要的目的.表５闪蒸实例ModelNaviｇaｔｏr选择１D,COMＳOLMultｉｐｈyｓics:ＰDＥmｏdes:PＤＥ，geneｒal形式Draｗ菜单Specifyobｊｅcts:Linｅ，Ｃoｏrdiｎａｔes弹出菜单,x：01名称：ｉnｔerｖal，完成Ｏptiｏｎｓ菜单：Consｔaｎts输入表4中数据X１,０．0079Ｘ2，0．1281X3,0．０849X４,0.２6９0X5,0．0589X6,0.1３6１X７，0.3151Oｐｔiｏnｓ菜单:ＳcalarExpressｉons为（8)式中右侧各项定义表达式t１-Ｘ1／（１－u＊（1—１/K１））t2－X2/(1－ｕ*(１－1／Ｋ２）)ｔ3-X３／（1－u＊（１—1/K3))ｔ4-X４／（1—u＊（１－１/K４))ｔ５—Ｘ5／（1-ｕ＊（1-1/K5）)t6-Ｘ６/(１—ｕ*（1-1/K6）)t7－X7／(1－u＊(１－1／Ｋ7)）Physiｃｓ菜单：Boｕnｄarｙsettｉｎｇs选择域：１和２(按住Cｔｒｌ键)选择Ｎeumanｎ边界条件采纳默认设置q=０,ｇ=0完成Ｐhｙｓｉcｓ菜单:Ｓｕbdomaｉnseｔtinｇs选择域:１设定F=1＋t１＋ｔ2+…+t7；=0选择Iｎiｔ选项卡，设定u（t0）=０.5Mesh菜单：Mｅsｈparａmｅｔeｒs设定最大基元尺寸为1点击Rｅmesｈ，完成Solｖe菜单：Solverｐaraｍeteｒｓ选择Stationaｒynｏｎｌineａr求解器求解,完成Post—ｐｒoｃeｓsｉｎg：PointEvaｌｕation选择边界：1,表达式:u,完成Reｐoｒt窗口,值：０。45８５0９，表达式：ｕ启动COＭＳOLＭｕlｔｉｐhyｓｉcs,打开模型导航栏.如果你已经打开COMＳOLMultｉｐhyｓics对话框,那么将你的工作空间保存为mph文件或者将命令保存为m文件，然后打开“Fｉle”菜单选择“Ｎeｗ”。依据表5的步骤建立闪蒸模型。注意本例中有两个新增内容--“Optｉons:Constａｎｔs”和“Optｉons：Expresｓｉons”.在这里可以定义常数，然后在ＣOMＳOLＭｕltiｐhyｓｉｃｓ中任何可以使用纯数字的输入框中使用定义的常数。表达式和常数类似，可以在任何ＣOMＳOLＭulｔiｐhｙｓics数字输入框中使用,但是不同之处在于表达式依靠于因变量。定义的常数和表达式同样可以在后处理过程中使用。后处理过程显示t1和t2为：Ｖalue:—0．01３865,Ｅxｐrｅｓｓiｏｎ:ｔ1，Bouｎdaｒy：1Vａlue:—０．２0344１,Eｘpｒeｓsｉoｎ:t2，Ｂouｎｄary：１另外求解器日志中常常包含有用的信息。打开“solｖeｒ"菜单，在最底部选择“ＶｉｅｗＬog”，会弹出求解器日志对话框。它显示了求解器何时开头计算，执行了什么命令，如何从初始状态得到计算结果。这里经过三次迭代，肯定误差达到１0—9。如果你的解不是不收敛或收敛很慢的话,很容易得到这样的信息。练习计算方程的根.由于该函数是三次多项式，所以存在解析的无理数解，。ueｘｐ（u）—1计算方程的根。这是个超越方程，意味着不存在解析的无理数解.如果你使用系数模式,设ｃ=1来帮助收敛。3．方法２:步进法数值积分数值积分是数值分析的支柱.在有计算机之前,科学计算的首要工作是制作特殊方程手册,其中几乎包括了全部特殊常微分方程的解。那么用什么方法计算得到的？就是用一维数值积分。一维数值积分分为两种：初值问题(IVP)和边界值问题（ＢVＰ）。我们下一节再介绍边界值问题.最容易积分的是初值问题，由于全部的初始条件都在一点定义,可以依据一阶导数直接步进积分。显然，如果常微分方程是一阶的，例如： ﻩ ﻩﻩ（10)那么当△ｔ→0时,（1０)中其次项严格成立。它满意欧拉法的条件,是积分一阶常微分方程最简洁的方法.对于一维情况，可以简洁地依据ｆ在（xn,yｎ)点的导数向前步进积分，这里n是指积分过程中的第n次离散化步骤。所以，ﻩﻩﻩ ﻩﻩ(１１）这里假设导数变化不超过步进量ｈ，这只对线性函数才严格成立。对于任何有曲率的函数,该假设都不成立。下面考虑在图1中，如果步进量h取得过大会造成什么错误.显然，欧拉法需要改进的一点就是增大步进量，由于它需要小步进量以保证精准性。欧拉法也被称作“一阶"精准性，由于误差只能降低到一阶h的程度.图1数值积分中的欧拉步进龙格库塔方法所以如果我们盼望增大步进量，那么就需要一个“更高阶的方法"，随着步进量的减小,误差快速降低.k阶方法误差大约在hk量级。假设我们忽视了曲率,可以通过计算斜率f（x）在ｘｎ和xn+１之间几个中间点的值来计算ｙ（ｘ）的曲率。可以首先依据初始导数计算间隔中点的值，然后再用整个间隔中点的微分值来获得二阶精准性。ﻩ ﻩﻩ ﻩ（１2)这样我们通过对两个函数的计算又得到一阶精准性.所以,如果使用一阶方法，Ｎ次计算得到Ｏ（1/N）误差,而二阶方法2Ｎ次计算得到O(1/４Ｎ２)误差，如果使用一阶方法的话,这需要Ｎ2次计算。更高阶龙格库塔方法我们可以做的更好吗?显然，如果用三个中间点可以获得三阶精准性，用四个中间点可以获得四阶精准性等等。不过更高阶的方法需要更多的编程工作，所以也要考虑我们的时间.但是从本质上讲，我们计算函数的k阶导数不肯定光滑。可能在增加“近似精准性”的同时,高阶导数的整体误差会变得很大,如果是这样，每次连续步进都可能导致误差的急剧增大。这表明高阶方法没有低阶方法稳定。常微分方程通常选用四阶龙格库塔方法积分，这样程序比较紧凑,精准性高，也有很好的稳定性。其它方法在数值积分编程的时候，还需要注意两个问题:数值稳定性：即使使用高阶精准性的方法，你的积分结果与检验值偏差仍然很大，这可能就是由于你的数值方法不稳定。你可以通过降低步进量来获得数值稳定性，但是这也意味着更长的计算时间。如果你的计算量很大，计算时间太长,可以选用半隐式方法,例如猜测校正法等。刚性系统：物理机理起作用时，刚性系统通常会有两个完全不同的长度或时间尺度。刚性系统可能会消灭上面提到的“系统非稳定性”现象，或者非物理现象的振荡。可以参考Geａr[2］的书来解决这个问题。３.１数值积分简洁实例高阶常微分方程通常通过降阶和步进法求解，考虑如下方程:ﻩﻩﻩﻩ ﻩ （13)除非q(ｘ）和r(x)是常数,那么你很幸运能够从大多数分析方法教材中找到答案。也有一些特殊的q（x)和r(x)可以求得解析解,但是最好能计算出各种情况下的数值解。为什么呢？由于为了搞清楚y（x）,你需要画出曲线,得到解析解时你也需要花费肯定的精力在作图上.那么应该怎么做?首先将上面的二阶系统降低为两个一阶系统：上式的常微分方程能够仿照（11)和（１2）进行时间步进法数值积分。举个简洁例子：ﻩﻩ ﻩﻩﻩﻩ（１4）降阶产生了两个一阶常微分方程：ﻩﻩﻩ ﻩ ﻩ(15）假设初始条件为u１=1，u2=０，可以建立零维空间系统来积分这一对常微分方程。打开ＣOMＳOＬMｕltiphysｉcs，等待模型导航栏窗口。如果已经打开一个CＯMＳOＬMulｔiphｙｓｉｃｓ窗口，把工作空间保存为mｐh文件或者把命令保存为ｍ文件，然后在“Fｉｌe"菜单中选择“Neｗ”选项。依据表6的步骤进行设置。表6时间步进积分的简洁实例ＭｏｄｅlＮaｖigａtor选择１D,COＭSOＬMｕltiphｙｓicｓ：PＤEｍoｄeｓ:ＰDE，genｅｒal模式，插入因变量：u1u２Drａw菜单Speｃｉｆyｏｂｊｅcts:Lｉｎe。Coordinaｔes弹出菜单,x:01名称:intｅrvａl,完成Physｉcｓ菜单：Bｏuｎｄaｒysetｔinｇs选择域：１和２（按住Ｃｔrl键)选择Ｎｅｕmann边界条件保持默认选项ｑ＝０，ｇ=0完成Phｙｓiｃs菜单：Subdomaｉnsｅttiｎｇs选择域：１设定Ｆ１=—ｕ２,F2=u1选择Iｎiｔ选项卡，设定u1（ｔw)=1Mｅsｈ菜单:Meｓｈpａrameterｓ设置最大基元尺寸为1点击Ｒemesh，完成Ｓoｌｖe菜单：Sｏｌｖｅrpａｒameters选择ｔｉｍｅdｅpeｎdent求解器,在Geneｒａl选项卡中输入Ｔimｅｓ：ｌinspaｃe(0,２*ｐｉ，5０)求解，完成Post－pｒoｃessing:Crosｓ-secｔioｎｐｌoｔparａｍeｔersＰｏinｔ选项卡，接受ｕ１的默认设置Geｎerａl选项卡，接受全部时间的默认设置完成这个例子给了我们两个独立变量u１、u２和空间坐标x。注意“Sｕbdoｍaiｎsｅｔtｉnｇs”窗口左上角方程中的矢量标志.由于是矢量变量,全部输入选项卡都是矢量(F）或者矩阵（dａ)输入。MＡTＬＡB命令linspace（0，2*pｉ，５0）生成50个基元的矢量，给出从0到２π的数据,u1(t=2π）=1.0０441４.假设解析解是ｕ１（t=2π)＝1,结果不够精准（0．４％误差)。之前FＥMＬAB允许用户在MＡＴＬAＢ和FEＭＬAＢ中选择针对时间积分的预置的求解器。也允许用户在求解器变量对话框的通用选项卡中调整容错值(相对和肯定)。将相对误差调整为0。0０１，肯定误差调整为0.00０1,得到一个相对较好的计算结果u１（t＝2π）＝0.998０27.注意累积的全局误差与求解方法精准性０。0０1的量级全都.图２一个周期的u１(ｔ)图3一个周期的ｕ２(t)图2和图3表明如果设定足够小的步进时间，FEＭＬAB能够完成高精度的正弦和余弦函数数值积分。虽然我们一般认为正弦和余弦是“解析函数”，但是通过比较，还是能清楚地看到解析函数和需要数值积分的函数是似是而非的——没有比Ｂｅｓsel函数和椭圆方程更有解析性的函数。练习使用物理场菜单中的常微分方程设置和相同的初时条件,试着积分方程（１5）,假设状态变量为v1和v２。基于时间的常微分方程和代数方程的唯一区分就是必须使用符号代替ｄv１/dt和v1ｔ等。４．管式反应器数值积分本节介绍在化工管式反应器设计时，如何同时求解一阶非线性常微分方程组。通常管式反应器设计的关键要素是反应器的长度。一个气态酒精脱水管式反应器，工作在2bar和１５0℃下。反应方程式为:已有试验表明反应速率可以表示如下，其中CA是酒精浓度(mｏl/L)，R是酒精的消耗速率（ｍoｌ/s／m3)：反应器直径为0.0５ｍ，入口酒精流速为１0ｇ／s。问：如何确定反应器的长度，进而获得各种酒精转化率？我们盼望确定酒精出口摩尔分数分别为0．５,０．4，0.３,0．２和0．1时的反应器长度。化工设计理论假设反应热很小，反应器内为柱状流，抱负气体,可以确定反应流体由四个常微分方程掌握，其中独立变量为CＡ，CＷ（水浓度),Ｖ（速度）和ｘ（反应器长度）: ﻩﻩ ﻩﻩ (16)最后一个方程说明表面速度使得反应器长度和反应物在反应器内的存在时间之间存在一个平衡.初始入口条件为：ﻩﻩﻩﻩ ﻩ (１7)方法从初时条件和化学计量系数可以看出,CW＝CE(酒精浓度），由于温度和压力设为常数，所以浓度C也为常数,可以看作是抱负气体，有：ﻩ ﻩﻩﻩﻩﻩ（１8）初始进口速度Ｖ0可以由给定流速、入口密度(酒精分子量为46kg/kmol）和管的横截面积计算出。该方程需要对时间t数值积分到需要的酒精摩尔分数。使用较为简洁的欧拉方法或者龙格库塔方法积分。读者也许注意到,积分ＣA可能不需要其它变量，但是CＷ，V和x与时间ｔ严格相关。但是由于我们需要求出对应摩尔分数的x值，所以最好同时求解四个变量。部分结果计算得到的数值解为：ﻩ ﻩﻩ ﻩﻩﻩ(19）其中CA/Ｃ如图4所示。图4酒精浓度随时间t变化曲线ＣOMSOLＭultｉphysics计算我们盼望再次建立虚拟的零维模拟空间，这次使用四个独立变量。打开CＯMSOＬMultiphyｓｉｃs，等待模型导航栏窗口。依据表7的步骤建立管式反应器模型。该模型给出四个独立变量u1、u２、ｕ3、u４和一个空间坐标x。由于变量比较多,建议用变量名称进行常数设定.进一步来说,由于使用了速度定律，用矢量表达式更为便利.表７COMＳＯLMｕlｔiphyｓicｓ中管式反应器设计建模步骤ＭoｄelNａｖｉｇatoｒ选择１D空间维数,COＭSＯLMｕltｉphysｉcｓ：PDEmodes：PDＥ，ｇeneral形式设定独立变量:u1u2u3ｕ4选择Elｅmｅnt:Lagrangｅ—Lineａr，完成Ｄｒａw菜单Ｓpeciｆyoｂjects：Lｉne，Cｏordinaｔes弹出菜单，x：01名称：intｅrval，完成Opｔiｏｎｓ菜单:Ｃonstanｔs如下填写变量表格：名称EｘprｅssioｎP20０00０T423R8３14ＭＭ46Fｌowrate0。01Dia0。0５CP/（RT)arｅapi＊Dｉａ^2／4rhoMＭ＊CvｅｌFlowraｔe/rho/ａera完成Options菜单：ScａlarExｐｒessions定义ratｅ=5２.7＊ｕ1^2／(１+０.０１３/u１）Physiｃs菜单：Ｂouｎdaｒysettinｇs选择域:1和2（按住Cｔrl键）选择Neｕmａnｎ边界条件默认设置ｑ=0g＝0,完成Phyｓｉcs菜单：Sｕｂdomainsｅtｔinｇs选择域１F选项卡；设定F1=-ｒate*（1+u１/C);F2=rａｔe＊(1-u２/C)；Ｆ3=ｒaｔe*u3/Ｃ；F4=u3/ＣＩnit选项卡；设定ｕ1（ｔ0）＝Ｃ;ｕ３(t0)=veｌ，完成Ｍｅsh菜单:Ｍｅsｈpａramｅters设定最大基元尺寸为1点击Ｒemｅsh，完成Ｓolve菜单:Solvｅrｐarameterｓ选择timedependenｔ求解器，在Gｅneｒaｌ选项卡中输入Tiｍｅｓ:linspａｃｅ（0，10，１00）求解，完成Post—ｐroｃessingCｒosｓ－ｓecｔionｐｌotpaｒａｍｅtｅrsＰｏint选项卡,接受u１的默认设置Ｇeneral选项卡，接受全部时间的默认设置完成试着在整个时间范围内画出u１,ｕ2,u3和u4的曲线，与图4是否保持全都？与之前计算结果是否相符？练习计算以下方程组中（ｘ=1)的值，并画出随x变化曲线，ｘ取0到３区间。连续搅拌反应器中的一阶可逆反应系统可以用线性常微分方程组表示。例如,考虑如下异构反应:正向反应速率为k1和k3,相应逆向反应速率为k２和ｋ4。一阶动力学常微分方程组为:ﻩﻩﻩﻩﻩ（２0）你也许会感到奇怪，由于该方程组为线性,它有通用的解析解。虽然是通用解析解,但是对动力学系统的反映并不多。假设刚开头时，CA＝１,ｋ1＝1Hz,k2=0Ｈz，ｋ3＝2Hz,ｋ４=３Ｈｚ。针对该初值问题，画出浓度随时间变化的曲线。5。方法3：常微分方程数值积分前一节介绍了步进法数值积分，通常也称作“时间步进",但是在反应器设计中，多是空间积分问题。步进法是挨次求解未知项，而其它常见积分方法是同时求解某一网格点的全部未知独立变量。步进法只能求解初值问题（ＩＶP),初值个数必须严格等于常微分方程组的阶数。但是对于二阶或者更高阶的系统,可能会用到其次类边界条件-—边界值问题，在一维情况下，这些边界条件只消灭在域的起点和终点，这就是两点边界值问题。步进法也能够牵强求解边界值问题——人为设定初值，通过尝试和误差修正找到满意边界值问题的初始条件。对于更高维数的偏微分方程，边界值问题发生在整个域的每一个边界上。有限元方法的一个最主要优点就是能够很容易地求解两点边界值问题。例如，一维反应和集中方程：ﻩﻩﻩﻩﻩ ﻩﻩ(2１）这里u是组分浓度，D是集中系数,L是域的长度,R(u）是反应消耗速率,ｘ是无量纲空间坐标。如果未知函数u(ｘ）在x=xj=j△x网格点能够用不连续值uj=u（xj）近似,那么依据中心差分，方程变为：ﻩ ﻩﻩﻩ（22）其中Ｒj=R（uj）,Mij是对角元素为－２，上、下对角元素为１的三对角矩阵：ﻩ ﻩﻩﻩﻩﻩ（23）。通过矩阵转置和对ｕiｎ进行积分可以求解该方程,这里n指第n次推测:ﻩﻩﻩﻩﻩﻩ (24）Rj=Ｒ（uiｎ－1)。无论是初值问题还是边界值问题，都可以将(23)式中矩阵Ｍ的行向量变得符合边界条件.(2３）式假设在ｘ=0和ｘ=１处都有ｕ＝0。这是狄利克雷边界条件，也是有限微分法的自然边界条件-—说它自然是由于如果在设定边界条件时没有转变（23)式中的行向量,那么就会消灭这样的边界条件。下面介绍如何用CＯMSOLMｕltipｈysics在一维域上求解方程（21),对于一阶反应R（u)=ku,典型的无量纲变量Dａｍkｏhler数为:ﻩﻩﻩﻩﻩ(２5)边界条件为ｘ=０时u＝1,在u=1处没有流量。下面结合MATLＡＢ来做这个练习,进而了解ＣＯMＳOLMulｔｉphyｓiｃs如何表示计算数据和模型输出的结构.在wｉndows中,结合MＡTLAB的COMSOＬＭｕltiphyｓics是一个可选的桌面图标，在UＮIX/ｌinux中,该功能可以通过在命令栏中执行以下语句实现：Coｍsｏlmaｔlａbpath－ｍlnodesktop－ｍlnosplaｓｈ其中“mａtlａb”告知femlab打开一个ｍaｔlaｂ命令窗口,“ｐａtｈ"依据COMＳOＬ命令库建立maｔlab命令窗口。首先运行CＯＭSOLMｕｌｔｉｐhｙｓｉｃs和模型导航栏窗口，然后依据表8中的步骤建立模型.该实例给了一个独立变量ｕ和一个一维空间坐标x。ｈ和ｒ是系数模式中狄利克雷边界条件的两个句柄。如果你盼望边界速度u为常数Ｕ0,可以通过设定ｈ＝1,ｒ=U０来实现。点击三角形按钮生成默认网格（１５个），然后点击两个嵌套三角形按钮加密网格（30个).表８反应－集中系统中两点边界值问题的常微分方程实例MｏdeｌNaｖigaｔor选择１D空间维数，CＯＭSOLMuｌtｉphysics：PDEｍｏdｅs:PDＥ,coｅｆｆicieｎt形式设置因变量:u选择Ｅlement：Ｌaｇraｎge-Linｅａr，完成Draw菜单Speｃiｆyobjects:Ｌinｅ，Cooｒｄinatｅs弹出菜单,ｘ：０1名称:ｉnteｒvaｌ，完成Phyｓiｃs菜单：Ｂｏｕndaryｓｅttｉngs选择域1选择Diｒicｈｌet边界条件，设定ｈ=1;ｒ=1选择域2选择Ｎeuｍann边界条件,完成Phｙｓｉcs菜单：Suｂdomａiｎsettiｎｇｓ选择域1设定Ｃ＝—１;f=0。83３*ｕ;=0选择Ｉｎｉt选项卡;设定u（ｔ０）＝1-x完成Ｍeshinｇ点击三角形按钮进行网格划分Solve菜单:ＳｏlveｒｐａrametｅrsStationarylineａｒ求解器默认设置,完成Ｇenｅｒal选项卡，设定解的形式为“Ｃｏｅfｆicienｔ"求解(＝）Pｏｓt－processｉng：Pｏｉntｅｖalution选择边界２,输入表达式ｕ。Rｅｐorｔs窗口显示:值：０.86116７,表达式：ｕ，边界:2依据计算结果可以画出图５，显然边界条件需要满意：x＝0处u=１，在ｘ=１处斜率减为零。但是COＭSＯLMuｌtipｈｙsics有没有依据我们设想的求解呢？图5稳态下的浓度曲线现在用稳态非线性求解器再计算一次。首先查看求解器菜单中的日志，迭代13次收敛.你注意到最终值从0。8６降到0.６9了吗?为什么会这样呢？线性（系数型）求解器只在初始时刻u（t0)=1—x时计算一次R（u）,所以方程（2４）只需要一次迭代。而非线性求解器在每次迭代时计算一次Ｒ(u）,计算时使用上次迭代的u值.所以非线性求解器在向收敛靠近时，经过几次迭代可能完全“忘记”了初时推测。检验一下这种解释对不对。转变初值将会转变稳定的线性解.返回子域设置对话框，设定初值为ｕ（t0)＝1,最终得到的u（x=1）是多少？再试一下稳态非线性求解器，能得到和其它初值一样的结果吗？该例说明为方程选择正确求解器的重要性.如果函数f依靠于任何非独立变量，就应该选择非线性求解器。线性求解器更快,但是它假设偏微分方程的系数不依靠于非独立变量ｕ(否则为非线性问题）。如果不确定，还是应该用非线性求解器。别忘了,(21）式满意R（ｕ）=ｋu，是线性问题，但是COMSＯLＭultiphysicｓ使用非线性求解器才得到正确解!收敛比较慢也是模型形式导致的结果——选中雅克比求解器选项，非线性求解器只需要两次迭代就可以得出正确结果.我们一开头认为（2２)式是该问题有限微分矩阵方程，但是实际上最后用(２４）式描述CＯMSＯLMuｌtiphｙsics的有限元问题。由于我们这里使用拉格朗日线性微元,在这种特殊情况下有限元和有限微分矩阵全都。依据这一点,我们看一下MＡTＬＡＢ/ＣＯＭSOLSｃｒｉpｔ如何计算该ＣOＭＳOLＭｕｌtｉphysics问题.打开“Fｉle”菜单,选择“EｘporｔFEMStｒuctuｒｅasｆem”。这就把现在的计算结果转化为MATＬAB/ＣOＭSOLScｒiｐｔ工作空间的数据结构,然后用MAＴLAＢ/ＣOＭSOLScript预置的函数和命令来计算它。在MAＴLＡB/CＯMＳOLScｒｉpｔ工作空间输入以下命令：〉〉x=feｍ.mesｈ。p；〉＞u＝fｅm．sｏl.u;＞〉plｏt(x，ｕ)将会弹出一个包含网格节点u值的MAＴLＡB／CＯMＳOLSｃript图形窗口。不要怀疑你的图形失真。这是由于CＯＭSＯLMultiphyｓics存储网格点和对应结果的方式是为了使矩阵变得稀疏和紧凑。可以用MATLAB排序命令对网格点进行排序，进而使图形有意义。在ＭAＴLＡＢ工作空间输入以下命令：〉>［ｘx，idx］＝sort（ｘ）；>＞ploｔ（xｘ,ｕ(iｄx)）最后考虑MAＴLＡB对有限元矩阵的访问。Fem文件结构没有包含有限元刚度矩阵,但是它包含重建这个矩阵的ＦEMLAＢ函数的必要信息。这对有限元方法格外重要,这个FＥMLAＢ函数是ａsseｍble,输入以下命令：>＞［Ｋ，L,Ｍ，N］=assembｌｅ（fｅm）；>＞K’／15这张图包含上次COMSＯＬMｕltipｈysics的计算结果，类似图5。实际上,我们之所以能够快速搞清楚ｆem求解结果的结构，是由于这是只有一个独立变量的一维问题。多变量和多维产生的网格和结果的数据结构只能用COＭSOＬMuｌtiphｙｓiｃｓ的工具／函数来读取。图６给出了由MAＴLAＢ命令spｙ(K’)产生的松散结构。很有必要将它与（２3)式相互比较,由于能够清楚地看到具有统一网格划分的一维拉格朗日线性基元和生成矩阵方程的有限微分法是类似的。图6由ＭATLAB命令sｐy（)产生的松散结构下面看一下MＡＴＬAB中矩阵的松散结构,全部的元素只有1，-2和1，位于三条对角线上。这是有限元法的刚度矩阵，对于未知类型，它等于(23)式.如果返回“Ｓuｂｄｏmａｉnsｅtｔings”对话框的“ｅlｅment”选项卡,选择拉格朗日二次元素，再次求解并导出FEＭ，生成上面的K，你会发现虽然矩阵也很疏松,但是和拉格朗日线性元素完全不同。练习偏微分方程的系数形式有一项αｕ，令α=0．8３3，f=0，再次计算以上“反应—集中"的例子，比较稳态线性和非线性求解器得出的结果。你能解释为什么这样的表达式会导致这样的结果吗？这样的表达式对刚度矩阵Ｋ又有什么影响?如果把Da选为某特定值,例如Da△x２=2，使得对角元素几乎为零，你认为求解时会消灭什么困难？6。方法４:线性系统分析MＡTLAＢ和COＭSＯＬMｕｌｔiphysics的核心是线性系统分析.本节将简要回顾线性算法理论—-本科生工程数学中的“矩阵方程”就是这类典型问题。好消息是不需要你自己动手编程处理矩阵，这就是MATLAB的用途：供应一个工程矩阵计算子程序库的用户界面。应当注意到，对于本章的例子,ＣOMSＯLSｃｒiｐt也能很好地实现这个功能.以前科学计算主要集中在矩阵计算的高效和松散方法上。Golub和ＶａnLoａn[3］的书是一本很好的专家级矩阵计算指南,但是对于ＭAＴLAＢ入门水平，Hａnsｅｌmaｎ和Littlｅｆieｌd[4]的书是个不错的选择。简洁来讲,标准的矩阵方程如下：ﻩﻩﻩﻩ （26)这里有N个与Ｍ方程有关的未知变量ｘj。系数ａiｊ和右侧常数项ｂi都是未知数。在工程应用中，通常将模型转化为这样的线性方程系统.例如，质量和能量守恒通常都会生成这样的线性方程系统。求解可能性当N=M时,约束条件和未知数相同，所以可能求出唯一解ｘj。但是如果一个或者多个方程彼此线性相关(行退化）,或者全部的方程都只含由某些特定变量组成的未知数（列退化）,就不能得出唯一解。对于稀疏矩阵，行退化和列退化相同，退化方程会导致奇异。但是从数值求解讲，至少还有两点会导致错误：如果方程彼此不是严格线性相关,一些方程也可能在计算过程中由于舍入误差而导致变得格外接近线性相关。求解过程中累积的舍入误差可能会“沉没”真实解，这在N比较大的时候容易发生。这种情况下，计算过程没有出错，但是计算结果会不满意原始方程。线性系统指南没有“典型”的线性方程系统,但是一个粗略的想法就是舍入误差变得不行忽视：当Ｎ在20—５0之间时,可以用单精度常规方法求解,不需要对以上两种错误进行修正.当N在几百数量级时可以用双精度求解。当N在几千数量级时，如果系数矩阵稀疏（几乎都是零)，由于系数特性可以求解.对于稀疏矩阵,ＭAＴLAＢ有特殊的数据形式和一套针对性的函数.但是在工程实际和物理过程中,有很多问题本身就是奇异或者格外接近奇异的。你可能会发现N=１0也很难求解.用奇异值分解法有时可以解决奇异值问题，将其分解为非奇异值。数值线性代数的常规任务方程（2６)可以写为紧凑矩阵方程形式:ﻩﻩﻩﻩ ﻩ ﻩ（27)求解未知矢量x，这里Ａ是方阵，ｂ是已知矢量。A为常数,求解多个ｂ矢量时的解。计算方阵A的逆矩阵A—１.计算方阵A的行列式。如果Ｍ〈N，或者Ｍ＝Ｎ但是方程退化，这时有效方程数比未知量少-—欠定系统。这种情况下可能没有解或者有多个解。解空间由特殊解xp乘以任何线性相关的Ｎ－M维矢量构成，称为矩阵Ａ的化零空间.我们的任务是通过奇异值分解找出解空间。如果M〉N，通常没有符合(26）的矢量解x。常常会消灭这种过定系统,我们的任务主要是寻找最接近于满意方程的折衷解。通常接近程度以（２６）方程左右两侧差值的“最小二乘法"来表示。MＡＴLＡB矩阵计算通过ｉｎv(矩阵）命令可以很容易获得逆矩阵。矩阵方程的解以矩阵除法运算符“\”表示:>＞A=［3—10；—16—2；0－２10］;〉〉B＝[1；５；２6］；〉〉Ｘ=Ａ\BＸ=１.０0002．00０03。００００行列式行列式多用在稳定理论和评估矩阵奇异性上。为什么需要知道行列式？大多数时候你是想知道行列式是否为零。但是,当行列式为零、或者数值上格外接近于零时，由于前面提到的“舍入误差”缘由很难进行数值计算。这是奇异值分解的另一个实例.MATLＡB用det（A）命令计算奇异值.在MAＴＬAB命令行中手工输入以下矩阵，或者从maｔｒix２。daｔ文件中复制以下矩阵：〉＞A＝［0。4５，—０.24４1１1,－0．０1９３３7３，0.32397２,-0。１１８829;－0．2４４11１，０．６８403６,－０．103427,０。2055６9,0。0029２382;-０.0１9３３73,—0。103427,0．8２９5，０.0１８9674，—0。０１1169；0.3２３972，０。20５56９，０.01８9674,0.６５9479,0。1９７３88；-０。１18829，０．０02923８2，-０.011１69，0.19738８,０．7769８5]用det命令可以计算出行列式:＞>dｅt(Ａ）ａns=—1。9682e—０08主轴理论：特征值和特征向量MATLＡB有预置的计算矩阵特征值和特征向量的函数:＞〉eiｇ(Ａ)aｎs＝-0.００00０。70０00。8000０.9０001.０00０当使用以下形式调用eig（)函数时，它会以矩阵V的列向量形式返回特征向量:>＞［V，Ｄ]=eiｇ（A)Ｖ＝—0．6836０。０00０-０。54６9－０。４7８5—0。０６８４—０.418１0。61６20.18310。4５3０—0.4５４７－０．0８３70.400３０。６189-0。6２320。2４７90.５４０90．2582－0。2４15－0.4042-０.647４—０.２416—0。62720．47５5—0.１1９０—０。5５50D＝－0。00０0０000０0.７0０0００0０００。８0000000０0．9000０0０0０1。0０00ｅｉｇs()函数是eig（）函数的变形，它能够计算稀疏矩阵的特征值和特征向量。在下一节中将结合COMSOLMultｉphyｓｉcｓ介绍它的应用。如果矩阵A的行列式格外接近于零，那么一个特征值也格外接近于零。特征向量是矩阵A的化零空间——给出映射到零的方向：＞〉Ａ＊Ｖ（:,1）ans＝1．0ｅ-0０7＊0。2669０．16３30.０32７－0．２1120.０９43全部特征向量可以通过乘以特征值等于矩阵Ａ的特性来证明这一点，例如:＞〉Ａ*Ｖ（：，2）.／V（:，2）aｎｓ=0．７００00。700０0．7０000。７0０00.7000在ＭＡＴLAＢ中,“./"符号表示元素对元素的除法，上面的冒号表示矩阵全部的列。由于系统接近奇异,全部包含该矩阵的解都很差，对这一点我们不应该感到奇怪。例如:>>B=［0；１；0；1;０];>〉A\Banｓ=1．０e+０0６＊２．14８7１.３14２０.２6３1-1。70０１0.７593由于A的全部元素都是个位数，被除数B矢量元素也是个位数，方程（２7)的解也应该由个位数组成,而不是百万量级。对于化学工程师，这意味着对于一个质量守恒系统，入口流速为1ｋg/ｈr，由于质量守恒限制出口也为1ｋg／hr，但是反应器中流速为100０00０kg／hｒ。这是不行能的。但这就是由近奇异矩阵给出的计算结果。奇异值分解（SVD)奇异值分解在很多情况下可以给出一个较好的解。类似于特征值和特征向量的主轴理论,全部矩阵都只有唯一的分解：ﻩ ﻩﻩﻩ ﻩ(28）这里的U和V都是正交方阵,diaｇ是包含奇异值的对角矩阵.在解出Ｕ，V和ｄｉag的基础上，(27）式可以轻易解得ﻩ ﻩ ﻩﻩﻩ(29)U和V为正交矩阵意味着其转置矩阵也是其逆矩阵.对角阵的逆矩阵为对角元素取倒数.所以求解的唯一问题就是何时一个或多个奇异值(ｄiａgj)接近于零.这时（1／diagj)是一个格外大的数,导致求解结果错误.一个比较好的近似方法是将这些奇异值的(1／ｄｉagj）设为零!ﻩﻩﻩﻩ ﻩ （３０）几乎为零元素的替代向量应该在数量级上最小,近似满意方程。在我们矩阵A的例子中,如果只有一个输出,ＭATLAB命令ｓvｄ()给格外异值，如果有三个输出,命令sｖd(）给出:＞>［Ｕ,D,Ｖ]=svｄ（A）Ｕ＝－０。06８４-0。4785０。５46９０．０000-0。6836-0。4５４70．４５3０－０。1831－0．６162—０．４1８10.2４79－0.6232-０.6189—0．４0０３-0.0837－0。6４7４-0．4０４２0。2415—0．25820．540９－0.5550—0．１１90—０.４75５0。６2７2－0。241６D=1.00０0000000．９00000０0０0。８００000０00０.700０000０－０．00001。００0０V=－0。０６84-0．478５０．546９0.00０0－0。68３6—0．454７0．45３0—0.1831－0.6１6２-0。4１８10.２4７９-0．6２32—0.61８9－０.４003-0。0837—０.6474—0.4０４２０。2４15—0.25820.５409－０．5550—0。119０－0.47550。6２７２-０.2416由于Ａ为对称矩阵,必定有Ｕ=V。最小数量级的奇异值分解命令如下：＞>ss=[1．1．／０。９1./0.８１./0．７０］;＞〉diｎｖ=diag（ss)；〉>V＊ｄinv*U’＊Ｂanｓ=０。０８９３１．2８２00。1４79１。０３1７-0.213０这是一个物理上更容易接受的答案,例如，前面商量的质量守恒系统中的内部气体流速。由于有限元方法基于矩阵,所以线性系统理论对于COMSOＬＭultｉphｙsｉcs建模格外重要。当产生的刚度矩阵变得接近奇异时,COMSOLＭｕｌtiｐhｙsｉcｓ可能不能给出满意要求的结果.MATＬＡB中的矩阵计算和稀疏化方法可轻易地用于检查COMSOLMｕltiphｙsics结果的合理性，也可以用于对系统自然动力学的特征分析。这些想法将在下一节中作为ＣOMSOＬＭｕｌtiphｙsｉcs的一个模型例子加以介绍。6.1非均匀介质中的热传导常常遇到的稳态热传导方程，是可以分析求解的最简洁的偏微分方程:泊松方程。但是对于简洁几何情况,会比较难处理。作者认为某些通解很难推导和收敛［5]。所以，数值解比通解更好,进一步讲，传热过程的任何转变都可能会破坏分析结果.本节我们考虑非均匀厚板中的一维传热问题,在板中有分布热源：ﻩ ﻩﻩﻩﻩﻩ（3１)打开耦合MATLAB/COＭＳOＬＳｃｒipt的COＭSOＬMuｌｔipｈysiｃs模型导航栏：依据表9中的步骤建立有分布热源的介质传热问题。依据计算结果可以画出一条几乎线性的曲线,斜率为－１。由于这是个线性问题,求解器日志显示两步迭代收敛。计算结果得T｜x＝0。5=０.4７4097,而解析解为Ｔ｜x=0.5～0。475：ﻩﻩﻩﻩﻩﻩ （３2）下面尝试ｋ=１－ｘ／２。这种情况下也有解析解,但是由于数值简洁性,需要采纳对数和分支切割。解析解为Ｔ|x＝０。5~0。550，你的计算结果呢?表9非均匀介质中的热传导——分布热源情况ModelNaviｇatｏr选择1D空间维数，CＯMSOLMuｌｔiphｙsics:ＨｅaｔTransfｅｒ:Conｄｕctioｎ:Steadyｓtａteanａlysiｓ,设定因变量：ｕ选择Elｅｍeｎt:Ｌagｒanｇｅ—Linｅar，完成Ｄraw菜单Ｓpｅｃiｆｙobjｅcts：Linｅ,Ｃｏoｒdinaｔes弹出菜单,x:0１名称：inｔeｒｖaｌ，完成Pｈysiｃs菜单：Boundarysettingｓ选择边界1设定边界条件为:温度；T0＝1选择边界2设定边界：T1＝0,完成Phｙｓicｓ菜单：Sｕｂdomａｉnｓｅｔtings选择域1设定ｋ=1;Ｑ=－ｘ＊(１—x）选择Iｎit选项卡；设定Ｔ(ｔ0）=１-x,完成Ｍesｈinｇ点击三角形按钮,划分网格Ｓolver点击求解按钮（＝)Post－ｐrocｅsｓｉngDａtａdisｐlay设定x=０。５值:0．４74０97［K],表达式:Ｔ，位置：(０。5）下面来看线性系统理论。打开“Ｆile"菜单，选择“ＥxｐoｒtFEＭstｒuctureaｓfem”。这是我们其次次导出为FEＭ结构，有必要简洁介绍一下。对于任何MATLAB/ＣOMSOLＳcript变量,在命令行中输入变量名字,ＭATLＡB/CＯMSOＬScrｉpt将会给出变量值或其数据结构。试一下：>〉ｆｅmfeｍ=version：［１x1ｓｔruｃt]ａpｐl:｛［１x１ｓｔrucｔ］｝geom：［1x1gｅoｍ1]meｓh:［１x１fｅｍmesｈ]bｏｒdｅr：１oｕｔｆorm：’general’form：’generaｌ’unitｓ:’SI’equ：[1x1struct]bnd:［1ｘ1strucｔ]dｒａｗ：[１x1struct]xmeｓｈ:[1x１cｏｍ.fｅmｌaｂ。xmesｈ.Xmｅsh］soｌ:[1x1ｆｅmsoｌ]当然,你也可以进一步查看FEM结构,了解整个树的干、枝、叶.我们已经剪除了fｅm.soｌ。u，fem。ｍesh。p1不同的分支包含了完整的COＭSOLMｕlｔiphｙｓicｓ模型-—为fem.gｅom中保存的几何结构定义了方程(fem.appl{1｝．equ）和边界条件（fｅm．appl｛1}．eｑｕ）.与其它变量相同,这些结构可以被ＭAＴLAＢ/CO