COMSOL MULTIPHYSICS和数值分析基础

W. B. J. ZIMMERMAN，B. N. HEWAKANDAMBY

Department of Chemical and Process Engineering, University of Sheffield,

Newcastle Street, Sheffield S1 3JD United Kingdom

E-mail: w.zimmerman@shef.ac.uk

本章主要介绍COMSOL Multiphysics 在零维和一维模型数值分析方面的几个关键内容。这些内容包括求根、步进式数值积分、常微分方程数值积分和线性系统分析。这几乎是所有的化工过程数学分析方法。下面通过COMSOL Multiphysics中的一些常见化工过程应用实例来介绍这些方法，包括：闪蒸、管式反应器设计、扩散反应系统和固体中热传导。

1．简介

本章内容很多，可以分为几个不同的目标。首先介绍了COMSOL Multiphysics的主要工作特性；其次介绍了如何使用这些特性来分析一些简单的，位于零维空间、一维空间或“空间－时间”系统中的化工问题。本章还希望通过展示COMSOL Multiphysics和MATLAB工具在化工过程分析中的强大功能，激发读者对使用COMSOL Multiphysics进行建模与仿真的兴趣。

由于COMSOL Multiphysics不是一个通用的问题求解工具，所以一些目标需要迂回实现。作者在使用FORTRAN、Mathematica和MATLAB解决化工问题方面有着丰富的教学经验，并用这些工具实现过这里所有的例子。而且，扩展化工问题的数值分析也已经在POLYMATH[1]中实现，这似乎只在化工委员会的CACHE项目中使用过。

本书前一版已经介绍过在零维空间中求解非线性代数方程和与时间有关的常微分方程的内容。从概念上讲，零维域就是一个简单的有限元。通过研究某一特定有限元中的变化对理解有限元方法非常有用。但是，COMSOL Multiphysics通过独立对话框设置，使得零维几何方程和与时间相关的常微分方程求解变得非常简单。所以本章将同时采用这两种方法求解这些例子。

2．方法1：求根

典型的数值分析课程会讲解多种求根方法，但是从实际经验来看，只有两种算法非常有用——二分法和牛顿法。我们这里没有列出所有方法，而是重点考虑为什么求根是最有效的数值分析工具。在线性系统中求根非常简单，但是对于非线性系统这就是一个挑战，而所有感兴趣的动力学问题几乎都是非线性系统。对非线性系统的求根起源于对反函数的描述。为什么呢？因为对于大多数非线性函数，“正向”y=f(u)很好表示，但是它的反函数u=f-－1(y)可能不能显式表示、多值（无意义）或根本不存在。如果反函数存在的话，求解反函数其实就是求根的过程——求解满足F(u)=0的u等价于求解F(u)=f(u)-y=0。因为大多数数值分析的目标是在系统约束下计算求解，所以这也等价于对所有的约束取反。COMSOL Multiphysics拥有求解非线性问题的核心函数——femnlin，本节主要介绍用它求解零维非线性问题。

femnlin函数使用牛顿方法求解，由于只有一个变量u，牛顿法通过对一阶倒数F'(u)迭代来求根。该方法首先估计函数的斜率范围，然后再逼近根。该斜

率可以通过理论分析（牛顿－拉夫逊方法）和数值（正割法）方法求得。如果能用任何一种方法求得斜率，就可以用泰勒定律来逼近根。其基本思想就是使用目前猜测值u0的泰勒展开式：

f(u)f(u0)(uu0)f'(u0) (1)

该公式可以化简，忽略(u-u0)的高阶项，计算根如下：

uu0f(u0) f'(u0) (2)

这个方法可以快速地扩展到多维求解空间，例如将u看作未知矢量，“被

f'(u0)除”看作“乘以f的雅克比矩阵的逆”。下一节介绍COMSOL Multiphysics

中的求根过程。

2.1 求根：COMSOL Multiphysics非线性求解器的应用实例

如上节所述，求根本身是一个“零维”活动，至少对于“空间-时间”系统多维未知矢量u来说是这样的。COMSOL多物理场没有零维模式，所以我们临时采用一维模式。这在方面增加了我们不需要的冗余功能。但是由于问题规模较小，COMSOL Multiphysics编码效率高，且现代微处理器的运算速度快，这点就不成为问题了。

启动MATLAB并在命令窗口键入COMSOL Multiphysics。屏幕闪烁几秒后，会出现一个模型导航窗口。按照表1所示步骤，建立一个零维应用模式来解决非线性多项式方程：

u3u24u20

(3)

通过在“Physics”菜单、“Subdomain settings”选项中的设定，使得表1中的每个子域都满足该方程。注意左上角，这是以矢量符号给出的方程。在一维模式下，这个方程可以简化为：

dauuucuauf txxx (4)

显然，α、γ和β在一维简化里是多余的。既然我们是求零维的根，所有公式4左侧的系数都可以设为0。通过重新组合多项式，我们发现a=4和f=u+ u2+2。注意我们是通过设定最大基元尺寸为1、将域离散化为只有一个基元的域，从而

得到零维空间的！

表1 系数模式下的求根实例。文件名称：rootfinder.mph. 选择1D空间维数，COMSOL Multiphysics: PDE modes: Model Navigator PDE, coefficient形式 Specify objects：Line。Coordinates弹出菜单。x：0 1 Draw菜单 名称：interval，完成 Physics菜单： 选择域：1和2（按住Ctrl键同时选中） Boundary settings 选择Neumann边界条件 3采用默认设置q=0 g=0，完成 选择域：1 Physics菜单： Subdomain settings Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing Point Evalution 设定初始猜测值为u(t0)=-2，下面我们来寻找最接近这个值的根。你可能对为什么设置a=4而不是把所有有关参数放进f中感到奇怪，那是因为对公式(4)右侧有限元离散后得不到一个奇异刚度矩阵。

后处理阶段在输出窗口中显示出结果： Value: -2.732051, Expression: u, Boundary: 1.

解析解得出的根为-1-3，数值解能够很好的与其吻合。根据代数中二次公式的解知道，另外一个根是-1+3，第三个根是1。回到子域设置中来，如果最初设置的猜测值为u(t0)=-0.5，则COMSOL Multiphysics解出u=0.762051，又一个很好的近似数值解。如果u(t0)=1.2作为最初猜想时，得到u=1。

该练习体现了非线性求解器和问题的两个特性——(1)非线性问题可以有多个解；(2) 最初猜测值对于求解很关键。牛顿法通常收敛到最接近的解，但是在高阶非线性问题中，这点可能不成立。这些特性在多维求解空间和“空间-时间”相关系统中同样适用。

