第一章 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所示步骤，建立一个零维应用模式来解决非线性多项式方程：
u?u?4u?2?0
32 (3)
通过在“Physics”菜单、“Subdomain settings”选项中的设定，使得表1中的每个子域都满足该方程。注意左上角，这是以矢量符号给出的方程。在一维模式下，这个方程可以简化为：
da?u?t????u?u?c??u??au???f???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中奇异刚度矩阵的问题可以通过一般模式来避免，它可以求解下式：
ea?u?t22设定：c=0; a=4;f=u^3+u^2+2; da=0 应用，选择Init选项卡：设定u(t0)=-2，完成 设定最大基元尺寸1 点击remesh，完成 求解器选择Stationary nonlinear，求解，完成 选择边界1，表达式：u，完成 ?da?u?t????x?F (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?1?1???1??Ki?? (7)
由于各yi相加为1，所以我们有?的非线性方程：
nf????1??i?1Xi?1?1???1??Ki???0 (8)
其中n是组分数量。该函数f(?)可用root(s)?求解，将计算结果回代就可以得到
所有产物的莫尔分数。牛顿-拉夫逊方法需要知道当前导数f'(?)来决定计算方向，在COMSOL Multiphysics中将解析计算该值。通过下式很容易得到牛顿-拉
夫逊迭代方法：