COMSOL MULTIPHYSICS和数值分析基础.docx
《COMSOL MULTIPHYSICS和数值分析基础.docx》由会员分享，可在线阅读，更多相关《COMSOL MULTIPHYSICS和数值分析基础.docx（42页珍藏版）》请在冰点文库上搜索。
COMSOLMULTIPHYSICS和数值分析基础
第一章COMSOLMULTIPHYSICS及数值分析基础
W.B.J.ZIMMERMAN，B.N.HEWAKANDAMBY
DepartmentofChemicalandProcessEngineering,UniversityofSheffield,
NewcastleStreet,SheffieldS13JDUnitedKingdom
E-mail:
w.zimmerman@shef.ac.uk
本章主要介绍COMSOLMultiphysics在零维和一维模型数值分析方面的几个关键内容。
这些内容包括求根、步进式数值积分、常微分方程数值积分和线性系统分析。
这几乎是所有的化工过程数学分析方法。
下面通过COMSOLMultiphysics中的一些常见化工过程应用实例来介绍这些方法，包括：
闪蒸、管式反应器设计、扩散反应系统和固体中热传导。
1．简介
本章内容很多，可以分为几个不同的目标。
首先介绍了COMSOLMultiphysics的主要工作特性；其次介绍了如何使用这些特性来分析一些简单的，位于零维空间、一维空间或“空间－时间”系统中的化工问题。
本章还希望通过展示COMSOLMultiphysics和MATLAB工具在化工过程分析中的强大功能，激发读者对使用COMSOLMultiphysics进行建模与仿真的兴趣。
由于COMSOLMultiphysics不是一个通用的问题求解工具，所以一些目标需要迂回实现。
作者在使用FORTRAN、Mathematica和MATLAB解决化工问题方面有着丰富的教学经验，并用这些工具实现过这里所有的例子。
而且，扩展化工问题的数值分析也已经在POLYMATH[1]中实现，这似乎只在化工委员会的CACHE项目中使用过。
本书前一版已经介绍过在零维空间中求解非线性代数方程和与时间有关的常微分方程的内容。
从概念上讲，零维域就是一个简单的有限元。
通过研究某一特定有限元中的变化对理解有限元方法非常有用。
但是，COMSOLMultiphysics通过独立对话框设置，使得零维几何方程和与时间相关的常微分方程求解变得非常简单。
所以本章将同时采用这两种方法求解这些例子。
2．方法1：
求根
典型的数值分析课程会讲解多种求根方法，但是从实际经验来看，只有两种算法非常有用——二分法和牛顿法。
我们这里没有列出所有方法，而是重点考虑为什么求根是最有效的数值分析工具。
在线性系统中求根非常简单，但是对于非线性系统这就是一个挑战，而所有感兴趣的动力学问题几乎都是非线性系统。
对非线性系统的求根起源于对反函数的描述。
为什么呢？
因为对于大多数非线性函数，“正向”y=f（u）很好表示，但是它的反函数u=f-－1（y）可能不能显式表示、多值（无意义）或根本不存在。
如果反函数存在的话，求解反函数其实就是求根的过程——求解满足F（u）=0的u等价于求解F（u）=f（u）-y=0。
因为大多数数值分析的目标是在系统约束下计算求解，所以这也等价于对所有的约束取反。
COMSOLMultiphysics拥有求解非线性问题的核心函数——femnlin，本节主要介绍用它求解零维非线性问题。
femnlin函数使用牛顿方法求解，由于只有一个变量u，牛顿法通过对一阶倒数
迭代来求根。
该方法首先估计函数的斜率范围，然后再逼近根。
该斜率可以通过理论分析（牛顿－拉夫逊方法）和数值（正割法）方法求得。
如果能用任何一种方法求得斜率，就可以用泰勒定律来逼近根。
其基本思想就是使用目前猜测值u0的泰勒展开式：
（1）
该公式可以化简，忽略（u-u0）的高阶项，计算根如下：
（2）
这个方法可以快速地扩展到多维求解空间，例如将u看作未知矢量，“被
除”看作“乘以f的雅克比矩阵的逆”。
下一节介绍COMSOLMultiphysics中的求根过程。
2.1求根：
COMSOLMultiphysics非线性求解器的应用实例
如上节所述，求根本身是一个“零维”活动，至少对于“空间-时间”系统多维未知矢量u来说是这样的。
COMSOL多物理场没有零维模式，所以我们临时采用一维模式。
这在方面增加了我们不需要的冗余功能。
但是由于问题规模较小，COMSOLMultiphysics编码效率高，且现代微处理器的运算速度快，这点就不成为问题了。
启动MATLAB并在命令窗口键入COMSOLMultiphysics。
屏幕闪烁几秒后，会出现一个模型导航窗口。
按照表1所示步骤，建立一个零维应用模式来解决非线性多项式方程：
（3）
通过在“Physics”菜单、“Subdomainsettings”选项中的设定，使得表1中的每个子域都满足该方程。
注意左上角，这是以矢量符号给出的方程。
在一维模式下，这个方程可以简化为：
（4）
显然，α、γ和β在一维简化里是多余的。
既然我们是求零维的根，所有公式4左侧的系数都可以设为0。
通过重新组合多项式，我们发现a=4和f=u
+u
+2。
注意我们是通过设定最大基元尺寸为1、将域离散化为只有一个基元的域，从而得到零维空间的！
表1系数模式下的求根实例。
文件名称：
rootfinder.mph.
ModelNavigator
选择1D空间维数，COMSOLMultiphysics:
PDEmodes:
PDE,coefficient形式
Draw菜单
Specifyobjects：
Line。
Coordinates弹出菜单。
x：
01
名称：
interval，完成
Physics菜单：
Boundarysettings
选择域：
1和2（按住Ctrl键同时选中）
选择Neumann边界条件
采用默认设置q=0g=0，完成
Physics菜单：
Subdomainsettings
选择域：
1
设定：
c=0;a=4;f=u^3+u^2+2;
=0
应用，选择Init选项卡：
设定u（t0）=-2，完成
Mesh菜单：
Meshparameters
设定最大基元尺寸1
点击remesh，完成
Solve菜单：
Solverparameters
求解器选择Stationarynonlinear，求解，完成
Post-processing
PointEvalution
选择边界1，表达式：
u，完成
设定初始猜测值为u（t
）=-2，下面我们来寻找最接近这个值的根。
你可能对为什么设置a=4而不是把所有有关参数放进f中感到奇怪，那是因为对公式（4）右侧有限元离散后得不到一个奇异刚度矩阵。
后处理阶段在输出窗口中显示出结果：
Value:
-2.732051,Expression:
u,Boundary:
1.
解析解得出的根为-1-
，数值解能够很好的与其吻合。
根据代数中二次公式的解知道，另外一个根是-1+
，第三个根是1。
回到子域设置中来，如果最初设置的猜测值为u（t0）=-0.5，则COMSOLMultiphysics解出u=0.762051，又一个很好的近似数值解。
如果u（t0）=1.2作为最初猜想时，得到u=1。
该练习体现了非线性求解器和问题的两个特性——
（1）非线性问题可以有多个解；
（2）最初猜测值对于求解很关键。
牛顿法通常收敛到最接近的解，但是在高阶非线性问题中，这点可能不成立。
这些特性在多维求解空间和“空间-时间”相关系统中同样适用。
COMSOLMultiphysics模型mph文件（rootfinder.mph）中包含了所有的MATLAB/COMSOLScript源代码和FEMLAB的扩展功能，以便再次恢复到FEMLAB当前图形用户界面。
该文件在网页http:
//eyrie.shef.ac.uk/femlab中可以得到。
只需打开“Menu”菜单，选择“Openmodelm-file”选项，在弹出的“Openfile”对话框中选中它，你就可以在“Subdomainsettings”中快速设定非线性函数，给定初始猜测值，然后用非线性求解器得到收敛解。
但是如果函数中没有线性项可以放在公式（4）左侧怎么办？
例如，FEMLAB将tan（u）-u2+5=0放在公式（4）左侧时，会生成奇异刚度矩阵。
建议在子域中设定一个u的二次派生系数，c=1。
通过与Neumann边界条件相结合，这个人为扩散因素不能改变基元中解为常数的事实，但是它可以避免刚度矩阵奇异。
在一般模式下求根
tan（u）-u2+5=0中奇异刚度矩阵的问题可以通过一般模式来避免，它可以求解下式：
（5）
其中Г（u,ux）和（4）式中的系数项功能类似，但是被求解器处理方式不同。
在系数模式中，认为系数与u无关，除非使用数值雅可比矩阵，这会产生一些非线性依赖关系——迭代次数增加。
一般模式下构建刚度矩阵时，精确雅可比矩阵会将Г和F对u进行微分。
通常一般模式比使用数值雅可比矩阵的系数模式需要的迭代次数少。
下面使用的雅可比矩阵，不需要任何特殊处理来避免刚度矩阵奇异。
总的来说，一般模式在处理非线性问题时比系数模式更加实用。
从我的个人观点看来，系数模式是COMSOLMultiphysics的一个遗传特性——MATLAB的偏微分方程工具箱，它在很多方面是FEMLAB和COMSOLMultiphysics的先驱，它广泛使用了系数形式。
此外，数值雅克比矩阵系数公式是一个长期存在的标准有限元方法，所以相对于其它代码，它是一种非常有用的公式。
表2给出了应用一般模式的步骤——对我们刚才的步骤做了小小改动。
表2一般模式下的求根实例。
文件名：
rtfingen.mph
ModelNavigator
1D,COMSOLMultiphysics:
PDEmodes,general模式
Options
设定Axes/Grid为[0,1]
Draw
名称：
Interval；Start＝0；Stop＝1
Physics菜单：
BoundarySettings
设定两个端点（域）为Neumann边界条件
Physics菜单：
SubdomainSettings
设定Г=0;
=0;F=u3+u2-4*u+2
初值u（t0）=-0.5
Mesh模式
设定最大基元尺寸，general＝1；点击Remesh
Solve
使用默认设置（nonlinearsolver，exactJacobian）
Post-Processing
经过5次迭代，结果收敛。
点击图片读出u＝0.732051。
改变初始条件找到另外两个根
虽然这里给出了单变量简单函数求根的模板（rtfindgen.mph），但是实际上MATLAB有更简单的计算子程序，通过内置函数fzero和函数声明来求根，现在COMSOLMultiphysics图形用户界面可以使用辅助方程中的辅助变量来求解几何约束，该变量称为状态变量。
作为对一般模式求根模型的补充，使用状态变量v按照表3的步骤来求解一个非线性方程。
表3在ODE设置中求根
Physics菜单：
ODEsettings
名称：
v。
方程：
tanh（v）-v^2+5，完成
Solve菜单：
Solverparameters
选择Stationarynonlinear求解器，求解，完成
Post-processing
选择Pointevaluation，边界2，表达式：
v
Report窗口
值：
－2.008819，表达式：
v，边界：
2
下一节我们将非线性求根方法应用到一个常见的化工实例中——闪蒸，它用到了COMSOLMultiphysics图形用户界面的更多新特性。
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.4bar状态。
液相组成和每种组分闪蒸“K”值见表4。
我们希望知道闪蒸过程中气相和液相产物的成分和液相的比例。
表4给出了原始成分。
组分i的质量平衡满足以下关系
（6）
Xi是闪蒸前液体的莫尔分数，xi是闪蒸后液体莫尔分数，yi是蒸发产物莫尔分数，
是液相产物与总供给摩尔流率的比值。
平衡系数的定义为：
Ki=yi/xi。
将该方程消去质量平衡中的xi，得到一个关于yi和Xi的方程：
（7）
由于各yi相加为1，所以我们有
的非线性方程：
（8）
其中n是组分数量。
该函数
可用root（s）
求解，将计算结果回代就可以得到所有产物的莫尔分数。
牛顿-拉夫逊方法需要知道当前导数
来决定计算方向，在COMSOLMultiphysics中将解析计算该值。
通过下式很容易得到牛顿-拉夫逊迭代方法：
（9）
现在再回到COMSOLMultiphysics求根计算。
作为一个练习，我们用一般PDE模式来计算求解。
我们可以只载入rootfinder.mph或rtfindgen.mph并加以修改来计算这个问题，但是熟悉COMSOLMultiphysics功能也是非常重要的目的。
表5闪蒸实例
ModelNavigator
选择1D，COMSOLMultiphysics:
PDEmodes:
PDE,general形式
Draw菜单
Specifyobjects:
Line，Coordinates弹出菜单，x：
01
名称：
interval，完成
Options菜单：
Constants
输入表4中数据
X1，0.0079
X2，0.1281
X3，0.0849
X4，0.2690
X5，0.0589
X6，0.1361
X7，0.3151
Options菜单：
ScalarExpressions
为（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菜单：
Boundarysettings
选择域：
1和2（按住Ctrl键）
选择Neumann边界条件
采用默认设置q=0,g=0完成
Physics菜单：
Subdomainsettings
选择域：
1
设定F=1+t1+t2+…+t7;
=0
选择Init选项卡，设定u（t0）=0.5
Mesh菜单：
Meshparameters
设定最大基元尺寸为1
点击Remesh，完成
Solve菜单：
Solverparameters
选择Stationarynonlinear求解器求解，完成
Post-processing：
PointEvaluation
选择边界：
1，表达式：
u，完成
Report窗口，值：
0.458509，表达式：
u
启动COMSOLMultiphysics，打开模型导航栏。
如果你已经打开COMSOLMultiphysics对话框，那么将你的工作空间保存为mph文件或者将命令保存为m文件，然后打开“File”菜单选择“New”。
按照表5的步骤建立闪蒸模型。
注意本例中有两个新增内容——“Options：
Constants”和“Options：
Expressions”。
在这里可以定义常数，然后在COMSOLMultiphysics中任何可以使用纯数字的输入框中使用定义的常数。
表达式和常数类似，可以在任何COMSOLMultiphysics数字输入框中使用，但是不同之处在于表达式依赖于因变量。
定义的常数和表达式同样可以在后处理过程中使用。
后处理过程显示t1和t2为：
Value:
-0.013865,Expression:
t1,Boundary:
1
Value:
-0.203441,Expression:
t2,Boundary:
1
另外求解器日志中经常包含有用的信息。
打开“solver”菜单，在最底部选择“ViewLog”，会弹出求解器日志对话框。
它显示了求解器何时开始计算，执行了什么命令，如何从初始状态得到计算结果。
这里经过三次迭代，绝对误差达到10－9。
如果你的解不是不收敛或收敛很慢的话，很容易得到这样的信息。
练习
1.1计算方程
的根。
由于该函数是三次多项式，所以存在解析的无理数解，
。
uexp（u）-1
1.2计算方程
的根。
这是个超越方程，意味着不存在解析的无理数解。
如果你使用系数模式，设c=1来帮助收敛。
3.方法2：
步进法数值积分
数值积分是数值分析的支柱。
在有计算机之前，科学计算的首要工作是制作特殊方程手册，其中几乎包括了所有特殊常微分方程的解。
那么用什么方法计算得到的？
就是用一维数值积分。
一维数值积分分为两种：
初值问题（IVP）和边界值问题（BVP）。
我们下一节再介绍边界值问题。
最容易积分的是初值问题，因为所有的初始条件都在一点定义，可以根据一阶导数直接步进积分。
显然，如果常微分方程是一阶的，例如：
（10）
那么当△t→0时，（10）中第二项严格成立。
它满足欧拉法的条件，是积分一阶常微分方程最简单的办法。
对于一维情况，可以简单地根据f在（xn,yn）点的导数向前步进积分，这里n是指积分过程中的第n次离散化步骤。
所以，
（11）
这里假设导数变化不超过步进量h，这只对线性函数才严格成立。
对于任何有曲率的函数，该假设都不成立。
下面考虑在图1中，如果步进量h取得过大会造成什么错误。
显然，欧拉法需要改进的一点就是增大步进量，因为它需要小步进量以保证准确性。
欧拉法也被称作“一阶”准确性，因为误差只能降低到一阶h的程度。
图1数值积分中的欧拉步进
龙格库塔方法
所以如果我们希望增大步进量，那么就需要一个“更高阶的方法”，随着步进量的减小，误差快速降低。
k阶方法误差大约在hk量级。
假设我们忽略了曲率，可以通过计算斜率f（x）在xn和xn+1之间几个中间点的值来计算y（x）的曲率。
可以首先根据初始导数计算间隔中点的值，然后再用整个间隔中点的微分值来获得二阶准确性。
（12）
这样我们通过对两个函数的计算又得到一阶准确性。
所以，如果使用一阶方法，N次计算得到O（1/N）误差，而二阶方法2N次计算得到O（1/4N2）误差，如果使用一阶方法的话，这需要N2次计算。
更高阶龙格库塔方法
我们可以做的更好吗？
显然，如果用三个中间点可以获得三阶准确性，用四个中间点可以获得四阶准确性等等。
不过更高阶的方法需要更多的编程工作，所以也要考虑我们的时间。
但是从本质上讲，我们计算函数的k阶导数不一定光滑。
可能在增加“近似准确性”的同时，高阶导数的整体误差会变得很大，如果是这样，每次连续步进都可能导致误差的急剧增大。
这表明高阶方法没有低阶方法稳定。
常微分方程通常选用四阶龙格库塔方法积分，这样程序比较紧凑，准确性高，也有很好的稳定性。
其它方法
在数值积分编程的时候，还需要注意两个问题：
数值稳定性：
即使使用高阶准确性的方法，你的积分结果与检验值偏差仍然很大，这可能就是因为你的数值方法不稳定。
你可以通过降低步进量来获得数值稳定性，但是这也意味着更长的计算时间。
如果你的计算量很大，计算时间太长，可以选用半隐式方法，例如预测校正法等。
刚性系统：
物理机理起作用时，刚性系统通常会有两个完全不同的长度或时间尺度。
刚性系统可能会出现上面提到的“系统非稳定性”现象，或者非物理现象的振荡。
可以参考Gear[2]的书来解决这个问题。
3.1数值积分简单实例
高阶常微分方程通常通过降阶和步进法求解，考虑如下方程：
（13）
除非q（x）和r（x）是常数，那么你很幸运能够从大多数分析方法教材中找到答案。
也有一些特殊的q（x）和r（x）可以求得解析解，但是最好能计算出各种情况下的数值解。
为什么呢？
因为为了搞清楚y（x），你需要画出曲线，得到解析解时你也需要花费一定的精力在作图上。
那么应该怎么做？
首先将上面的二阶系统降低为两个一阶系统：
上式的常微分方程能够仿照（11）和（12）进行时间步进法数值积分。
举个简单例子：
（14）
降阶产生了两个一阶常微分方程：
（15）
假设初始条件为u1=1，u2=0，可以建立零维空间系统来积分这一对常微分方程。
打开COMSOLMultiphysics，等待模型导航栏窗口。
如果已经打开一个COMSOLMultiphysics窗口，把工作空间保存为mph文件或者把命令保存为m文件，然后在“File”菜单中选择“New”选项。
按照表6的步骤进行设置。
表6时间步进积分的简单实例
ModelNavigator
选择1D，COMSOLMultiphysics:
PDEmodes:
PDE,general模式，插入因变量：
u1u2
Draw菜单
Specifyobjects:
Line。
Coordinates弹出菜单，x：
01
名称：
interval，完成
Physics菜单：
Boundarysettings
选择域：
1和2（按住Ctrl键）
选择Neumann边界条件
保持默认选项q=0,g=0完成
Physics菜单：
Subdomainsettings
选择域：
1
设定F1＝－u2，F2＝u1
选择Init选项卡，设定u1（tw）=1
Mesh菜单：
Meshparameters
设置最大基元尺寸为1
点击Remesh，完成
Solve菜单：
Solverparameters
选择timedependent求解器，在General选项卡中输入
Times:
linspace（0,2*pi,50）
求解，完成
Post-processing：
Cross-sectionplot
parameters
Point选项卡，接受u1的默认设置
General选项卡，接受所有时间的默认设置
完成
这个例子给了我们两个独立变量u1、u2和空间坐标x。
注意“Subdomainsettings”窗口左上角方程中的矢量标志。
由于是矢量变量，所有输入选项卡都是矢量（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π）
升级会员 