COMSOL Multiphysics模型mph文件（rootfinder.mph）中包含了所有的MATLAB/COMSOL Script源代码和FEMLAB的扩展功能，以便再次恢复到FEMLAB当前图形用户界面。该文件在网页 http://eyrie.shef.ac.uk/femlab 中可以得到。只需打开“Menu”菜单，选择“Open model m-file”选项，在弹出的“Open file”对话框中选中它，你就可以在“Subdomain settings”中快速设定非线性函数，给定初始猜测值，然后用非线性求解器得到收敛解。但是如果函数中没有线性项可以放在公式(4)左侧怎么办？例如，FEMLAB将tan(u)-u2+5=0放在公式(4)左侧时，会生成奇异刚度矩阵。建议在子域中设定一个u的二次派生系数，c=1。通过与Neumann边界条件相结合，这个人为扩散因素不能改变基元中解为常数的事实，但是它可以避免刚度矩阵奇异。 在一般模式下求根

tan(u)-u2+5=0中奇异刚度矩阵的问题可以通过一般模式来避免，它可以求解下式：

设定：c=0; a=4;f=u^3+u^2+2; da=0 应用，选择Init选项卡：设定u(t0)=-2，完成 设定最大基元尺寸1 点击remesh，完成 求解器选择Stationary nonlinear，求解，完成 选择边界1，表达式：u，完成 2uuea2daF ttx (5)

其中Г(u,ux)和(4)式中的系数项功能类似，但是被求解器处理方式不同。在系数模式中，认为系数与u无关，除非使用数值雅可比矩阵，这会产生一些非线性依赖关系——迭代次数增加。一般模式下构建刚度矩阵时，精确雅可比矩阵会将Г和F对u进行微分。通常一般模式比使用数值雅可比矩阵的系数模式需要的迭代次数少。下面使用的雅可比矩阵，不需要任何特殊处理来避免刚度矩阵奇异。总的来说，一般模式在处理非线性问题时比系数模式更加实用。从我的个人观点看来，系数模式是COMSOL Multiphysics的一个遗传特性——MATLAB的偏微分方程工具箱，它在很多方面是FEMLAB和COMSOL Multiphysics的先驱，它广泛使用了系数形式。此外，数值雅克比矩阵系数公式是一个长期存在的标准有限元方法，所以相对于其它代码，它是一种非常有用的公式。

表2给出了应用一般模式的步骤——对我们刚才的步骤做了小小改动。

表2 一般模式下的求根实例。文件名：rtfingen.mph Model Navigator 1D, COMSOL Multiphysics: PDE modes, general模式 Options 设定Axes/Grid为[0,1] Draw 名称：Interval；Start＝0；Stop＝1 Physics菜单： 设定两个端点（域）为Neumann边界条件 Boundary Settings Physics菜单： Subdomain Settings Mesh模式 Solve Post-Processing 设定Г=0; da=0; F=u3+u2-4*u+2 初值u(t0)=-0.5 设定最大基元尺寸，general＝1；点击Remesh 使用默认设置（nonlinear solver，exact Jacobian） 经过5次迭代，结果收敛。点击图片读出u＝0.732051。改变初始条件找到另外两个根 虽然这里给出了单变量简单函数求根的模板(rtfindgen.mph)，但是实际上MATLAB有更简单的计算子程序，通过内置函数fzero和函数声明来求根，现在COMSOL Multiphysics图形用户界面可以使用辅助方程中的辅助变量来求解几何约束，该变量称为状态变量。作为对一般模式求根模型的补充，使用状态变量v按照表3的步骤来求解一个非线性方程。

表3 在ODE设置中求根 Physics菜单： 名称：v。方程：tanh(v)-v^2+5，完成 ODE settings Solve菜单： 选择Stationary nonlinear求解器，求解，完成 Solver parameters Post-processing 选择Point evaluation，边界2，表达式：v Report窗口 值：－2.008819，表达式：v，边界：2 下一节我们将非线性求根方法应用到一个常见的化工实例中——闪蒸，它用到了COMSOL Multiphysics图形用户界面的更多新特性。

2.2 求根：闪蒸实例

化学热力学中广泛存在求根法的应用实例，由于化学平衡和质量守衡的约束条件很充分，再加上状态方程，就产生了和问题中未知变量个数相同的约束条件。在本节中，我们以闪蒸为例来介绍简单的求根方法，构建只有一个自由度的系统，这里以相分率作为未知变量。

表4 闪蒸单元的初始组分以及平衡时的分配系数K 组分 XI 65℃, 3.4bar下的Ki 乙烷 0.0079 16.2 丙烷 0.1281 5.2 i-丁烷 0.0849 2.6 n-丁烷 0.2690 1.98 i-正戊烷 0.0589 0.91 n-正戊烷 0.1361 0.72 己烷 0.3151 0.28

液相碳氢化合物混合物通过闪蒸到65℃、3.4 bar状态。液相组成和每种组

分闪蒸“K”值见表4。我们希望知道闪蒸过程中气相和液相产物的成分和液相的比例。表4给出了原始成分。

组分i的质量平衡满足以下关系

Xi(1)yixi

(6)

Xi是闪蒸前液体的莫尔分数，xi是闪蒸后液体莫尔分数，yi是蒸发产物莫尔分数，

是液相产物与总供给摩尔流率的比值。平衡系数的定义为：Ki= yi/ xi。将该方程消去质量平衡中的xi，得到一个关于yi和Xi的方程：

yiXi111Ki (7)

由于各yi相加为1，所以我们有的非线性方程：

f1Xi1i111Kin0 (8)

其中n是组分数量。该函数f()可用root(s)求解，将计算结果回代就可以得到所有产物的莫尔分数。牛顿-拉夫逊方法需要知道当前导数f'()来决定计算方向，在COMSOL Multiphysics中将解析计算该值。通过下式很容易得到牛顿-拉夫逊迭代方法：

f'i1n1Xi1Ki111Ki2 (9)

现在再回到COMSOL Multiphysics求根计算。作为一个练习，我们用一般PDE

模式来计算求解。我们可以只载入rootfinder.mph或rtfindgen.mph并加以修改来计算这个问题，但是熟悉COMSOL Multiphysics功能也是非常重要的目的。

