1、第一章 COMSOL MULTIPHYSICS 及数值分析基础W. B. J. ZIMMERMAN，B. N. HEWAKANDAMBYDepartment of Chemical and Process Engineering, University of Sheffield,Newcastle Street, Sheffield S1 3JD United KingdomE-mail: w.zimmermanshef.ac.uk本章主要介绍 COMSOL Multiphysics 在零维和一维模型数值分析方面的几个关键内容。这些内容包括求根、步进式数值积分、常微分方程数值积分和线性系统分析
2、。这几乎是所有的化工过程数学分析方法。下面通过 COMSOL Multiphysics 中的一些常见化工过程应用实例来介绍这些方法，包括：闪蒸、管式反应器设计、扩散反应系统和固体中热传导。1简介本章内容很多，可以分为几个不同的目标。首先介绍了 COMSOL Multiphysics 的主要工作特性；其次介绍了如何使用这些特性来分析一些简单的，位于零维空间、一维空间或“空间时间”系统中的化工问题。本章还希望通过展示 COMSOL Multiphysics 和 MATLAB 工具在化工过程分析中的强大功能，激发读者对使用 COMSOL Multiphysics 进行建模与仿真的兴趣。由于 COMS
3、OL Multiphysics 不是一个通用的问题求解工具，所以一些目标需要迂回实现。作者在使用 FORTRAN、Mathematica 和 MATLAB 解决化工问题方面有着丰富的教学经验，并用这些工具实现过这里所有的例子。而且，扩展化工问题的数值分析也已经在 POLYMATH1中实现，这似乎只在化工委员会的 CACHE 项目中使用过。本书前一版已经介绍过在零维空间中求解非线性代数方程和与时间有关的常微分方程的内容。从概念上讲，零维域就是一个简单的有限元。通过研究某一特定有限元中的变化对理解有限元方法非常有用。但是，COMSOL Multiphysics 通过独立对话框设置，使得零维几何方程
4、和与时间相关的常微分方程求解变得非常简单。所以本章将同时采用这两种方法求解这些例子。2方法 1：求根典型的数值分析课程会讲解多种求根方法，但是从实际经验来看，只有两种算法非常有用二分法和牛顿法。我们这里没有列出所有方法，而是重点考虑为什么求根是最有效的数值分析工具。在线性系统中求根非常简单，但是对于非线性系统这就是一个挑战，而所有感兴趣的动力学问题几乎都是非线性系统。对非线性系统的求根起源于对反函数的描述。为什么呢？因为对于大多数非线性函数， “正向”y= f(u)很好表示，但是它的反函数 u=f-1 (y)可能不能显式表示、多值（无意义）或根本不存在。如果反函数存在的话，求解反函数其实就是求
5、根的过程求解满足 F(u)=0 的 u 等价于求解 F(u)=f(u)-y=0。因为大多数数值分析的目标是在系统约束下计算求解，所以这也等价于对所有的约束取反。COMSOL Multiphysics 拥有求解非线性问题的核心函数femnlin ，本节主要介绍用它求解零维非线性问题。femnlin 函数使用牛顿方法求解，由于只有一个变量 u，牛顿法通过对一阶倒数 迭代来求根。该方法首先估计函数的斜率范围，然后再逼近根。该斜()Fu率可以通过理论分析（牛顿拉夫逊方法）和数值（正割法）方法求得。如果能用任何一种方法求得斜率，就可以用泰勒定律来逼近根。其基本思想就是使用目前猜测值 u0 的泰勒展开式：
6、(1)00()()(fuuf该公式可以化简，忽略(u-u 0)的高阶项，计算根如下：(2)0()f这个方法可以快速地扩展到多维求解空间，例如将 u 看作未知矢量， “被除”看作“乘以 f 的雅克比矩阵的逆” 。下一节介绍 COMSOL 0()fuMultiphysics 中的求根过程。2.1 求根：COMSOL Multiphysics 非线性求解器的应用实例如上节所述，求根本身是一个“零维”活动，至少对于“空间-时间”系统多维未知矢量 u 来说是这样的。COMSOL 多物理场没有零维模式，所以我们临时采用一维模式。这在方面增加了我们不需要的冗余功能。但是由于问题规模较小，COMSOL Mul
7、tiphysics 编码效率高，且现代微处理器的运算速度快，这点就不成为问题了。启动 MATLAB 并在命令窗口键入 COMSOL Multiphysics。屏幕闪烁几秒后，会出现一个模型导航窗口。按照表 1 所示步骤，建立一个零维应用模式来解决非线性多项式方程：(3)3240u通过在“Physics ”菜单、 “Subdomain settings”选项中的设定，使得表 1 中的每个子域都满足该方程。注意左上角，这是以矢量符号给出的方程。在一维模式下，这个方程可以简化为：(4)auudcuaftxx显然，、 和 在一维简化里是多余的。既然我们是求零维的根，所有公式 4左侧的系数都可以设为 0
8、。通过重新组合多项式，我们发现 a=4 和 f=u + 3u +2。注意我们是通过设定最大基元尺寸为 1、将域离散化为只有一个基元的2域，从而得到零维空间的！表 1 系数模式下的求根实例。文件名称：rootfinder.mph.Model Navigator 选择 1D 空间维数，COMSOL Multiphysics: PDE modes: PDE, coefficient 形式Draw 菜单 Specify objects：Line。Coordinates 弹出菜单。x：0 1名称：interval，完成Physics 菜单：Boundary settings选择域：1 和 2（按住 Ct
9、rl 键同时选中）选择 Neumann 边界条件采用默认设置 q=0 g=0，完成Physics 菜单：Subdomain settings选择域：1设定：c =0; a=4;f=u3+u2+2; =0ad应用，选择 Init 选项卡：设定 u(t0)=-2，完成Mesh 菜单：Mesh parameters设定最大基元尺寸 1点击 remesh，完成Solve 菜单：Solver parameters 求解器选择 Stationary nonlinear，求解，完成Post-processingPoint Evalution 选择边界 1，表达式：u，完成设定初始猜测值为 u(t )=-2，
10、下面我们来寻找最接近这个值的根。你可能对为0什么设置 a=4 而不是把所有有关参数放进 f 中感到奇怪，那是因为对公式(4)右侧有限元离散后得不到一个奇异刚度矩阵。后处理阶段在输出窗口中显示出结果：Value: -2.732051, Expression: u, Boundary: 1.解析解得出的根为-1- ，数值解能够很好的与其吻合。根据代数中二次3公式的解知道，另外一个根是-1+ ，第三个根是 1。回到子域设置中来，如果最初设置的猜测值为 u(t0)=-0.5，则 COMSOL Multiphysics 解出 u=0.762051，又一个很好的近似数值解。如果 u(t0)=1.2 作为最
11、初猜想时，得到 u=1。该练习体现了非线性求解器和问题的两个特性(1)非线性问题可以有多个解；(2) 最初猜测值对于求解很关键。牛顿法通常收敛到最接近的解，但是在高阶非线性问题中，这点可能不成立。这些特性在多维求解空间和“空间-时间”相关系统中同样适用。COMSOL Multiphysics 模型 mph 文件（rootfinder.mph）中包含了所有的MATLAB/COMSOL Script 源代码和 FEMLAB 的扩展功能，以便再次恢复到FEMLAB 当前图形用户界面。该文件在网页 http:/eyrie.shef.ac.uk/femlab 中可以得到。只需打开“Menu”菜单，选择“
12、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 中奇异刚度矩阵的问题可以通过一般模式来
13、避免，它可以求解下式：(5)2aauedFttx其中 (u,ux)和(4) 式中的系数项功能类似，但是被求解器处理方式不同。在系数模式中，认为系数与 u 无关，除非使用数值雅可比矩阵，这会产生一些非线性依赖关系迭代次数增加。一般模式下构建刚度矩阵时，精确雅可比矩阵会将 和 F 对 u 进行微分。通常一般模式比使用数值雅可比矩阵的系数模式需要的迭代次数少。下面使用的雅可比矩阵，不需要任何特殊处理来避免刚度矩阵奇异。总的来说，一般模式在处理非线性问题时比系数模式更加实用。从我的个人观点看来，系数模式是 COMSOL Multiphysics 的一个遗传特性MATLAB 的偏微分方程工具箱，它在很多
14、方面是 FEMLAB 和 COMSOL Multiphysics 的先驱，它广泛使用了系数形式。此外，数值雅克比矩阵系数公式是一个长期存在的标准有限元方法，所以相对于其它代码，它是一种非常有用的公式。表 2 给出了应用一般模式的步骤对我们刚才的步骤做了小小改动。表 2 一般模式下的求根实例。文件名：rtfingen.mphModel Navigator 1D, COMSOL Multiphysics: PDE modes, general 模式Options 设定 Axes/Grid 为0,1Draw 名称：Interval；Start0；Stop1Physics 菜单：Boundary Se
15、ttings 设定两个端点（域）为 Neumann 边界条件Physics 菜单：Subdomain Settings设定 =0; =0; F=u3+u2-4*u+2ad初值 u(t0)=-0.5Mesh 模式 设定最大基元尺寸，general 1；点击 RemeshSolve 使用默认设置（nonlinear solver ，exact Jacobian）Post-Processing 经过 5 次迭代，结果收敛。点击图片读出 u0.732051。改变初始条件找到另外两个根虽然这里给出了单变量简单函数求根的模板(rtfindgen.mph)，但是实际上MATLAB 有更简单的计算子程序，通过
16、内置函数 fzero 和函数声明来求根，现在 COMSOL Multiphysics 图形用户界面可以使用辅助方程中的辅助变量来求解几何约束，该变量称为状态变量。作为对一般模式求根模型的补充，使用状态变量 v 按照表 3 的步骤来求解一个非线性方程。表 3 在 ODE设置中求根Physics 菜单：ODE settings 名称：v 。方程：tan h(v)-v2+5，完成Solve 菜单：Solver parameters 选择 Stationary nonlinear 求解器，求解，完成Post-processing 选择 Point evaluation，边界 2，表达式：vReport
17、 窗口 值：2.008819，表达式：v，边界：2下一节我们将非线性求根方法应用到一个常见的化工实例中闪蒸，它用到了 COMSOL Multiphysics 图形用户界面的更多新特性。2.2 求根：闪蒸实例化学热力学中广泛存在求根法的应用实例，由于化学平衡和质量守衡的约束条件很充分，再加上状态方程，就产生了和问题中未知变量个数相同的约束条件。在本节中，我们以闪蒸为例来介绍简单的求根方法，构建只有一个自由度的系统，这里以相分率 作为未知变量。表 4 闪蒸单元的初始组分以及平衡时的分配系数 K组分 XI 65, 3.4bar 下的 Ki乙烷 0.0079 16.2丙烷 0.1281 5.2i-丁烷
18、 0.0849 2.6n-丁烷 0.2690 1.98i-正戊烷 0.0589 0.91n-正戊烷 0.1361 0.72己烷 0.3151 0.28液相碳氢化合物混合物通过闪蒸到 65、3.4 bar 状态。液相组成和每种组分闪蒸“K ”值见表 4。我们希望知道闪蒸过程中气相和液相产物的成分和液相的比例。表 4 给出了原始成分。组分 i 的质量平衡满足以下关系(6)(1)iiiXyxXi 是闪蒸前液体的莫尔分数，x i 是闪蒸后液体莫尔分数，y i 是蒸发产物莫尔分数，是液相产物与总供给摩尔流率的比值。平衡系数的定义为：K i= yi/ xi。将该方程消去质量平衡中的 xi，得到一个关于 y
19、i 和 Xi 的方程：(7)1ii iK由于各 yi 相加为 1，所以我们有 的非线性方程：(8)10nii iXf K其中 n 是组分数量。该函数 可用 root(s) 求解，将计算结果回代就可以得()f到所有产物的莫尔分数。牛顿-拉夫逊方法需要知道当前导数 来决定计算()f方向，在 COMSOL Multiphysics 中将解析计算该值。通过下式很容易得到牛顿 -拉夫逊迭代方法：(9)21inii iXKf现在再回到 COMSOL Multiphysics 求根计算。作为一个练习，我们用一般 PDE模式来计算求解。我们可以只载入 rootfinder.mph 或 rtfindgen.mp
20、h 并加以修改来计算这个问题，但是熟悉 COMSOL Multiphysics 功能也是非常重要的目的。表 5 闪蒸实例Model Navigator 选择 1D，COMSOL Multiphysics: PDE modes:PDE, general 形式Draw 菜单 Specify objects: Line，Coordinates 弹出菜单，x：0 1名称：interval，完成Options 菜单：Constants输入表 4 中数据X1，0.0079X2，0.1281X3，0.0849X4，0.2690X5，0.0589X6，0.1361X7，0.3151Options 菜单：Sca
21、lar Expressions为(8)式中右侧各项定义表达式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 菜单：Boundary settings选择域：1 和 2（按住 Ctrl 键）选择 Neumann 边界条件采用默认设置 q=0, g=0 完成Physics 菜单：Subdomain settings选择域：1设定 F=1+t1+t2+t7; =0ad选
22、择 Init 选项卡，设定 u(t0)=0.5Mesh 菜单：Mesh parameters设定最大基元尺寸为 1点击 Remesh，完成Solve 菜单：Solver parameters选择 Stationary nonlinear 求解器求解，完成Post-processing：Point Evaluation选择边界：1，表达式：u，完成Report 窗口，值：0.458509，表达式：u启动 COMSOL Multiphysics，打开模型导航栏。如果你已经打开 COMSOL Multiphysics 对话框，那么将你的工作空间保存为 mph 文件或者将命令保存为m 文件，然后打开“
23、File”菜单选择“New” 。按照表 5 的步骤建立闪蒸模型。注意本例中有两个新增内容“Options：Constants ”和“Options：Expressions” 。在这里可以定义常数，然后在 COMSOL Multiphysics中任何可以使用纯数字的输入框中使用定义的常数。表达式和常数类似，可以在任何 COMSOL Multiphysics 数字输入框中使用，但是不同之处在于表达式依赖于因变量。定义的常数和表达式同样可以在后处理过程中使用。后处理过程显示 t1 和 t2 为：Value: -0.013865, Expression: t1, Boundary:1Value: -0
24、.203441, Expression: t2, Boundary:1另外求解器日志中经常包含有用的信息。打开“solver”菜单，在最底部选择“View Log”，会弹出求解器日志对话框。它显示了求解器何时开始计算，执行了什么命令，如何从初始状态得到计算结果。这里经过三次迭代，绝对误差达到 109 。如果你的解不是不收敛或收敛很慢的话，很容易得到这样的信息。练习1.1 计算方程的根。由3251()0fuu于该函数是三次多项式，所以存在解析的无理数解，。 12uuu exp(u)-11.2 计算方程 的根。这是个超越方程，意味着不存在解析的无()10ufe理数解。如果你使用系数模式，设 c=1
25、 来帮助收敛。3. 方法 2：步进法数值积分数值积分是数值分析的支柱。在有计算机之前，科学计算的首要工作是制作特殊方程手册，其中几乎包括了所有特殊常微分方程的解。那么用什么方法计算得到的？就是用一维数值积分。一维数值积分分为两种：初值问题（IVP）和边界值问题（BVP） 。我们下一节再介绍边界值问题。最容易积分的是初值问题，因为所有的初始条件都在一点定义，可以根据一阶导数直接步进积分。显然，如果常微分方程是一阶的，例如：(10)()()ttdyfyft那么当 t0 时，(10) 中第二项严格成立。它满足欧拉法的条件，是积分一阶常微分方程最简单的办法。对于一维情况，可以简单地根据 f 在( xn
26、,yn)点的导数向前步进积分，这里 n 是指积分过程中的第 n 次离散化步骤。所以，(11)1(,)nnyhfxy这里假设导数变化不超过步进量 h，这只对线性函数才严格成立。对于任何有曲率的函数，该假设都不成立。下面考虑在图 1 中，如果步进量 h 取得过大会造成什么错误。显然，欧拉法需要改进的一点就是增大步进量，因为它需要小步进量以保证准确性。欧拉法也被称作“一阶”准确性，因为误差只能降低到一阶 h 的程度。图 1 数值积分中的欧拉步进龙格库塔方法所以如果我们希望增大步进量，那么就需要一个“更高阶的方法” ，随着步进量的减小，误差快速降低。k 阶方法误差大约在 hk 量级。假设我们忽略了曲率
27、，可以通过计算斜率 f(x)在 xn 和 xn+1 之间几个中间点的值来计算 y(x)的曲率。可以首先根据初始导数计算间隔中点的值，然后再用整个间隔中点的微分值来获得二阶准确性。(12)12 1312(,)()nnkhfyxkyO这样我们通过对两个函数的计算又得到一阶准确性。所以，如果使用一阶方法，N 次计算得到 O(1/N)误差，而二阶方法 2N 次计算得到 O(1/4N2)误差，如果使用一阶方法的话，这需要 N2 次计算。更高阶龙格库塔方法我们可以做的更好吗？显然，如果用三个中间点可以获得三阶准确性，用四个中间点可以获得四阶准确性等等。不过更高阶的方法需要更多的编程工作，所以也要考虑我们的
28、时间。但是从本质上讲，我们计算函数的 k 阶导数不一定光滑。可能在增加“近似准确性”的同时，高阶导数的整体误差会变得很大，如果是这样，每次连续步进都可能导致误差的急剧增大。这表明高阶方法没有低阶方法稳定。常微分方程通常选用四阶龙格库塔方法积分，这样程序比较紧凑，准确性高，也有很好的稳定性。其它方法在数值积分编程的时候，还需要注意两个问题：数值稳定性：即使使用高阶准确性的方法，你的积分结果与检验值偏差仍然很大，这可能就是因为你的数值方法不稳定。你可以通过降低步进量来获得数值稳定性，但是这也意味着更长的计算时间。如果你的计算量很大，计算时间太长，可以选用半隐式方法，例如预测校正法等。刚性系统：物理
29、机理起作用时，刚性系统通常会有两个完全不同的长度或时间尺度。刚性系统可能会出现上面提到的“系统非稳定性”现象，或者非物理现象的振荡。可以参考 Gear2的书来解决这个问题。3.1 数值积分简单实例高阶常微分方程通常通过降阶和步进法求解，考虑如下方程：(13)2()()dyqxr除非 q(x)和 r(x)是常数，那么你很幸运能够从大多数分析方法教材中找到答案。也有一些特殊的 q(x)和 r(x)可以求得解析解，但是最好能计算出各种情况下的数值解。为什么呢？因为为了搞清楚 y(x)，你需要画出曲线，得到解析解时你也需要花费一定的精力在作图上。那么应该怎么做？首先将上面的二阶系统降低为两个一阶系统：
30、 ()dzxrq上式的常微分方程能够仿照(11)和(12) 进行时间步进法数值积分。举个简单例子：(14)20dut降阶产生了两个一阶常微分方程：(15)121udt假设初始条件为 u1=1，u 2=0，可以建立零维空间系统来积分这一对常微分方程。打开 COMSOL Multiphysics，等待模型导航栏窗口。如果已经打开一个COMSOL Multiphysics 窗口，把工作空间保存为 mph 文件或者把命令保存为 m文件，然后在“File”菜单中选择“New”选项。按照表 6 的步骤进行设置。表 6 时间步进积分的简单实例Model Navigator 选择 1D， COMSOL Mul
31、tiphysics: PDE modes: PDE, general 模式，插入因变量：u1 u2Draw 菜单 Specify objects: Line。Coordinates 弹出菜单，x：0 1名称：interval，完成Physics 菜单：Boundary settings选择域：1 和 2（按住 Ctrl 键）选择 Neumann 边界条件保持默认选项 q=0, g=0 完成Physics 菜单：Subdomain settings选择域：1设定 F1u 2，F 2u 1选择 Init 选项卡，设定 u1(tw)=1Mesh 菜单：Mesh parameters设置最大基元尺寸为
32、 1点击 Remesh，完成Solve 菜单：Solver parameters选择 time dependent 求解器，在 General 选项卡中输入Times: linspace(0,2*pi,50)求解，完成Post-processing：Cross-section plotparametersPoint 选项卡，接受 u1 的默认设置General 选项卡，接受所有时间的默认设置完成这个例子给了我们两个独立变量 u1、u 2 和空间坐标 x。注意“Subdomain settings”窗口左上角方程中的矢量标志。由于是矢量变量，所有输入选项卡都是矢量( F)或者矩阵(d a)输入。
33、MATLAB 命令 linspace(0,2*pi,50)生成 50 个基元的矢量，给出从 0 到 2 的数据，u 1(t=2)=1.004414。假设解析解是 u1(t=2)=1，结果不够准确（0.4误差）。之前 FEMLAB 允许用户在 MATLAB 和 FEMLAB 中选择针对时间积分的预置的求解器。也允许用户在求解器变量对话框的通用选项卡中调整容错值（相对和绝对） 。将相对误差调整为 0.001，绝对误差调整为 0.0001，得到一个相对较好的计算结果 u1(t=2)=0.998027。注意累积的全局误差与求解方法准确性0.001 的量级一致。图 2 一个周期的 u1(t)图 3 一个
34、周期的 u2(t)图 2 和图 3 表明如果设定足够小的步进时间，FEMLAB 能够完成高精度的正弦和余弦函数数值积分。虽然我们一般认为正弦和余弦是“解析函数” ，但是通过比较，还是能清楚地看到解析函数和需要数值积分的函数是似是而非的没有比 Bessel 函数和椭圆方程更有解析性的函数。练习1.3 使用物理场菜单中的常微分方程设置和相同的初时条件，试着积分方程(15)，假设状态变量为 v1 和 v2。基于时间的常微分方程和代数方程的唯一区别就是必须使用符号代替 dv1/dt 和 v1t 等。4管式反应器数值积分本节介绍在化工管式反应器设计时，如何同时求解一阶非线性常微分方程组。通常管式反应器设
35、计的关键要素是反应器的长度。一个气态酒精脱水管式反应器，工作在 2bar 和 150下。反应方程式为：25242CHO已有试验表明反应速率可以表示如下，其中 CA 是酒精浓度(mol/L)，R 是酒精的消耗速率(mol/s/m 3)： 2.7013AR反应器直径为 0.05m，入口酒精流速为 10g/s。问：如何确定反应器的长度，进而获得各种酒精转化率？我们希望确定酒精出口摩尔分数分别为0.5，0.4，0.3，0.2 和 0.1 时的反应器长度。化工设计理论假设反应热很小，反应器内为柱状流，理想气体，可以确定反应流体由四个常微分方程控制，其中独立变量为 CA，C W（水浓度） ，V （速度）和
36、 x（反应器长度）：(16)1WdRttCdVxt最后一个方程说明表面速度使得反应器长度和反应物在反应器内的存在时间之间存在一个平衡。初始入口条件为：(17)0(0)()AWCVx方法从初时条件和化学计量系数可以看出，C WC E（酒精浓度） ，由于温度和压力设为常数，所以浓度 C 也为常数，可以看作是理想气体，有：(18)J8314kmolKpCT初始进口速度 V0 可以由给定流速、入口密度（酒精分子量为 46kg/kmol）和管的横截面积计算出。该方程需要对时间 t 数值积分到需要的酒精摩尔分数。使用较为简单的欧拉方法或者龙格库塔方法积分。读者也许注意到，积分 CA 可能不需要其它变量，但
37、是 CW，V 和 x 与时间 t严格相关。但是由于我们需要求出对应摩尔分数的 x 值，所以最好同时求解四个变量。部分结果计算得到的数值解为：(19)0.15628.43ACtx其中 CA/C 如图 4 所示。图 4 酒精浓度随时间 t变化曲线COMSOL Multiphysics 计算我们希望再次建立虚拟的零维模拟空间，这次使用四个独立变量。打开COMSOL Multiphysics，等待模型导航栏窗口。根据表 7 的步骤建立管式反应器模型。该模型给出四个独立变量 u1、u 2、u 3、u 4 和一个空间坐标 x。由于变量比较多，建议用变量名称进行常数设定。进一步来说，由于使用了速度定律，用矢
38、量表达式更为方便。表 7 COMSOL Multiphysics中管式反应器设计建模步骤Model Navigator 选择 1D 空间维数，COMSOL Multiphysics: PDE modes: PDE, general 形式设定独立变量：u 1 u2 u3 u4 选择 Element：Lagrange-Linear，完成Draw 菜单 Specify objects: Line，Coordinates 弹出菜单，x ：0 1名称：interval，完成Options 菜单：Constants如下填写变量表格：名称 ExpressionP 200000T 423R 8314MM 46
39、Flowrate 0.01Dia 0.05C P/(RT)area pi*Dia2/4rho MM*Cvel Flowrate/rho/aera完成Options 菜单：Scalar Expressions定义 rate=52.7*u12/(1+0.013/u1)Physics 菜单：Boundary settings选择域：1 和 2（按住 Ctrl 键）选择 Neumann 边界条件默认设置 q0 g=0，完成Physics 菜单：Subdomain settings选择域 1F 选项卡；设定 F1=-rate*(1+u1/C); F2=rate*(1-u2/C);F3=rate*u3/C
40、; F4=u3/CInit 选项卡；设定 u1(t0)=C; u3(t0)=vel，完成Mesh 菜单：Mesh parameters设定最大基元尺寸为 1点击 Remesh，完成Solve 菜单：Solver parameters选择 time dependent 求解器，在 General 选项卡中输入Times: linspace(0,10,100)求解，完成Post-processingCross-section plotparametersPoint 选项卡，接受 u1 的默认设置General 选项卡，接受所有时间的默认设置完成试着在整个时间范围内画出 u1，u 2，u 3 和 u
41、4 的曲线，与图 4 是否保持一致？与之前计算结果是否相符？练习1.4 计算以下方程组中 (x=1)的值，并画出 随 x 变化曲线，x 取 0 到 3 区间。yy20()1yx1.5 连续搅拌反应器中的一阶可逆反应系统可以用线性常微分方程组表示。例如，考虑如下异构反应： ABC正向反应速率为 k1 和 k3，相应逆向反应速率为 k2 和 k4。一阶动力学常微分方程组为：(20)1234ABBCCBdckctkdckct你也许会感到奇怪，由于该方程组为线性，它有通用的解析解。虽然是通用解析解，但是对动力学系统的反映并不多。假设刚开始时，CA 1，k 1=1Hz，k 2=0Hz，k 3=2Hz，k
42、 4=3Hz。针对该初值问题，画出浓度随时间变化的曲线。5. 方法 3：常微分方程数值积分前一节介绍了步进法数值积分，通常也称作“时间步进” ，但是在反应器设计中，多是空间积分问题。步进法是顺序求解未知项，而其它常见积分方法是同时求解某一网格点的所有未知独立变量。步进法只能求解初值问题（IVP） ，初值个数必须严格等于常微分方程组的阶数。但是对于二阶或者更高阶的系统，可能会用到第二类边界条件边界值问题，在一维情况下，这些边界条件只出现在域的起点和终点，这就是两点边界值问题。步进法也能够勉强求解边界值问题人为设定初值，通过尝试和误差修正找到满足边界值问题的初始条件。对于更高维数的偏微分方程，边界
43、值问题发生在整个域的每一个边界上。有限元方法的一个最主要优点就是能够很容易地求解两点边界值问题。例如，一维反应和扩散方程：(21)2()DuRLx这里 u 是组分浓度，D 是扩散系数，L 是域的长度，R (u)是反应消耗速率，x 是无量纲空间坐标。如果未知函数 u(x)在 x=xj=jx 网格点能够用不连续值 uj=u(xj)近似，那么根据中心差分，方程变为：(22)21Nij ijLMD其中 Rj=R(uj)，M ij 是对角元素为-2 ，上、下对角元素为 1 的三对角矩阵：(23)201 。通过矩阵转置和对 uin 进行积分可以求解该方程，这里 n 指第 n 次()jjRu猜测：(24)2
44、1Nni ijjLxuMRDRj=R(uin-1)。无论是初值问题还是边界值问题，都可以将(23)式中矩阵 M 的行向量变得符合边界条件。(23)式假设在 x0 和 x1 处都有 u0。这是狄利克雷边界条件，也是有限微分法的自然边界条件说它自然是因为如果在设定边界条件时没有改变(23) 式中的行向量，那么就会出现这样的边界条件。下面介绍如何用 COMSOL Multiphysics 在一维域上求解方程 (21)，对于一阶反应 R(u)=ku，典型的无量纲变量 Damkohler 数为：(25)231312900.83./smkLDas边界条件为 x0 时 u1 ，在 u1 处没有流量。下面结合
45、 MATLAB 来做这个练习，进而了解 COMSOL Multiphysics 如何表示计算数据和模型输出的结构。在 windows 中，结合 MATLAB 的 COMSOL Multiphysics 是一个可选的桌面图标，在 UNIX/linux 中，该功能可以通过在命令栏中执行以下语句实现：Comsol matlab path-ml nodesktop-ml nosplash其中“matlab”告诉 femlab 打开一个 matlab 命令窗口， “path”根据 COMSOL命令库建立 matlab 命令窗口。首先运行 COMSOL Multiphysics 和模型导航栏窗口，然后根据
46、表 8 中的步骤建立模型。该实例给了一个独立变量 u 和一个一维空间坐标 x。h 和 r 是系数模式中狄利克雷边界条件的两个句柄。如果你希望边界速度 u 为常数 U0，可以通过设定 h1，rU 0 来实现。点击三角形按钮生成默认网格（15 个） ，然后点击两个嵌套三角形按钮加密网格（30 个） 。表 8 反应扩散系统中两点边界值问题的常微分方程实例Model Navigator选择 1D 空间维数，COMSOL Multiphysics: PDE modes: PDE, coefficient 形式设置因变量：u选择 Element： Lagrange-Linear，完成Draw 菜单 Spe
47、cify objects: Line， Coordinates 弹出菜单，x：0 1名称：interval，完成Physics 菜单：Boundary settings选择域 1选择 Dirichlet 边界条件，设定 h=1; r=1选择域 2选择 Neumann 边界条件，完成Physics 菜单：Subdomain settings选择域 1设定 C=-1; f=0.833*u; =0ad选择 Init 选项卡；设定 u(t0)=1-x完成Meshing 点击三角形按钮进行网格划分Solve 菜单：Solver parametersStationary linear 求解器默认设置，完成
48、General 选项卡，设定解的形式为“Coefficient”求解（）Post-processing：Point evalution选择边界 2，输入表达式 u。Reports 窗口显示：值：0.861167，表达式：u，边界：2根据计算结果可以画出图 5，显然边界条件需要满足：x0 处 u1，在x1 处斜率减为零。但是 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 依赖于任何非独立变量，就应该选择非线性求解器。线性
