图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)：
R?52.7CA1?0.013CA2
反应器直径为0.05m，入口酒精流速为10g/s。问：如何确定反应器的长度，进而获得各种酒精转化率？我们希望确定酒精出口摩尔分数分别为0.5，0.4，0.3，
0.2和0.1时的反应器长度。
化工设计理论
假设反应热很小，反应器内为柱状流，理想气体，可以确定反应流体由四个常微分方程控制，其中独立变量为CA，CW（水浓度），V（速度）和x（反应器长度）：
dCAC????R?1?A?dtC??C???R?1?W?dtC?? dVRV?dtCdxdt?VdCW (16)
最后一个方程说明表面速度使得反应器长度和反应物在反应器内的存在时
间之间存在一个平衡。初始入口条件为：
CA(0)?CCW(0)?0V(0)?V0x(0)?0 (17)
方法
从初时条件和化学计量系数可以看出，CW＝CE（酒精浓度），由于温度和压力设为常数，所以浓度C也为常数，可以看作是理想气体，有：
C?p??T?8314?kmolK??J (18)
初始进口速度V0可以由给定流速、入口密度（酒精分子量为46kg/kmol）和管的
横截面积计算出。该方程需要对时间t数值积分到需要的酒精摩尔分数。使用较为简单的欧拉方法或者龙格库塔方法积分。
读者也许注意到，积分CA可能不需要其它变量，但是CW，V和x与时间t严格相关。但是由于我们需要求出对应摩尔分数的x值，所以最好同时求解四个变量。
部分结果
计算得到的数值解为：
CAC?0.1t?5.65225x?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'?y?0y(x?0)?1y'(x?0)?02选择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。一阶动力学常微分方程组为：
dcAdtdcBdt??k1cA?k2cB?k1cA?k2cB?k3cB?k4cC dcCdt?k3cB?k4cC (20)
你也许会感到奇怪，由于该方程组为线性，它有通用的解析解。虽然是通用解析解，但是对动力学系统的反映并不多。假设刚开始时，CA＝1，k1=1Hz，k2=0Hz，k3=2Hz，k4=3Hz。针对该初值问题，画出浓度随时间变化的曲线。
5. 方法3：常微分方程数值积分
前一节介绍了步进法数值积分，通常也称作“时间步进”，但是在反应器设计中，多是空间积分问题。步进法是顺序求解未知项，而其它常见积分方法是同时求解某一网格点的所有未知独立变量。步进法只能求解初值问题（IVP），初
值个数必须严格等于常微分方程组的阶数。但是对于二阶或者更高阶的系统，可能会用到第二类边界条件——边界值问题，在一维情况下，这些边界条件只出现在域的起点和终点，这就是两点边界值问题。步进法也能够勉强求解边界值问题——人为设定初值，通过尝试和误差修正找到满足边界值问题的初始条件。对于更高维数的偏微分方程，边界值问题发生在整个域的每一个边界上。
有限元方法的一个最主要优点就是能够很容易地求解两点边界值问题。例如，一维反应和扩散方程：
D?uL?x222?R(u) (21)
这里u是组分浓度，D是扩散系数，L是域的长度，R(u)是反应消耗速率，x是无量纲空间坐标。如果未知函数u(x)在x=xj=j△x网格点能够用不连续值uj=u(xj)近似，那么根据中心差分，方程变为：
N?Mj?1ijuj?L?xD22Ri
(22)
其中Rj=R(uj)，Mij是对角元素为-2，上、下对角元素为1的三对角矩阵：
??2?1?M??0??0???Rj?R(uj)。通过矩阵转置和对ui
n
1?210?01???00?????????? ?????? (23)
进行积分可以求解该方程，这里n指第n次猜
测：
u?niL?xD22N?Mj?1?1ijRj
(24)
Rj=R(uin-1)。无论是初值问题还是边界值问题，都可以将(23)式中矩阵M的行向量变得符合边界条件。(23)式假设在x＝0和x＝1处都有u＝0。这是狄利克雷边界条件，也是有限微分法的自然边界条件——说它自然是因为如果在设定边界条件时没有改变(23)式中的行向量，那么就会出现这样的边界条件。
下面介绍如何用COMSOL Multiphysics在一维域上求解方程(21)，对于一阶反应R(u)=ku，典型的无量纲变量Damkohler数为：
Da?kLD2??10s?3?1??10m?3?92?1?21.2?10m/s?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