表5 闪蒸实例 Model Navigator 选择1D，COMSOL Multiphysics: PDE modes:PDE, general形式 Draw菜单 Specify objects: Line，Coordinates弹出菜单，x：0 1 名称：interval，完成 Options菜单： 输入表4中数据 Constants X1，0.0079 X2，0.1281 X3，0.0849 X4，0.2690 X5，0.0589 X6，0.1361 X7，0.3151 Options菜单： 为(8)式中右侧各项定义表达式 Scalar Expressions t1-X1/(1-u*(1-1/K1)) t2-X2/(1-u*(1-1/K2)) t3-X3/(1-u*(1-1/K3)) t4-X4/(1-u*(1-1/K4)) t5-X5/(1-u*(1-1/K5)) t6-X6/(1-u*(1-1/K6)) t7-X7/(1-u*(1-1/K7)) Physics菜单： 选择域：1和2（按住Ctrl键） Boundary settings 选择Neumann边界条件 采用默认设置q=0, g=0完成 Physics菜单： 选择域：1 Subdomain settings 设定F=1+t1+t2+…+t7; da=0 Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing： Point Evaluation 选择Init选项卡，设定u(t0)=0.5 设定最大基元尺寸为1 点击Remesh，完成 选择Stationary nonlinear求解器求解，完成 选择边界：1，表达式：u，完成 Report窗口，值：0.458509，表达式：u

启动COMSOL Multiphysics，打开模型导航栏。如果你已经打开COMSOL Multiphysics对话框，那么将你的工作空间保存为mph文件或者将命令保存为m文件，然后打开“File”菜单选择“New”。按照表5的步骤建立闪蒸模型。注意本例中有两个新增内容——“Options：Constants”和“Options：Expressions”。在这里可以定义常数，然后在COMSOL Multiphysics中任何可以使用纯数字的输入框中使用定义的常数。表达式和常数类似，可以在任何COMSOL Multiphysics数字输入框中使用，但是不同之处在于表达式依赖于因变量。定义的常数和表达式同样可以在后处理过程中使用。后处理过程显示t1和t2为：

Value: -0.013865, Expression: t1, Boundary:1 Value: -0.203441, Expression: t2, Boundary:1

另外求解器日志中经常包含有用的信息。打开“solver”菜单，在最底部选择“View Log”，会弹出求解器日志对话框。它显示了求解器何时开始计算，执行了什么命令，如何从初始状态得到计算结果。这里经过三次迭代，绝对误差达到10－9。如果你的解不是不收敛或收敛很慢的话，很容易得到这样的信息。 练习

511.1 计算方程f(u)u33u2u022的根。由于该函数是三次多项式，所以存在解析的无理数解，

u1u1

12u11。 2u exp(u)-1

1.2 计算方程f(u)ueu10的根。这是个超越方程，意味着不存在解析的无理数解。如果你使用系数模式，设c=1来帮助收敛。

3. 方法2：步进法数值积分

数值积分是数值分析的支柱。在有计算机之前，科学计算的首要工作是制作特殊方程手册，其中几乎包括了所有特殊常微分方程的解。那么用什么方法计算得到的？就是用一维数值积分。

一维数值积分分为两种：初值问题（IVP）和边界值问题（BVP）。我们下一节再介绍边界值问题。最容易积分的是初值问题，因为所有的初始条件都在一点定义，可以根据一阶导数直接步进积分。显然，如果常微分方程是一阶的，例如：

dyf(t)yttyttf(t) (10) dt那么当△t→0时，(10)中第二项严格成立。它满足欧拉法的条件，是积分一阶常微分方程最简单的办法。对于一维情况，可以简单地根据f在(xn,yn)点的导数向前步进积分，这里n是指积分过程中的第n次离散化步骤。所以，

yn1ynhf(xn,yn)xn1xnh (11)

这里假设导数变化不超过步进量h，这只对线性函数才严格成立。对于任何有曲率的函数，该假设都不成立。下面考虑在图1中，如果步进量h取得过大会造成什么错误。显然，欧拉法需要改进的一点就是增大步进量，因为它需要小步进量以保证准确性。欧拉法也被称作“一阶”准确性，因为误差只能降低到一阶h的程度。

图1 数值积分中的欧拉步进

龙格库塔方法

所以如果我们希望增大步进量，那么就需要一个“更高阶的方法”，随着步进量的减小，误差快速降低。k阶方法误差大约在hk量级。假设我们忽略了曲率，可以通过计算斜率f(x)在xn和xn+1之间几个中间点的值来计算y(x)的曲率。可以首先根据初始导数计算间隔中点的值，然后再用整个间隔中点的微分值来获得二阶准确性。

k1hf(xn,yn)11k2hfxnh,ynk1

22yn1ynk2O(h3) (12)

这样我们通过对两个函数的计算又得到一阶准确性。所以，如果使用一阶方法，N次计算得到O(1/N)误差，而二阶方法2N次计算得到O(1/4N2)误差，如果使用一阶方法的话，这需要N2次计算。 更高阶龙格库塔方法

我们可以做的更好吗？显然，如果用三个中间点可以获得三阶准确性，用四个中间点可以获得四阶准确性等等。不过更高阶的方法需要更多的编程工作，所以也要考虑我们的时间。但是从本质上讲，我们计算函数的k阶导数不一定光滑。可能在增加“近似准确性”的同时，高阶导数的整体误差会变得很大，如果是这样，每次连续步进都可能导致误差的急剧增大。这表明高阶方法没有低阶方法稳定。

常微分方程通常选用四阶龙格库塔方法积分，这样程序比较紧凑，准确性高，也有很好的稳定性。 其它方法

在数值积分编程的时候，还需要注意两个问题： 数值稳定性：即使使用高阶准确性的方法，你的积分结果与检验值偏差仍然很大，这可能就是因为你的数值方法不稳定。你可以通过降低步进量来获得数值稳定性，但是这也意味着更长的计算时间。如果你的计算量很大，计算时间太长，可以选用半隐式方法，例如预测校正法等。

刚性系统：物理机理起作用时，刚性系统通常会有两个完全不同的长度或时

间尺度。刚性系统可能会出现上面提到的“系统非稳定性”现象，或者非物理现象的振荡。可以参考Gear[2]的书来解决这个问题。

3.1 数值积分简单实例

高阶常微分方程通常通过降阶和步进法求解，考虑如下方程：

d2ydyq(x)r(x) 2dxdx (13)

除非q(x)和r(x)是常数，那么你很幸运能够从大多数分析方法教材中找到答案。

也有一些特殊的q(x)和r(x)可以求得解析解，但是最好能计算出各种情况下的数值解。为什么呢？因为为了搞清楚y(x)，你需要画出曲线，得到解析解时你也需要花费一定的精力在作图上。那么应该怎么做？首先将上面的二阶系统降低为两个一阶系统：

dyz(x)dxdzr(x)q(x)z(x)dx

上式的常微分方程能够仿照(11)和(12)进行时间步进法数值积分。举个简单例子：

d2uu0 2dt降阶产生了两个一阶常微分方程：

(14)

du1u2dt

du2u1dt (15)

假设初始条件为u1=1，u2=0，可以建立零维空间系统来积分这一对常微分方程。

打开COMSOL Multiphysics，等待模型导航栏窗口。如果已经打开一个COMSOL Multiphysics窗口，把工作空间保存为mph文件或者把命令保存为m文件，然后在“File”菜单中选择“New”选项。按照表6的步骤进行设置。

表6 时间步进积分的简单实例 选择1D，COMSOL Multiphysics: PDE modes: PDE, generalModel Navigator 模式，插入因变量：u1 u2 Specify objects: Line。Coordinates弹出菜单，x：0 1 Draw菜单 名称：interval，完成 选择域：1和2（按住Ctrl键） Physics菜单： 选择Neumann边界条件 Boundary settings 保持默认选项q=0, g=0完成 选择域：1 Physics菜单： 设定F1＝－u2，F2＝u1 Subdomain settings 选择Init选项卡，设定u1(tw)=1 Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing： Cross-section plot parameters 这个例子给了我们两个独立变量u1、u2和空间坐标x。注意“Subdomain settings”窗口左上角方程中的矢量标志。由于是矢量变量，所有输入选项卡都是矢量(F)或者矩阵(da)输入。

MATLAB命令linspace(0,2*pi,50)生成50个基元的矢量，给出从0到2π的数据，u1(t=2π)=1.004414。假设解析解是u1(t=2π)=1，结果不够准确（0.4％误差）。之前FEMLAB允许用户在MATLAB和FEMLAB中选择针对时间积分的预置的求解器。也允许用户在求解器变量对话框的通用选项卡中调整容错值（相对和绝对）。将相对误差调整为0.001，绝对误差调整为0.0001，得到一个相对较好的计算结果u1(t=2π)=0.998027。注意累积的全局误差与求解方法准确性0.001的量级一致。

设置最大基元尺寸为1 点击Remesh，完成 选择time dependent求解器，在General选项卡中输入 Times: linspace(0,2*pi,50) 求解，完成 Point选项卡，接受u1的默认设置 General选项卡，接受所有时间的默认设置 完成

图2 一个周期的u1(t)

图3 一个周期的u2(t)

图2和图3表明如果设定足够小的步进时间，FEMLAB能够完成高精度的正弦和余弦函数数值积分。虽然我们一般认为正弦和余弦是“解析函数”，但是通过比较，还是能清楚地看到解析函数和需要数值积分的函数是似是而非的——没有比Bessel函数和椭圆方程更有解析性的函数。 练习

1.3 使用物理场菜单中的常微分方程设置和相同的初时条件，试着积分方程(15)，

假设状态变量为v1和v2。基于时间的常微分方程和代数方程的唯一区别就是必须使用符号代替dv1/dt和v1t等。

4．管式反应器数值积分

本节介绍在化工管式反应器设计时，如何同时求解一阶非线性常微分方程组。通常管式反应器设计的关键要素是反应器的长度。

一个气态酒精脱水管式反应器，工作在2bar和150℃下。反应方程式为：

C2H5OHC2H4H2O

已有试验表明反应速率可以表示如下，其中CA是酒精浓度(mol/L)，R是酒精的消耗速率(mol/s/m3)：

252.7CAR

0.0131CA反应器直径为0.05m，入口酒精流速为10g/s。问：如何确定反应器的长度，进而获得各种酒精转化率？我们希望确定酒精出口摩尔分数分别为0.5，0.4，0.3，

0.2和0.1时的反应器长度。

化工设计理论

假设反应热很小，反应器内为柱状流，理想气体，可以确定反应流体由四个常微分方程控制，其中独立变量为CA，CW（水浓度），V（速度）和x（反应器长度）：

dCACR1AdtCdCWCR1WdtC  (16) dVRVdtCdxVdt最后一个方程说明表面速度使得反应器长度和反应物在反应器内的存在时间之间存在一个平衡。初始入口条件为：

CA(0)CV(0)V0CW(0)0x(0)0 (17)

方法

从初时条件和化学计量系数可以看出，CW＝CE（酒精浓度），由于温度和压力设为常数，所以浓度C也为常数，可以看作是理想气体，有：

CpJT8314kmolK (18)

初始进口速度V0可以由给定流速、入口密度（酒精分子量为46kg/kmol）和管的

横截面积计算出。该方程需要对时间t数值积分到需要的酒精摩尔分数。使用较为简单的欧拉方法或者龙格库塔方法积分。

读者也许注意到，积分CA可能不需要其它变量，但是CW，V和x与时间t严格相关。但是由于我们需要求出对应摩尔分数的x值，所以最好同时求解四个变量。

部分结果

计算得到的数值解为：

CA0.1Ct5.65225 x18.5435 (19)

其中CA/C如图4所示。

图4 酒精浓度随时间t变化曲线

COMSOL Multiphysics计算

我们希望再次建立虚拟的零维模拟空间，这次使用四个独立变量。打开COMSOL Multiphysics，等待模型导航栏窗口。根据表7的步骤建立管式反应器模型。该模型给出四个独立变量u1、u2、u3、u4和一个空间坐标x。由于变量比较多，建议用变量名称进行常数设定。进一步来说，由于使用了速度定律，用矢量表达式更为方便。

表7 COMSOL Multiphysics中管式反应器设计建模步骤 Model Navigator 选择1D空间维数，COMSOL Multiphysics: PDE modes: PDE, general形式 设定独立变量：u1 u2 u3 u4 选择Element：Lagrange-Linear，完成 Draw菜单 Specify objects: Line，Coordinates弹出菜单，x：0 1 名称：interval，完成 Options菜单： 如下填写变量表格： Constants 名称 Expression P 200000 T 423 R 8314 MM 46 Flowrate 0.01 Dia 0.05 C P/(RT) area pi*Dia^2/4 rho MM*C vel Flowrate/rho/aera 完成 Options菜单： 定义rate=52.7*u1^2/(1+0.013/u1) Scalar Expressions Physics菜单： 选择域：1和2（按住Ctrl键） Boundary settings Physics菜单： Subdomain settings Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing Cross-section plot parameters 试着在整个时间范围内画出u1，u2，u3和u4的曲线，与图4是否保持一致？与之前计算结果是否相符？

练习 1.4 计算以下方程组中y'(x=1)的值，并画出y'随x变化曲线，x取0到3区间。

y''y'y20y(x0)1 y'(x0)0选择Neumann边界条件 默认设置q＝0 g=0，完成 选择域1 F选项卡；设定F1=-rate*(1+u1/C); F2=rate*(1-u2/C); F3=rate*u3/C; F4=u3/C Init选项卡；设定u1(t0)=C; u3(t0)=vel，完成 设定最大基元尺寸为1 点击Remesh，完成 选择time dependent求解器，在General选项卡中输入 Times: linspace(0,10,100) 求解，完成 Point选项卡，接受u1的默认设置 General选项卡，接受所有时间的默认设置 完成 1.5 连续搅拌反应器中的一阶可逆反应系统可以用线性常微分方程组表示。例

如，考虑如下异构反应：

ABC

正向反应速率为k1和k3，相应逆向反应速率为k2和k4。一阶动力学常微分方程组为：

dcAk1cAk2cBdtdcBk1cAk2cBk3cBk4cC dtdcCk3cBk4cCdt (20)

你也许会感到奇怪，由于该方程组为线性，它有通用的解析解。虽然是通用解析解，但是对动力学系统的反映并不多。假设刚开始时，CA＝1，k1=1Hz，k2=0Hz，k3=2Hz，k4=3Hz。针对该初值问题，画出浓度随时间变化的曲线。

5. 方法3：常微分方程数值积分

前一节介绍了步进法数值积分，通常也称作“时间步进”，但是在反应器设计中，多是空间积分问题。步进法是顺序求解未知项，而其它常见积分方法是同时求解某一网格点的所有未知独立变量。步进法只能求解初值问题（IVP），初

值个数必须严格等于常微分方程组的阶数。但是对于二阶或者更高阶的系统，可能会用到第二类边界条件——边界值问题，在一维情况下，这些边界条件只出现在域的起点和终点，这就是两点边界值问题。步进法也能够勉强求解边界值问题——人为设定初值，通过尝试和误差修正找到满足边界值问题的初始条件。对于更高维数的偏微分方程，边界值问题发生在整个域的每一个边界上。

有限元方法的一个最主要优点就是能够很容易地求解两点边界值问题。例如，一维反应和扩散方程：

D2uR(u) L2x2 (21)

这里u是组分浓度，D是扩散系数，L是域的长度，R(u)是反应消耗速率，x是无量纲空间坐标。如果未知函数u(x)在x=xj=j△x网格点能够用不连续值uj=u(xj)近似，那么根据中心差分，方程变为：

L2xMijujRi Dj1N2 (22)

其中Rj=R(uj)，Mij是对角元素为-2，上、下对角元素为1的三对角矩阵： 21001210 (23) M0100RjR(uj)。通过矩阵转置和对uin进行积分可以求解该方程，这里n指第n次猜

测：

L2xnuiD2Mj1N1ijRj

(24)

Rj=R(uin-1)。无论是初值问题还是边界值问题，都可以将(23)式中矩阵M的行向量变得符合边界条件。(23)式假设在x＝0和x＝1处都有u＝0。这是狄利克雷边界条件，也是有限微分法的自然边界条件——说它自然是因为如果在设定边界条件时没有改变(23)式中的行向量，那么就会出现这样的边界条件。

下面介绍如何用COMSOL Multiphysics在一维域上求解方程(21)，对于一阶反应R(u)=ku，典型的无量纲变量Damkohler数为：

DakLD2103s1103m11.2109m2/s20.833 (25)

边界条件为x＝0时u＝1，在u＝1处没有流量。

下面结合MATLAB来做这个练习，进而了解COMSOL Multiphysics如何表示计算数据和模型输出的结构。在windows中，结合MATLAB的COMSOL Multiphysics是一个可选的桌面图标，在UNIX/linux中，该功能可以通过在命令栏中执行以下语句实现：

Comsol matlab path-ml nodesktop-ml nosplash

其中“matlab”告诉femlab打开一个matlab命令窗口，“path”根据COMSOL命令库建立matlab命令窗口。

首先运行COMSOL Multiphysics和模型导航栏窗口，然后根据表8中的步骤建立模型。该实例给了一个独立变量u和一个一维空间坐标x。h和r是系数模式中狄利克雷边界条件的两个句柄。如果你希望边界速度u为常数U0，可以通过设定h＝1，r＝U0来实现。点击三角形按钮生成默认网格（15个），然后点击两个嵌套三角形按钮加密网格（30个）。

表8 反应－扩散系统中两点边界值问题的常微分方程实例 选择1D空间维数，COMSOL Multiphysics: PDE modes: PDE, coefficient形式 Model Navigator 设置因变量：u 选择Element：Lagrange-Linear，完成 Specify objects: Line，Coordinates弹出菜单，x：0 1 Draw菜单 名称：interval，完成 选择域1 Physics菜单： 选择Dirichlet边界条件，设定h=1; r=1 Boundary settings 选择域2 选择Neumann边界条件，完成 选择域1 Physics菜单： Subdomain settings Meshing Solve菜单： Solver parameters 设定C=-1; f=0.833*u; da =0 选择Init选项卡；设定u(t0)=1-x 完成 点击三角形按钮进行网格划分 Stationary linear求解器默认设置，完成 General选项卡，设定解的形式为“Coefficient” 求解（＝） 选择边界2，输入表达式u。Reports窗口显示： 值：0.861167，表达式：u，边界：2 Post-processing： Point evalution 根据计算结果可以画出图5，显然边界条件需要满足：x＝0处u＝1，在x＝1处斜率减为零。但是COMSOL Multiphysics有没有按照我们设想的求解呢？

图5 稳态下的浓度曲线

现在用稳态非线性求解器再计算一次。首先查看求解器菜单中的日志，迭代13次收敛。你注意到最终值从0.86降到0.69了吗？为什么会这样呢？线性（系数型）求解器只在初始时刻u(t0)=1-x时计算一次R(u)，所以方程(24)只需要一次迭代。而非线性求解器在每次迭代时计算一次R(u)，计算时使用上次迭代的u值。所以非线性求解器在向收敛靠近时，经过几次迭代可能完全“忘记”了初时猜测。

检验一下这种解释对不对。改变初值将会改变稳定的线性解。返回子域设置对话框，设定初值为u(t0)=1，最终得到的u(x=1)是多少？再试一下稳态非线性求解器，能得到和其它初值一样的结果吗？

该例说明为方程选择正确求解器的重要性。如果函数f依赖于任何非独立变量，就应该选择非线性求解器。线性求解器更快，但是它假设偏微分方程的系数不依赖于非独立变量u（否则为非线性问题）。如果不确定，还是应该用非线性求解器。别忘了，(21)式满足R(u)=ku，是线性问题，但是COMSOL Multiphysics使用非线性求解器才得到正确解！收敛比较慢也是模型形式导致的结果——选中雅克比求解器选项，非线性求解器只需要两次迭代就可以得出正确结果。

我们一开始认为(22)式是该问题有限微分矩阵方程，但是实际上最后用(24)式描述COMSOL Multiphysics的有限元问题。因为我们这里使用拉格朗日线性微元，在这种特殊情况下有限元和有限微分矩阵一致。根据这一点，我们看一下MATLAB/COMSOL Script如何计算该COMSOL Multiphysics问题。

打开“File”菜单，选择“Export FEM Structure as fem”。这就把现在的计算结果转化为MATLAB/COMSOL Script工作空间的数据结构，然后用MATLAB/COMSOL Script预置的函数和命令来计算它。

在MATLAB/COMSOL Script工作空间输入以下命令：

>>x=fem.mesh.p; >>u=fem.sol.u; >>plot(x,u)

将会弹出一个包含网格节点u值的MATLAB/COMSOL Script图形窗口。不要怀疑你的图形失真。这是因为COMSOL Multiphysics存储网格点和对应结果的方式是为了使矩阵变得稀疏和紧凑。可以用MATLAB排序命令对网格点进行排序，进而使图形有意义。

在MATLAB工作空间输入以下命令：

>>[xx, idx]=sort(x); >>plot (xx, u(idx))

最后考虑MATLAB对有限元矩阵的访问。Fem文件结构没有包含有限元刚度矩阵，但是它包含重建这个矩阵的FEMLAB函数的必要信息。这对有限元方法非常重要，这个FEMLAB函数是assemble，输入以下命令：

>>[K,L,M,N]= assemble(fem); >>K’/15

这张图包含上次COMSOL Multiphysics的计算结果，类似图5。实际上，我们之所以能够快速搞清楚fem求解结果的结构，是因为这是只有一个独立变量的一维问题。多变量和多维产生的网格和结果的数据结构只能用COMSOL Multiphysics的工具/函数来读取。图6给出了由MATLAB命令spy(K’)产生的松散结构。很有必要将它与(23)式相互比较，因为能够清楚地看到具有统一网格划分的一维拉格朗日线性基元和生成矩阵方程的有限微分法是类似的。

图6 由MATLAB命令spy(K')产生的K'松散结构

下面看一下MATLAB中矩阵的松散结构，所有的元素只有1，-2和1，位于三条对角线上。这是有限元法的刚度矩阵，对于未知类型，它等于(23)式。如果返回“Subdomain settings”对话框的“element”选项卡，选择拉格朗日二次元素，再次求解并导出FEM，生成上面的K，你会发现虽然矩阵也很疏松，但是和拉格朗日线性元素完全不同。

练习

1.6 偏微分方程的系数形式有一项αu，令α＝0.833，f＝0，再次计算以上“反应

－扩散”的例子，比较稳态线性和非线性求解器得出的结果。你能解释为什么这样的表达式会导致这样的结果吗？这样的表达式对刚度矩阵K又有什么影响？如果把Da选为某特定值，例如Da△x2＝2，使得对角元素几乎为零，你认为求解时会出现什么困难？

6．方法4：线性系统分析

MATLAB和COMSOL Multiphysics的核心是线性系统分析。本节将简要回顾线性算法理论——本科生工程数学中的“矩阵方程”就是这类典型问题。好消息是不需要你自己动手编程处理矩阵，这就是MATLAB的用途：提供一个工程矩阵计算子程序库的用户界面。应当注意到，对于本章的例子，COMSOL Script也能很好地实现这个功能。以前科学计算主要集中在矩阵计算的高效和松散方法上。Golub和Van Loan[3]的书是一本很好的专家级矩阵计算指南，但是对于MATLAB入门水平，Hanselman和Littlefield[4]的书是个不错的选择。

简单来讲，标准的矩阵方程如下：

a11x1a12x2a13x3a21x1a22x2a23x3a31x1a32x2a33x3aM1x1aM2x2aM3x3a1NxNb1a2NxNb2a3NxNb3aMNxNbM (26)

这里有N个与M方程有关的未知变量xj。系数aij和右侧常数项bi都是未知数。在工程应用中，通常将模型转化为这样的线性方程系统。例如，质量和能量守恒通常都会生成这样的线性方程系统。

求解可能性

当N＝M时，约束条件和未知数相同，所以可能求出唯一解xj。但是如果一个或者多个方程彼此线性相关（行退化），或者所有的方程都只含由某些特定变量组成的未知数（列退化），就不能得出唯一解。对于稀疏矩阵，行退化和列退化相同，退化方程会导致奇异。但是从数值求解讲，至少还有两点会导致错误：

 如果方程彼此不是严格线性相关，一些方程也可能在计算过程中由于舍

入误差而导致变得非常接近线性相关。

 求解过程中累积的舍入误差可能会“淹没”真实解，这在N比较大的时

候容易发生。这种情况下，计算过程没有出错，但是计算结果会不满足原始方程。 线性系统指南

没有“典型”的线性方程系统，但是一个粗略的想法就是舍入误差变得不可忽略：

 当N在20－50之间时，可以用单精度常规方法求解，不需要对以上两

种错误进行修正。

 当N在几百数量级时可以用双精度求解。

 当N在几千数量级时，如果系数矩阵稀疏（几乎都是零），由于系数特

性可以求解。对于稀疏矩阵，MATLAB有特殊的数据形式和一套针对性的函数。

但是在工程实际和物理过程中，有很多问题本身就是奇异或者非常接近奇异的。你可能会发现N＝10也很难求解。用奇异值分解法有时可以解决奇异值问题，将其分解为非奇异值。

数值线性代数的常规任务

方程(26)可以写为紧凑矩阵方程形式：

Axb (27)

 求解未知矢量x，这里A是方阵，b是已知矢量。  A为常数，求解多个b矢量时的解。  计算方阵A的逆矩阵A－1。  计算方阵A的行列式。

 如果M欠定系统。这种情况下可能没有解或者有多个解。解空间由特殊解xp乘以任何线性相关的N－M维矢量构成，称为矩阵A的化零空间。我们的任务是通过奇异值分解找出解空间。

 如果M>N，通常没有符合(26)的矢量解x。经常会出现这种过定系统，

我们的任务主要是寻找最接近于满足方程的折衷解。通常接近程度以(26)

方程左右两侧差值的“最小二乘法”来表示。 MATLAB矩阵计算

通过inv(矩阵)命令可以很容易获得逆矩阵。矩阵方程的解以矩阵除法运算符“\\”表示：

>> A=[ 3 -1 0; -1 6 -2; 0 -2 10]; >> B=[1; 5; 26]; >> X=A\\B X = 1.0000 2.0000 3.0000

行列式

行列式多用在稳定理论和评估矩阵奇异性上。为什么需要知道行列式？大多数时候你是想知道行列式是否为零。但是，当行列式为零、或者数值上非常接近于零时，由于前面提到的“舍入误差”原因很难进行数值计算。这是奇异值分解的另一个实例。

MATLAB用det(A)命令计算奇异值。在MATLAB命令行中手工输入以下矩阵，或者从matrix2.dat文件中复制以下矩阵：

>> A=[0.45, -0.244111, -0.0193373, 0.323972, -0.118829; -0.244111, 0.684036, -0.103427, 0.205569, 0.00292382; -0.0193373,-0.103427,0.8295, 0.0189674, -0.011169; 0.323972, 0.205569,0.0189674, 0.659479, 0.197388; -0.118829,0.00292382,-0.011169, 0.197388, 0.776985]

用det命令可以计算出行列式：

>> det(A) ans = -1.9682e-008

主轴理论：特征值和特征向量

MATLAB有预置的计算矩阵特征值和特征向量的函数：

>> eig(A) ans = -0.0000 0.7000 0.8000 0.9000 1.0000

当使用以下形式调用eig()函数时，它会以矩阵V的列向量形式返回特征向量：

>> [V,D]=eig(A) V =

-0.6836 0.0000 -0.5469 -0.4785 -0.0684 -0.4181 0.6162 0.1831 0.4530 -0.4547 -0.0837 0.4003 0.6189 -0.6232 0.2479 0.5409 0.2582 -0.2415 -0.4042 -0.6474 -0.2416 -0.6272 0.4755 -0.1190 -0.5550

D =

- 0.0000 0 0 0 0 0 0.7000 0 0 0 0 0 0.8000 0 0 0 0 0 0.9000 0 0 0 0 0 1.0000

eigs()函数是eig()函数的变形，它能够计算稀疏矩阵的特征值和特征向量。在下一节中将结合COMSOL Multiphysics介绍它的应用。

如果矩阵A的行列式 非常接近于零，那么一个特征值也非常接近于零。特征向量是矩阵A的化零空间——给出映射到零的方向：

>> A*V(:,1) ans = 1.0e-007 * 0.2669 0.1633 0.0327 -0.2112 0.0943

所有特征向量可以通过乘以特征值等于矩阵A的特性来证明这一点，例如：

>> A*V(:,2)./ V(:,2) ans = 0.7000 0.7000 0.7000 0.7000 0.7000

在MATLAB中，“./”符号表示元素对元素的除法，上面的冒号表示矩阵所有的列。

由于系统接近奇异，所有包含该矩阵的解都很差，对这一点我们不应该感到

奇怪。例如：

>> B=[0; 1; 0; 1; 0]; >> A\\B ans = 1.0e+006 * 2.1487 1.3142 0.2631 -1.7001 0.7593

由于A的所有元素都是个位数，被除数B矢量元素也是个位数，方程(27)的解也应该由个位数组成，而不是百万量级。对于化学工程师，这意味着对于一个质量守恒系统，入口流速为1kg/hr，由于质量守恒限制出口也为1kg/hr，但是反应器中流速为1000000kg/hr。这是不可能的。但这就是由近奇异矩阵给出的计算结果。

奇异值分解（SVD）

奇异值分解在很多情况下可以给出一个较好的解。类似于特征值和特征向量的主轴理论，所有矩阵都只有唯一的分解：

AUdiagVT

(28)

这里的U和V都是正交方阵，diag是包含奇异值的对角矩阵。在解出U，V和diag的基础上，(27)式可以轻易解得

TA1V1/diagU j (29)

U和V为正交矩阵意味着其转置矩阵也是其逆矩阵。对角阵的逆矩阵为对角元素

取倒数。所以求解的唯一问题就是何时一个或多个奇异值(diagj)接近于零。这时(1/diagj)是一个非常大的数，导致求解结果错误。一个比较好的近似方法是将这些奇异值的(1/diagj)设为零！

TxV1/diagUb j (30)

几乎为零元素的替代向量应该在数量级上最小，近似满足方程。

在我们矩阵A的例子中，如果只有一个输出，MATLAB命令svd()给出奇异值，如果有三个输出，命令svd()给出：

>>[U,D,V]=svd(A) U =

-0.0684 -0.4785 0.5469 0.0000 -0.6836 -0.4547 0.4530 -0.1831 -0.6162 -0.4181 0.2479 -0.6232 -0.6189 -0.4003 -0.0837 -0.6474 -0.4042 0.2415 -0.2582 0.5409 -0.5550 -0.1190 -0.4755 0.6272 -0.2416 D =

1.0000 0 0 0 0 0 0.9000 0 0 0 0 0 0.8000 0 0 0 0 0 0.7000 0 0 0 0 -0.0000 1.0000 V =

-0.0684 -0.4785 0.5469 0.0000 -0.6836 -0.4547 0.4530 -0.1831 -0.6162 -0.4181 0.2479 -0.6232 -0.6189 -0.4003 -0.0837 -0.6474 -0.4042 0.2415 -0.2582 0.5409 -0.5550 -0.1190 -0.4755 0.6272 -0.2416

由于A为对称矩阵，必然有U＝V。

最小数量级的奇异值分解命令如下：

>> ss=[1. 1./0.9 1./0.8 1./0.7 0]; >> dinv=diag(ss); >> V*dinv*U’*B ans = 0.0893 1.2820 0.1479 1.0317 -0.2130

这是一个物理上更容易接受的答案，例如，前面讨论的质量守恒系统中的内部气体流速。

由于有限元方法基于矩阵，所以线性系统理论对于COMSOL Multiphysics建模非常重要。当产生的刚度矩阵变得接近奇异时，COMSOL Multiphysics可能不能给出满足要求的结果。MATLAB中的矩阵计算和稀疏化方法可轻易地用于检查COMSOL Multiphysics结果的合理性，也可以用于对系统自然动力学的特征分析。这些想法将在下一节中作为COMSOL Multiphysics的一个模型例子加以介绍。

6.1 非均匀介质中的热传导

经常遇到的稳态热传导方程，是可以分析求解的最简单的偏微分方程：泊松方程。但是对于复杂几何情况，会比较难处理。作者认为某些通解很难推导和收敛[5]。所以，数值解比通解更好，进一步讲，传热过程的任何改变都可能会破坏分析结果。本节我们考虑非均匀厚板中的一维传热问题，在板中有分布热源：

ddTkf(x)dxdx Tx01Tx10 (31)

打开耦合MATLAB/COMSOL Script的COMSOL Multiphysics模型导航栏：

按照表9中的步骤建立有分布热源的介质传热问题。根据计算结果可以画出一条几乎线性的曲线，斜率为－1。由于这是个线性问题，求解器日志显示两步迭代收敛。计算结果得T|x=0.5=0.474097，而解析解为T|x=0.5～0.475：

13x3x4T1x

12612 (32)

下面尝试k＝1－x/2。这种情况下也有解析解，但是由于数值复杂性，需要采用

对数和分支切割。解析解为T|x=0.5～0.550，你的计算结果呢？

表9 非均匀介质中的热传导——分布热源情况 Model Navigator 选择1D空间维数，COMSOL Multiphysics: Heat Transfer: Conduction: Steady state analysis，设定因变量：u 选择Element：Lagrange-Linear，完成 Draw菜单 Specify objects: Line, Coordinates弹出菜单，x：0 1 名称：interval，完成 Physics菜单： 选择边界1 Boundary settings 设定边界条件为：温度；T0＝1 选择边界2 设定边界：T1＝0，完成 Physics菜单： 选择域1 Subdomain settings 设定k=1; Q=-x*(1-x) 选择Init选项卡；设定T(t0)=1-x，完成 Meshing 点击三角形按钮，划分网格 Solver 点击求解按钮（＝） Post-processing 设定x＝0.5 Data display 值：0.474097[K]，表达式：T，位置：(0.5) 下面来看线性系统理论。打开“File”菜单，选择“Export FEM structure as fem”。这是我们第二次导出为FEM结构，有必要简单介绍一下。对于任何MATLAB/COMSOL Script变量，在命令行中输入变量名字，MATLAB/COMSOL Script将会给出变量值或其数据结构。试一下：

>> fem fem =

version: [1x1 struct] appl: {[1x1 struct]} geom: [1x1 geom1] mesh: [1x1 femmesh] border: 1 outform: ’general’ form: ’general’ units: ’SI’ equ: [1x1 struct] bnd: [1x1 struct]

draw: [1x1 struct]

xmesh: [1x1 com.femlab.xmesh.Xmesh] sol: [1x1 femsol]

当然，你也可以进一步查看FEM结构，了解整个树的干、枝、叶。我们已经剪除了

fem.sol.u, fem.mesh.p1

不同的分支包含了完整的COMSOL Multiphysics模型——为fem.geom中保存的几何结构定义了方程(fem.appl{1}.equ)和边界条件(fem.appl{1}.equ)。与其它变量相同，这些结构可以被MATLAB/COMSOL Script编辑。以边界条件为例：

>> fem.appl{1}.bnd ans = type: ’T’ T0: {[0] [1]} ind: [2 1]

>> fem.appl{1}.bnd.T0{2} ans =

>> fem.appl{1}.bnd.T0{2}=10 ans = 10

现在可以用“File”菜单的导入功能将FEM结构载入到COMSOL Multiphysics中。之后检查一下“Physics”菜单边界设置——应当看到边界1为T＝10。显然，使用图形用户界面的下拉菜单可以更容易地实现该功能。但是在MATLAB中改变FEM结构是一个非常有用的功能。后面几章我们将学习如何使用该功能。但是在给定FEM数据结构复杂性的情况下，最好能够简化该功能。

在上一节中，已经学过建立刚度矩阵和使用eigs命令进行分析，现在你可以用MATLAB/COMSOL Script求解。由于K是系数矩阵，可以用eigs命令得到最小的六个特征值：

>>[K,L,M,N] = assemble(fem); >>dd=eig(K) >>dd=eigs(K, 6,-0.1) >>[V,D]= eigs(K, 6, -0.1)

注意到K有两个零特征值，并且所有特征值均非负。这点不需要担心，COMSOL Multiphysics通过块矩阵N和辅助向量M来建立其边界条件。可以通过取代K中行向量和L中元素来近似边界条件，但是在用N和M增大矩阵方程的同时，COMSOL Multiphysics允许使用多种通用型边界条件。实际上K奇异，因为对于使用Neumann边界条件（无流）的有限元方法，块矩阵是自然边界条件的结果。对于纯Neumann边界条件有无穷解，因为任意解加上任何值仍然是满足条件的解。我还是本科生时，第一次使用有限元方法时，K奇异就难住了我。

纯Neumann边界条件（本例中的零热流）是病态的，例如，如果把本例中的边界条件改为纯Neumann边界条件，你会发现得不到稳态解——导致矩阵奇异。为了得到唯一解必须定义一些参考值（Dirichlet边界条件）。但是通过奇异值分解或主轴理论，MATLAB能够求解这类问题。由于K是半正定矩阵，所有的特征值都是非负数。所以可以根据前几节介绍的伪逆方法计算含有零特征值K的伪逆：

>>ss=1./dd; >>ss(6)=0; >>dinv=diag(ss);

>>uneumann=V*dinv*V’*L

最后记得在FEMLAB的相同网格下进行求解。使用以下命令绘制计算结果曲线：

>>[xs, idx]=sort(fem.mesh.p); >>plot(xs, fem.sol.u(idx)) ;

类似的，投影在前5个最小非零特征值的特征向量得到的近似Neumann解如下：

>>plot(xs, uneumann(idx));

图7 非均匀和分布热源热传导问题的纯Neumann投影解

图7给出了使用伪逆方法投影在前6个模式下的解。你可以认为该曲线背离

了满足边界条件的线性函数，如1－x。进一步讲，特征向量也可以用相同方法获得：

>>col1=V(:,1); >>plot(xs,col1(idx)); >>col2=V(:,2); >>plot(xs,col2(idx));

我经常抱怨为什么在数学课上不讲解特征值和特征向量，很多同学不知道它们为什么要学这个，认为这个很难而放弃了。该问题中的特征向量表明了在很长时间内稳态结果的衰变“模式”。在这种情况下，由于COMSOL Multiphysics符号与通常情况相反，特征值就是（指数）衰变速率。负特征值就表示非稳态模式的增长速率。这不是该特性的所有分析——一些非稳态情况会引起有趣的物理现象（形成图案，爆炸等），但是也很难进行数学建模——计算模型不收敛。

图8 非均匀和分布热源热传导问题的最小特征值/特征向量

图9非均匀和分布热源热传导问题的倒数第二小特征值/特征向量

注意对于Neumann边界的解，加上任何常数仍是问题的解。特征向量没有归一化，乘以任何数仍为特征向量。图8和图9给出了不同数量级特征值和对应特征向量。这是最慢的衰减模式，所以预期的“驻波”随着稳态的到达而消失。

参考文献

[1] M. B. Cutlip and M. Shacham, Problem solving in chemical Engineering

with Numerical Methods (Prentice-Hall, Upper Saddle River, NJ, 1999). [2] W. C. Gear, Numerical Initial Value Problems in Ordinary Differential

Equations (1971).

[3] G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd edn. (Johns

Hopkins University Press, Baltimore, London, 1996).

[4] D. Hanselman and B. Littlefield, Mastering MATLAB 7: A Comprehensive

Tutorial and Referrence (Prentice-Hall, Saddle River, NJ, 2005).

[5] W. B. Zimmerman, On the resistance of a spherical particle settling in a tube

of viscous fluid, International Journal of Engineering Science 42(17-18) (2004) 1753-1778.