第二章发展方程有限元分析

时间：2020-10-16 05:04:09 来源：蒲公英阅读网 本文已影响人

相关热词搜索：第二章 方程 有限元

　第二章

　发展方程的有限元分析 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:

　科学研究和工程应用中的偏微分方程（PDE）多源自复杂的平衡方程。常见的偏微分方程主要来自质量守恒、动量守恒、组分守恒和能量守恒定律。由于这些守恒定律是整个域上的积分方程，所以在连续性假设下，偏微分方程很容易用有限元方法近似描述。本章介绍了 COMSOL Multiphysics 中典型的三种不同类型“时间－空间”系统偏微分方程——椭圆方程，抛物线方程和双曲线方程。本章还对有限元方法进行了总体介绍，结合应用实例讲解有限元方法精确计算的特性，更深层次的内容将在后续章节中引出。

　1. 简介 在科学研究和工程应用中常会遇到满足守恒定律约束的偏微分方程，通常以积分形式出现。所有由质量守恒、动量守恒、组分守恒和能量守恒控制的传递现象都会产生连续逼近的偏微分方程。相信化工人员对传热、传质和动量传递现象不会感到陌生。

　与前一章 COMSOL Multiphysics 化工实例中介绍的零维、一维空间系统相比，化工课程中通常不会出现超过二维或三维的偏微分方程计算。从文献[1]中找到一个非常珍贵的例子。实际上，很多常见的化工模型和公式都是由实际过程中更高维数的动力学过程简化而来的。流体动力学中的阻力系数，传热传质系数，多相催化的 Thiele 模型，精馏塔设计中的 McCabe-Thiele 图等许多描述高维数系统传递现象或非稳态动力学过程的技术，都是半经验性的方法，也许可以用偏微分方程来描述这些过程，但是由于基本物理、化学现象的复杂性，这些方程通常很难求解。所以对于初步的设计计算，这些快速计算的简化方法很受欢迎，但是对于细节设计、设计翻新、过程分析和优化过程，只有简化方法是不够的。在基础科学研究过程中，这些方法仍然不断地从化学工程师传给生物学家、材料学家等，逐渐应用到各个领域的工作中。但是计算流体动力学（CFD）的出现彻底改变了这些方法在传导模型中的地位。虽然 CFD 在传导现象可视化和量化方面具有独特优势，但是唯象方法对于描述分布系统模型仍然有非常重要的作用。

　COMSOL Multiphysics 不是一个“商业 CFD 软件”，但是也可以做一些 CFD计算。它包含一些通用的 CFD 软件包，在某些模型支持上具有其独特的优点。对于 CFD 方法，大多数过程工程师希望有对湍流和燃烧模型的支持。但是COMSOL Multiphysics 不同，擅长多物理场计算。除了 CFD 的传统传递现象，COMSOL Multiphysics 还包含了电动力学、磁动力学、结构力学等实例模式，可以对这些现象进行模拟计算。实际上 COMSOL Multiphysics 的最大优点在于——首先，用户自定义编程容易，可以轻易建立用户自定义的模型，通过改变变量系数、边界条件、初始条件，还可以同时耦合多个物理场（甚至不在同一个域上的物理场）；其次，基于 MATLAB（或 COMSOL Script），可以实现对复杂模型模拟的所有编程功能，可以把 COMSOL Multiphysics 看作是一个方便的高阶有限元编程和分析软件。上一章中，我们介绍了一些用户自定义编程分析的功能，本章将介绍 COMSOL Multiphysics 在更高维数偏微分方程系统有限元建模上的核心

　优势。多物理场、扩展多物理场和无 PDE 约束的处理将留给后续章节。

　2 ． 偏微分方程 偏微分方程可以根据阶数、边界条件类型和线性度（线性、非线性和准线性）进行分类。令人惊讶的是，大多数科学研究和工程应用中碰到的偏微分方程都是二阶的，即最高导数项是二阶偏导数。这是偶然的么？这一巧合往往是在各自领域分别解释的。但是，最近 Frieden[2]证明所有已知物理定律都可以由最小 Fisher信息原理推导出，而该原律通常将一个二阶项作为最高项引入相应的物理定律中——包括从 Schrodinger 波函数方程到经典电动力学等各种物理定律。所以二维和三维的二阶“空间－时间”系统求解和分类在科学研究和工程中有着广泛的应用，这非常重要！出于这个原因，且有限元方法本质上非常适合处理二阶系统，所以有限元是有着广泛应用的技术。

　本章我们将集中介绍二维和三维空间的二阶系统。这里有三个典型的例子，在课本中经常见到，它们是：

　2 22 2222 22 2Laplace ( ): 0 ( ):

　 ( ):

　 u ux yu ut xu ut x       方程 椭圆

　扩散方程 抛物线

　波方程 双曲线 下面用 COMSOL Multiphysics 对这几种偏微分方程进行建模，同时在每个例子中我们将引入一个新的或特殊的 COMSOL Multiphysics 功能。这些椭圆、抛物线和双曲线方程可以归结为一种通用的线性二阶偏微分方程，具有两个独立变量和一个非独立变量：

　2 2 22 22 0u u u u ua b c d e fu gx x y y x y               

　 (1) 方程中的系数只是独立变量 x 和 y 的函数，或者是常数。根据以下原则分为三种经典形式：

　222

　: 0

　: 0: 0b acb acb ac   椭圆抛物线双曲线 这种分类方式大致反应了域中的信息流。例如，对于椭圆方程，边界点的信息会快速地传到所有内部点。所以椭圆方程是“非本地的”，意味着很远处的信息会对给定位置有影响，而“本地的”则表示只有附近信息才对其有影响。抛物线系统“扩散”信息，即信息会向所有方向传播。双曲系统“传播”信息，也就是说已经接收到信息的区域和将要接受到信息的区域有一个明显的划分。如果系统是线性的或者准线性的（即一些系数依赖于非独立变量，或者乘以一个更低阶的偏导数），这种分类方式和信息传播方式可以作为对二阶系统的一个引导。但是在非线性系统中，非线性破坏了信息传递结构，信息可能发生“异常”，也就是说，不再传递、超出给定范围、来自噪声或者在给定窗口内因丢失初始条件而消失。

　2.1 泊松方程：

　椭圆偏微分方程

　拉普拉斯方程经过适当变化就成为泊松方程：

　2( ) u f x  

　(2) 我们在第一章描述非均匀介质中具有分布热源的热传导问题时见过该方程的 1一维形式，但是这里的热传导率是均匀的。下面给出(2)式的另一个侧面，应该注意到，该方程满足给定涡量曲线的流函数方程：

　2( ) x     

　(3) 有两种较为容易的常见涡量类型——兰金涡（某一区域内的涡量是常数）和点源涡（涡量快速降低，可以理想为一个点涡）。你可能会对这两种类型涡产生的流线感兴趣。

　下面我们来研究一下这些流线。打开 COMSOL Multiphysics 模型导航栏。按照表 1 的步骤求解具有单位分布源项（常涡量）的泊松方程。该模型具有一个独立变量 u，二维空间坐标 x 和 y。我们绘制一个实心圆作为研究域，COMSOL Multiphysics 将边界分为四部分，默认网格划分 762 个。等流函数线如图 1 所示。

　表 1 单位分布源项的泊松方程 Model Navigator 选择 2D 空间维数，COMSOL Multiphysics| PDE Modes| Classical PDEs| Poisson’s Equation 设定 Element：Lagrange-Quadratic，完成 Draw 菜单 Specify objects：Circle 默认设置（单位半径；圆形在原点），完成 Physics 菜单：

　Boundary settings Dirichlet 默认模式，h=1，r=0(u=0) 完成 Physics 菜单：

　Subdomain settings 选择域 1 接受默认选项，c＝1；f＝1。完成 Meshing 点击三角形符号绘制网格 Solver 点击求解按钮（＝）求解 Post-processing：

　Plot parameters 在 General 选项卡中取消选中 Surface 选项 Contour 选项卡：选中 contour 框，完成

　 图 1 光滑圆形域中的常涡量流线

　流线用等高线的形式显示，图 1 通过“Menu”菜单中“Export：Image”选项导出。

　显然流线是一组同心圆。整个边界也是条流线（  ＝0），最大流函数出现在圆心（  ＝0.250），所以由常涡量引入的体积流量为 0.25。不改变几何形状，将网格数加密到 3045 个单元。由于流函数等高线相减的差值就是流线间的体积流率，流线的切线方向就是流动方向，所以图 1 的含义就很清楚了——常涡量导致圆形域中的环形流动。涡量从边界扩散到每个地方，导致边界上产生最大流速，圆心速度最小。

　下面来看域中的点源情况。我们将在域中画一个非常小的圆形作为第二个域，来近似点源情况。在第二个域中设定 f=1/πR 2 ，这里 R 是小圆的半经，小圆外侧 f＝0。这样 f 积分后为单位 1，当 R→0 时，f 趋近于 Dirac delta 函数。但是极限将通过细化网格和使用弱形式来逼近。

　打开“Draw”下拉菜单，在域中定义一个点。在域中明确绘制一个容易区分的点，会出现两种情况——可能（但不是必须）会将点贡献或点约束施加在方程系统上，也可能点本身成为产生网格的节点，约束网格。在网格节点容易（可能计算的更准确）求解独立变量的值，因为不需要进行值的网格内插。在该点附近可能需要更精细（或粗糙）的网格划分。按照表 2 的步骤来实现。

　表 2 位于圆形域圆点处的点涡量产生的流线 Draw 菜单 Specify objects：Point

　设置点 x=0, y=0，完成 Physics：

　Point settings Point selection：设置为 3（原点）

　在 week term 项输入 u_test 完成 Physics 菜单：

　Subdomain settings 选择域 1 设置 c=1；f＝1 完成 Meshing 点击三角形按钮绘制网格 Solve 菜单：

　Solver parameters Advanced 选项卡：改变解的形式为 weak 完成 Post-processing 点击原点。报告栏显示：Value：0.805017，表达式：u， 位置：(0,0)

　 图 2 位于圆形域圆点处的点涡量产生的流线

　 图 3 位于圆形域圆点处的点涡量产生的流线

　图 2 给出了表 2 中定义的弱的点源项产生的流线。从域中流函数的值能够发现，体积流率要比均匀分布源项情况下大很多。从本质上讲这是一组围绕圆点的等高线，不太满足圆心对称，当然在趋近圆心处同心圆消失了。读者可能非常想知道“u_test”是什么。如何能够将“弱的点条件”转化为点源。关于弱条件和弱约束的讨论我们将留在本章后面，等我们更加详细的介绍完有限元方法以后再加以介绍。但是这里提一下，它需要满足有限元方法的一个强项，就是不连续函数，如 Dirac delta 函数、步进函数的半分析处理。这些问题起源于工程动力学，例如合成材料或普通材料的性能在表面会发生变化。

　将网格细化到 2986 个单元，并不能使圆形等高线看起来更为圆滑，但是最大流函数值从 0.805 增大到 0.915。这是由细化圆心位置网格导致的。打开“Mesh”菜单，选择“Parameters”选项。在“Point”选项卡中输入顶点 3 的最大基元尺寸为 0.001。重新划分网格，会产生 1658 个基元，但是原点处的网格更密，计算得到圆点处  ＝1.57。

　图 2 中最大流函数为 1.57。将网格增加到 6632，最大流函数为 1.68。类似增大网格求解，26528 个网格对应流函数 1.79，106112 个网格对应流函数 1.90。虽然不知道网格是否能够收敛，但是由于涡量随着离源点距离的增加快速降低，本质上流线分布已经收敛。

　习 练习 2.1 求解当涡量随半经指数降低时的流线，即：

　2 20 exp() x y     

　2.2 扩散方程：

　抛物线偏微分方程 一维非稳态扩散方程的形式能够很好地适用于以下三种常见的传递性现象：

　222222

　pc cDt xT k Tt c xt x        质量扩散：

　 热传导：二维涡量输运：

　这里 c，T 和 ω 分别是浓度、温度和二维流动中 z 方向的涡量，它们相应的扩散系数为 D，α 和 ν。该方程已经在本科课程中详细介绍过。它可以通过傅立叶和拉普拉斯变换求解，初始和边界条件下的相似解是否失效取决于变量Dtx  。这里没有有限元方法求解的机会——只有求解这类老问题的老办法，对吗？不对。COMSOL Multiphysics 仍然可以改进、计算这类问题，这一点常被忽略。COMSOL Multiphysics 非常适用于非常数系数情况，即传递性能取决于场变量。例如，对于低压和高温情况，气体必须满足理想气体定律：

　nM PMV RT  ＝

　(4) 这里 R 是气体常数，M 是组分相对原子质量。在这种情况下，很少有气体热容能够保持常数。例如，在一定温度范围内，CO 2 气体的热容可以用温度的二次函数近似(MJ/kg mol ℃)：

　5 236.11 0.04233 2.887 10pc T T   

　 (5) 需要满足：

　3 7 2( )36.11( 273)( )1 1.172 10 7.995 10pk kRf Tc PMTf TT T    

　(6) 对于时间和位置满足：

　236.11kR xt XPML L  

　(7) 带入热传导方程得到简化形式为：

　22( )T Tf TX   

　(8) 下面来考虑在 CO 2 气体层中的热传导，设气体层长度为 L，水平边界保持恒温 400℃和 500℃。

　表 3 非线性介质中的热量传递——扩散率随温度改变。稳态分析 Model Navigator 选择 1D 空间维数，COMSOL Multiphysics| PDE Modes|

　Classical PDEs| Heat Equation 对于因变量和基元类型保持默认设置，完成 Draw 菜单 Specify objects：Line。设定 x：0 1 名称：interval，完成 Physics 菜单：

　Boundary settings 选择域 1 选择 Dirichlet 条件，设定 h＝1；r＝500(T=500) 选择域 2 选择 Dirichlet 条件，设定 h＝1；r＝400(T=400) 完成 Options 菜单：

　Constants 名称

　表达式 a1

　1.172E-3 a2

　7.995E-7 F400

　421.5

　完成 Options 菜单：

　Scalar Expressions 名称

　表达式 FofT

　(u+273)/(1+a1*u+a2*u 2 /F400) 完成 Physics 菜单：

　Subdomain settings 选择域 1 设定 f=0; c=1; ad =1/FofT 选择 Init 选项卡；设定 u(tw)=500-100*x，完成 Meshing 点击工具栏上三角形按钮绘制网格 点击两次“Refine mesh”按钮 Solver 电解求解按钮（＝）计算 （稳态线性求解器）

　Post-processing 点击 0.5 值：450，表达式：u，位置：(0.5)

　打开 COMSOL Multiphysics 模型导航栏，按照表 3 的步骤建立模型。该模型给出一个独立变量 u 和一维空间坐标 x。点击工具栏上网格划分按钮，然后再次加密网格两次，使网格单元达到 60 个。点击“＝”按钮进行求解（稳态线性求解器）。稳态解不随初始条件发生变化，因为经过足够长的时间后，累计项逐渐消失。当然在稳态求解器下是这样。所以我们可以用依赖于温度的扩散系数来研究对瞬态解的影响。所以改变初始条件为 u(t 0 )＝400℃，即 t＝0 时左侧边界条件跃升至 u＝500℃。打开“Solver”菜单选择“Parameters”选项，会弹出求解器变量设置对话框。按表 4 内容进行瞬态分析设置。经过后处理生成图 4 和图 5。图 4 给出了不断逼近稳态的温度曲线（几乎线性），图 5 给出了整个区域中间点的温度变化。

　表 4 非线性介质中的热量传递——扩散系数随温度变化。瞬态分析 Physics 菜单：

　Subdomain settings Init 选项卡，设定 u(t 0 )=400 完成 Solve 菜单：

　Solver parameters General 选项卡，选择 time-dependent solver 设定 output times：0:0.01:0.2 Advanced 选项卡，设定解形式：General。完成

　求解（＝）

　Post-processing：

　Domain Plot 接受默认设定，完成 Post-processing：

　Cross-section plot parameters Point 选项卡，设定 x：0.5，完成

　 图 4 从 t ＝0 到 t ＝0.2 的温度曲线

　 图 5

　x ＝0.5 处的温度变化

　图 4 和图 5 给出了向前扩散的速率，问题是非线性介质在其中起什么作用。特别是由于扩散率随温度升高而增大，我们发现与常数扩散率相比，曲线能够更快达到稳态。由于表 4 中没有出现自相似项Dtx  ，所以高温比低温更快达到稳态线性曲线。图 5 给出了域中心点处温度升至稳态值的过程，曲线呈 S 形，比预计要快。

　习 练习 2.2 设常数扩散率 f(T)=1 求解图 4 和图 5，这些曲线有什么区别？ 2.3 波动方程：双曲偏微分方程 一维正则波动方程已经研究的非常透彻，尤其在化工实例中经常见到。最典型的就是化工课程中容易忽视的声波。这是否意味着化学和过程工业中波动方程不重要呢？实际上总是会遇到很多复杂的波动方程或者非线性发展方程。例如，前沿研究经常关注的化学波反应器，冷凝器表面的波，旋转喷油嘴，分裂蒸馏塔质量传递影响，声学，超声波和声化学等等。但是碰巧化学工程师很少学习波的知识，市面上也很难找到经典的课本，分析声学在化工单元中的应用。

　我们将一维正则波动方程留作练习，这里以 Korteweg-de Vries 波传播方程为例加以介绍：

　336 0u u uut x x      

　(9) 这是薄层表面波的波方程——常见的用石头打水飘就是这样一个过程——在物

　理学中广泛存在，因为它可以在任何波传播介质中产生，不会消散，从某种意义上讲是非常“微弱的”。作者最近在地球物理[3]和表面涂层应用[4]方面对这类波方程加以推导。在这两方面工作中，都使用特殊的线性化方法(MOL)模拟作为数值分析的指导。线性化方法使用有限差分方法离散偏微分项，所以在网格点产生了一套可以同时求解 u 值的常微分方程。COMSOL Multiphysics 用有限元方法求解时间依赖的偏微分方程系统，本章后续章节还会遇到类似情况。与使用常微分方程计算网格点的 u 值不同，从 COMSOL Multiphysics 求解常微分方程中产生的“自由度”，更容易得到任何位置独立变量 u 的值。出于技术原因，我们不考虑用 FEMLAB 建模计算从[3][4]中推导出的类似与(9)的波方程——我们希望严格控制数值积分过程，使用刚性系统或者全隐式进行时间积分，因为这类方程很容易导致数值非稳定性。但是 COMSOL Multiphysics 具有强有效的时间积分求解器，会自动选择刚性求解器。所以仅仅出于好奇，我们决定用 COMSOL Multiphysics 求解孤波传递的 Korteweg-de Vries 方程。

　我们选择 COMSOL Multiphysics 而不用 FEMLAB 的一个重要原因在于对高阶偏导数项的处理。据我所知，FEMLAB 不能处理类似33xu的项。在 COMSOL Multiphysics 中也没有定义 u xxx ，但是定义了 u xxt ，使得我们可以通过一些办法来构造 u xxx 。我的办法是引入一个辅助变量 v，生成一个与(9)等价的系统：

　222[3 ] 0uu vt xuvx    

　(10) 第一个方程是通用守恒形式，第二个方程是一个关于 v 的椭圆方程。所以整个系统是一个“微分－几何”方程，不是所有的变量都需要“发展”。但是初始系统具有双曲特性。

　由于假设波传播没有方向性，理想情况下应当在无限介质中模拟。因为没有哪台计算机有无限的内存，所以提供足够大介质使得边界效应可以忽略的情况就显得更有意义，于是引入周期性边界条件。这里我们将根据方程(9)模拟一个孤波的传播。如果这些孤立或者本地的结构足够紧凑，域足够长，波就不会和它的周期性反射波相互作用。通过结果分析我们将看到这一点。方程(9)的孤波解析解为：

　2 2 3( , ) 2 sech ( 4 ) u x t a ax a t  

　这里 a 是振幅。注意到传播速率（系数 x 与 t 的比值）依赖于 a 2 ，因此波振幅越大传播速度越快。这是由 Korteweg-de Vries 方程的非线性引起的。有必要知道初始条件并进行模拟，用二阶偏微分项定义 v：

　4 4( , 0) 4 (cosh(2 ) 2)sech ( ) v x t a ax ax   

　打开 COMSOL Multiphysics 模型导航栏。表 5 中包含建立孤波传递模型的具体步骤。

　表 5 非线性波方程举例——Korteweg-de Vries 孤波传递

　Model Navigator 选择 2D 空间维数，COMSOL Multiphysics| PDE Modes| General Mode。设定两个因变量：u v 设定 Element：Lagrange-Cubic，完成 Draw 菜单 Specify objects: Line x:-10 10，名称：interval，完成 Options 菜单：

　Constants 定义常数 a＝1。我们希望通过调节 a 研究振幅的影响 Options 菜单：

　Scalar Expressions 定义以下变量：

　initialu2*a^2*(sech(a*x))^2 initialv4*a^4*(-2+cosh(2*a*x))*(sech(a*x))^4 Physics 菜单：

　Boundary settings 设定两个端点都为 Neumann 边界条件，G 1 ＝G 2 ＝0 Physics 菜单：

　Boundary conditions 选择边界 1，输入表达式：u。点击约束名框自动生成约束名 pconstr1。Destination 选项卡，选中边界 2，选中“Use selected boundaries as destination”选项，输入表达式 u。Source vertices 选项卡，选择边界 1，添加>>。Destination选项卡，选择边界 2，添加>>。返回 Source 选项卡，对于pconstr2，将 u 替换为 v，重复所有步骤。

　Physics 菜单：

　Subdomain settings 选择域 1 设定 Г 1 ＝3*u^2+v；Г 2 =ux；F 1 =0；F 2 =v；ad (v 2 entry)=0 Init 选项卡，U(t 0 )=initialu; v(t 0 )=initialv，完成 Mesh 菜单：

　Mesh parameters Global 选项卡，设定最大基元尺寸 0.02 重绘网格，完成 Solve 菜单：

　Solver parameters General 选项卡，设定 output times：0:0.01:3 Advanced 选项卡，设定约束保持 Lagrange Multiphiers。

　设置没有正规化空间函数。设置解的形式为 weak。完成 点击求解按钮（＝）

　Post-processing：

　Plot parameters 选择 line 选项卡，设定 u（默认）和 v 作为绘图变量。

　Animate 选项卡，开始 animation。

　 图 6 KdV 方程孤波传递解的初始波形

　图 7 KdV 方程孤波传递解拉普拉斯算子 v 的初始波形

　图 6 和图 7 分别给出了孤波解和方程(10)的拉普拉斯算子的初始波形。表 5

　建立了一个周期性边界条件。应当注意到，对于 Korteweg-de Vries 系统(9)，正振幅的波向右移动，负振幅的波向左移动。由于这里只有时间的一阶偏导数项，所以 Korteweg-de Vries 方程不可能使波向正则波动方程那样同时向两个方向传递。

　习 练习 2.3 正则波动方程，2222xutu，具有相应的经典偏微分方程应用模式。该应用模式具有二阶空间和时间偏导数。改变边界条件（Neumann 边界条件），对于正则波动方程建立表 5 中的周期性边界条件。尝试初始条件 u(t 0 )=sin(20*Pi*x)和u_t(t 0 )=-10*Pi*cos(10*Pi*x)。你预计 u(t)和 u_t(t)的结果如何？计算结果是否与你预料的不同？注意 COMSOL Multiphysics 有预定义的 pi。

　习 练习 2.4 当 a＝2 时，求解 Korteweg-de Vries 模型。理论上讲，振幅因子会导致传播速度增大 a 2 倍。模拟结果如何？对克服这个问题有什么建议？ 3 ．有限元方法 到目前为止，你一定很想知道 COMSOL Multiphysics 是如何求解偏微分方程系统的。有限元分析已经存在了几十年，从 20 世纪 80 年代开始就出现了相应的商业软件，Reddy[5]的书对此给出了很好的介绍。这里不是想详细介绍有限元方法，也不是仔细讲解 COMSOL Multiphysics 求解是如何实现的，而是希望对有限元方法计算类型有一个大致印象，理解为什么说 COMSOL Multiphysics 是实现有限元方法的一个很方便的工具。

　有限元方法的本质是将任何场变量约束以弱形式来表示。理解什么是弱形式（以及为什么数学上要求必须是弱形式），就必须要理解系统约束的强形式是偏微分方程系统和适当的边界条件。它为什么强？由于场变量需要连续并且要有与方程阶数相同的连续偏导数。这就是强约束。而弱形式对函数满足的约束条件限制要相对弱一些——不连续但是可积。

　来看一个偏微分方程和它的等价弱形式。考虑在三维空间域 Ω 中的具有一个独立变量 u 的稳态偏微分方程，通用形式为：

　( ) ( ) u F u  

　(11) 我们假设测试函数 v 是域 Ω 中定义的一个任意函数，并且满足函数 v ∈V。方程(11)两侧同乘以 v，并在整个域内积分有：

　( ) ( ) v u dx vF u dx   

　 (12) 这里 dx 表示体积微元。根据散度定律，我们得到：

　( ) ( ) v ndS v u dx vF u dx        

　(13) 当偏微分方程被 Neumann 边界条件约束时，方程(13)中的边界项消失。这也是有限元方法有 Neumann 自然边界条件的一个原因。回顾第一章，我们发现有限微分法有自然的 Dirichlet 边界条件。这导致体积积分变为：

　[ ( ) ( )] 0 vF u v u dx  

　 (14)

　这必须对所有的 v ∈V 都满足。现在来看有限元和基底函数。假设 u 可以分解为一系列基底函数：

　( ) ( )i iiu x u x   

　 (15) 例如，如果i 是正弦、余弦等基本渐进谐波函数，那么(15)就是由一系列傅立叶函数构成的。相反，在有限元方法中，要选择那些只在某一基元中作用的基底函数，即除了某一基元，在其它基元中都为零。

　 图 8 一维空间临近基元的两段线性基底函数

　图 8 给出了一维空间中两个拉格朗日线性基底函数的例子。显然，任意函数u(x)都可以由分段线性基底函数和足够小的基元进行任意精确度的近似。基底函数可以使用更高阶数，但是相应的每个基元也需要更多的未知变量 u i 。所以未知变量个数随着基底函数阶数的增大而增大，基底函数的个数也随着阶数的增大而增大。对于拉格朗日线性基底函数，在任何单元中都是由两个基底函数组成的：

　1       

　 (16) 这里  是基元中的本地坐标。所以对于 N 个单元，就有 2N 个基底函数i Φ i 。二次拉格朗日单元有三个基底函数：

　(1 )(1 2 ) 4 (1 ) (2 1)                

　 (17) 所以对于拉格朗日二次单元，N 个单元有 3N 个基底函数。

　回顾前面讲的，方程(14)必须对所有的 v∈V 满足，下面我们考虑所有由基底函数线性组合构成的函数空间，即：

　( ) ( )i iiv x v x   

　 (18) 但是由于 v 线性加入(14)，这就说明，如果以任意基底函数i Φ i 作为 v，方程(14)

　都成立，那么对于所有基底函数的线性组合(18)也成立，即所有 v ∈V 都成立。所以(18)等价于一个(k+1)N 个方程（每个i Φ i 都有方程(14)）、(k+1)N 个未知变量（u i ）的系统，这里 k 是基元的阶数（k＝1 线性，k＝2 二阶）。

　下面解释为什么 COMSOL Multiphysics 有限元方法具有这样的作用。COMSOL Multiphysics 建立(k+1)N 个方程(14)的集合。首先，注意到 Г(u)和 F(u)都是 u 的通用非线性函数。所以一般没有闭合解。在第一章中，我们介绍了COMSOL Multiphysics 中预置的零维问题非线性求解器，即 f(u)=0，u 是唯一未知变量。非线性求解器基于牛顿方法。牛顿方法的 N 维矢量方程为：

　( ) 0 L U 

　 (19) 这里 U 是未知数 u i 的矢量，L(U)是取代(14)中 v 的基底函数i Φ i 的方程系统，即：

　0 0 0( )( ) ( ) K U U U L U  

　 (20) 这里的 K(U 0 )是刚性矩阵刚度矩阵，L(U 0 )是荷载向量。刚性矩阵刚度矩阵是 L 的负雅克比矩阵：

　0 0( ) ( )LK U UU 

　 (21) 所以(20)是给定近似解 U 0 情况下 U 的线性方程。因子此，如果 U 0 能够足够接近真实解，线性方程(20)会找出更加精确的解 U，该过程可以不断重复，直到最后解足够精确。显然，牛顿法非线性求解器是 COMSOL Multiphysics 偏微分求解器的核心。COMSOL Multiphyiscs 可以实现偏微分方程有限元分析的所有步骤。如果可以的话，它用符号组成非线性运算符 L(U)的雅克比矩阵，如果不行的话，采用数值方法建立雅克比矩阵。如果偏微分方程本身是线性的，就不会很复杂。但是建立刚性矩阵刚度矩阵是一个非常庞大的工作——这在有限元分析中很常见。不久以前，有限元分析，包括网格划分和建立刚性矩阵刚度矩阵还经常是很多博士的研究课题。对于新的偏微分方程组合，甚至变系数问题（例如 1.2 节介绍的准线性系统），建立刚性矩阵刚度矩阵的工作量大的吓人。进一步讲，方程(20)的稀疏求解方法和时间积分需要通过编程工作来调整使其适合于某一问题。好在 COMSOL Multiphysics 将这个工作作为一个子程序（MATLAB 函数）来完成，实现与复杂偏微分方程系统的无缝衔接。

　在介绍边界条件实现方法之前，我们先通过一个常微分方程练习来强调一下刚才讨论的概念。通过该练习了解一下弱形式，找出给定二阶常微分方程的近似数值解。

　习 练习 2.5 有限元详细计算实例。使用有限元方法核心变分原理求解一个简单的常微分方程，详细介绍计算过程中出现的概念。问题比较简单，二阶常微分方程如下：

　2224 8d uu xdx 

　(22) 边界条件为：

　(0) ( ) 04u u 

　使用了弱方程式。

　该二阶常微分方程解析解如下（试证明之！）：

　22( ) cos(2 ) 1 sin(2 ) 2 18u x x x x        

　(23) 因为已经知道了解析解，就可以比较近似解的误差。

　建立弱方程式的第一步是假设一个权重函数和试探函数。将 U(x)作为试探函数， ( ) x  Φ(x)作为权重函数。我们一会儿再讨论这些函数的准确形式。试探函数U(x)构成了常微分方程的解。所以如果将 U(x)带入(22)式，可以得到残差为：

　2224 8d UR U xdx  

　 (24) 接下来的工作就是最小化残差。要最小化残差首先要计算权重残差，将方程(24)两侧同乘以 ( ) x  Φ(x)，然后在整个域中进行积分（即 0≤x≤π/4）

　2/4220( ) 4 8d UR x U x dxdx       

　 (25) 使用分步积分法简化上式为：

　/4 /420 0( ) 4 8d dUR x U dx x dxdx dx        

　 (26) 为后续计算方便，我们需要做一些关键假设。因为只要满足边界条件，我们可以任意设定 U(x)和 Φ(x) ( ) x  ，我们令 U(x)＝ ( ) x  Φ(x)。这被称作 Galerkin 方法。如果 U(x)≠ ( ) x  Φ(x)，给出 Rayleigh-Ritz 公式。我们需要选择一个 x 的代数函数以满足边界条件 u(0)=u(π/4)=0。

　1 1 2 20( ) ( ) ( )NN N i iix U x x u u u u             

　 (27) 假设 N 个方程为：

　21 24 4 4NNx x x x x x                        

　(28) 所以 1( )4Niiix u x x    

　(29) 这种选择满足边界条件，且不需要考虑方程包含的项数。由于 U(x)＝Φ(x) ( ) x  ，权重残差为：

　2/4 /42 20 01( ) 4 82dR x dx x dxdx             

　 (30) 将(29)式带入(30)并积分，得到一个与 x 无关的 R 表达式。但是该表达式包含 N个未知变量 u i 。我们需要计算 u i 的值使得权重残差最小。

　由(29)式得：

　1( 1)4i iidu ix i xdx      

　 (31) 所以我们有 3141 2 1( ) 24 1 2 31 2 ( 1)( 1)2 4 1 11 18

　, 1,2,3, ,4 3 4i ji i ji jj jiiR u uui j i j i jij ij i j i ju ui j i j i ju i j Ni i                                             

　(32) 然后通过将 R 对 u j 取导数，最小化残差。对于预定义的 N，产生了 N 个需要同时求解的代数方程。

　( )0ijdR udu

　 (33) 当 N＝2 时只有两个未知变量，u 1 和 u 2 ，产生两个线性方程：

　1 21 20.122 0.048 0.1200.048 0.033 0.063u uu u    

　(34) 矩阵形式表示为：

　120.122 0.048 0.1200.048 0.033 0.063uu               

　(35) 通用形式为：

　[ ]{ } { } K u L 

　(36) 这里 K 是雅克比矩阵（刚性矩阵刚度矩阵），u 是未知向量，L 是限制向量（负载向量）。方程(35)的解为：

　120.5545791.112560uu   所以方程(22)的解为：

　2( ) 0.554579 1.112564 4u x x x x x               

　 (37) 如果进一步假设 N＝3，得到含三个未知变量 u 1 ，u 2 和 u 3 的代数方程，矩阵形式带 格 式 的: 字体: 加粗, 非倾斜带 格 式 的: 字体: 加粗, 非倾斜带 格 式 的: 字体: 非倾斜

　为：

　1230.1216 0.0477 0.0228 0.1200.0477 0.0328 0.0200 0.6300.0228 0.0200 0.0139 0.351uuu                       

　 (38) 方程(38)的解为：

　1 2 30.588 0.838 0.349 u u u      

　所以方程(22)的解变为：

　2 3( ) 0.588 0.838 0.3494 4 4u x x x x x x x                        

　 (39) 图 9 给出了(37)(39)式和解析解的比较。可以清楚看到，只要取的项数足够多，近似代数解能够与精确解很好的吻合。

　 图 9 (39)式曲线。三项几何方程能够与解析解很好的符合

　通过本例，你应该对解的弱形式有了清晰的概念。下一节我们将讨论如何在COMSOL Multiphysics 中实现边界条件。

　3.1 边界条件 正如前面例子中提到的，读者应当注意到刚性矩阵刚度矩阵 K 等价于Neumann 边界条件。在第一章中学过，纯 Neumann 条件将会导致刚性矩阵刚度

　矩阵奇异，这样 COMSOL Multiphysics 就不能直接求解，因为它会对用投影法求出的刚性矩阵刚度矩阵特征系统的解加上任意常数。有限元方法的一个特别之处就在于边界条件的处理。

　我们打算像有限微分法那样处理边界条件。矩阵方程中的某些行将会被对未知变量 u i 的直接约束所取代，以保证矩阵的阶数。但是这会破坏刚性矩阵刚度矩阵的稀疏性，对一些边界条件有不利影响，所以需要人为的全矩阵求解器。对于相同矩阵方程，全矩阵求解器与稀疏求解器相比计算精度和效率都会有所降低。通过拉格朗日因子，有限元方法可以得到一个很好的解，这是处理等式约束优化问题非常著名的方法。假设除了偏微分方程约束，我们还有一系列边界条件，对所有的 v ∈V 都满足弱形式。通过将每个基底函数展开并写出边界积分，类比偏微分方程约束(19)，可以得到边界约束的矢量方程：

　( ) 0 M U 

　(40) 该约束的剩余方程并不需要 N 个方程。通常它只是其中包含 N 个未知变量的少数几个方程，因为不是所有的基底向量都会对边界约束起到测试函数 v 的作用。方程(40)的线性形式类似与(20)：

　0 0 0( )( ) ( ) N U U U M U  

　(41) 这里 N 是 M 的负雅克比矩阵：

　0 0( ) ( )MN U UU 

　(42) 下面需要一定的处理技巧。刚性方程乘以未知向量 Λ，Λ 是拉格朗日因子乘以N T ，这里 T 表示转置：

　0 0 0 0( )( ) ( ) ( )TK U U U N U L U    

　 (43) 这里有什么技巧呢？如果约束(40)成立，将会有唯一的拉格朗日因子满足(43)式。但是(41)和(43)式不只适用与边界条件，任何约束——内部逐点，域积分，边或边界等所有能够以弱形式表示的约束都可以用拉格朗日因子法处理。

　拉格朗日因子

　那么拉格朗日因子如何满足 M(U)=0 呢？根据变分原理。当 Λ＝0 时，(43)式等价于对下式应用最小化原理：

　1min2T TUU KU U L   

　(44) 如果想确认约束(40)是否也同时满足，我们对(44)额外增加一个权重项，确保不满足 M(U)=0。权重项就称作拉格朗日因子 Λ。

　1min2T TUU KU U L M    

　 (45) 下面用线性化 M(U)=M(U 0 )+N(U 0 )(U-U 0 )来简化(45)。注意到常数项对最小化不起作用。

　1min2T TUU KU U L NU    

　 (46)

　最小化(46)等价于(43)，即(43)式的解能够最小化(46)。在有限元方法中，(46)式相当于“能量最小化”，即用刚性矩阵刚度矩阵 K 作为权重。注意不要与最小二乘法的最小化混淆，它与(45)式类似：

　2 1min ( )2UL U M   

　 (47) 线性化 L 和 M，并经过整理并忽略常数项为：

　01min ( )2T T T T T TUU K K U U U K L U N      

　(48) 所以根据线性代数理论，(48)式的解 U 就是正则方程[6]的解：

　10( ) ( )T T T T TK K U U K K N K L   

　 (49) 也是下式最小二乘法的解 10( ) ( )T TK U U K N L   

　(50) 所以从最小二乘误差意义上讲，方程(50)的解满足约束(40)，M(U)=0；同时从最小能量意义上讲，约束(43)也满足(40)。为使求解简单、计算精确，COMSOL Multiphysics 用(43)代替(50)。两者的差别非常明显，最小二乘误差方程(47)适用于通用非线性算子 L(U)，而刚性能量最小化方程(45)只适用于对称正定矩阵 K。否则，(43)式中的拉格朗日算子就只是出于方便，并不能保证约束(40)满足近似条件。(50)是强更强的条件，但是是以额外矩阵操作为代价的。(47)表明 M(U)不是惩罚性的约束，所以更强的条件必须明确将每个约束分别看作是一个惩罚[7] 22 1min ( ) ( )2j jUjL U M U    

　 (51) 该方法没有将计算结果简介表达为矩阵方程的形式，因为约束项含有 N，M，K和 Λ。

　弱条件

　下面来看如何处理 2.1 节中涡量的点源。它满足弱条件，只计算如下积分， 0( )xv x dx v 

　 (52) 并对刚性矩阵刚度矩阵和(44)中的荷载向量有适当贡献。

　3.2 基元 有限元方法的一个基本思想是将每个域看作由很多个各种形状的更小子域组成。这些子域就被称作有限元。基元的角点被称作节点，场变量在节点上进行计算。在角点之间也可以有节点，称作边节点。在 COMSOL Mulitiphyiscs 中划分网格时，将整个计算域分解为一个个可选基元并且构建相应的节点。在实际应用中可以找到一百多种类型的基元。作为初学者，可能会对使用哪种类型基元和基元的个数感到困惑。

　离散化过程设定了基元的类型和个数。基元的个数与解的准确性直接相关。基元个数越多，误差越小。但是大的基元个数也会增大计算成本，需要更大的内

　存空间和更长的计算时间。

　经常见到适当往多的设定基元个数，因为没有优化基元个数的公式。一个域中选择的基元个数完全靠经验来决定。虽然准确性随着基元个数 N 的增大而增大，但是存在某个数 N c ，当 N 大于 N c 时，精确度敏感性可以忽略。

　 图 10 归一化误差随基元个数变化曲线。随着 N 的增大准确性提高，但是当超过某特定值 N c 时，改进效果可以忽略。

　图 10 给出了归一化误差随基元个数 N 变化的曲线。在每次迭代时基元个数翻倍。可以看到，最后三个点在计算精度上并没有明显的改进。可以通过适当增加一些运算来得到使用的基元个数。在很多情况下，人们对某个区域更感兴趣，而不是整个计算空间。例如考虑液体流过圆筒的过程，希望研究圆筒壁面的边界层现象。在这种情况下就可以、并且应该增大圆筒壁面处的基元密度，降低远离壁面处的基元密度。这样就可以在不增大基元个数的情况下提高准确度。COMSOL Multiphysics 允许使用这种网格拉伸的方法。在这点上还需要一提的是，拉伸网格中基元发生偏斜。偏斜基元使得雅克比矩阵不可用。所以在使用拉伸网格时必须非常小心。基元的类型取决于需要求解的问题。而域的维数决定了基元的维数。最简单的一维基元就是在两个端点间划分的一段直线。最基础的二维基元是三角形基元，在三维中就是四面体基元。表 6 给出了实际应用中对应维数下的一些基本基元。虽然我们在节点间一般使用直线段，但是曲线段能够提供更为通用的基元形式。

　表 6 基本基元 维数 形状

　1-D

　2-D

　三角形 矩形 菱形 3-D

　四面体 正六面体 不规则六面体

　虽然基元的几何形状是重要的区分方式，但是基元常根据对应的插值多项式进行分类。根据这种分类方法可以将其分为三种类型的基元：单纯型，复合型和多元型。如果多项式在边角节点处只有线性项和常数项，就称为单纯型基元。而符合型基元在角节点、边节点和内部节点处使用更高阶的多项式（平方、立方、五次方等等）。而多元型基元使其边界与坐标轴平行，并使用更高阶的多项式。

　正如前面提到的，复合型基元使用更高阶的多项式。而多项式的组成和节点的结构则由所使用的 Pascal 三角形，Pascal 四面体和 Pascal 超立方体来决定。表7 给出了 0－5 阶多项式的 Pascal 三角形。而采用的多项式必须要完整，即包含从最低阶到最高阶的每一项。例如 c 1 +c 2 x+c 3 x 2 是完整的，而 c 1 +c 2 x 2 就不完整，因为缺少一次项。对于任意二维基元，只要包含了表 7 中对应阶数的所有项，就可以获得所需的完整多项式。通常对于多项式的每一项都应该有一个节点。例如立方形基元就应该在每个边上设定四个节点，但是也有更复杂的结合情况。表 8给出了二维三角形的线性、平方和立方型基元。

　表 7 二维基元的 Pascal 三角形 1 常数 x

　y 线性 x 2

　xy

　y 2

　平方 x 3

　x 2 y

　xy 2

　 y 3

　立方 x 4

　 xy 3

　x 2 y 2

　 x 3 y

　y 4

　四次方 x 5

　x 4 y

　x 3 y 2

　x 2 y 3

　 xy 4

　y 5

　五次方

　表 8 1-3 阶基底函数复杂基元的类型

　 线性 平方 立方

　COMSOL Multiphysics 提供了两种忽略几何形状的主要基元类型，称作拉格朗日基元和厄密基元。这种命名方式主要来自它们所使用的差值多项式。拉格朗日基元使用拉格朗日多项式作为基底函数。假设在一维基元中，使用拉格朗日多项式 L n 表示场变量 u(x)，则有：

　1 1 2 2( ) ( ) ( ) ( )N Nx L x u L x u L x u     ｕ

　(53) 这里 u n 是未知系数。L n (x)为：

　1,1 2 1 11 2 1 1( )( ) ( )( ) ( )( )( ) ( )( ) ( )nMnM M NN MN N nN N N N N N N nx xLx xx x x x x x x x x xx x x x x x x x x x           

　(54) 展开式给出了对应阶数的多项式。拉格朗日基元通常出现在计算流体力学中。它们可以提供节点处的变量值。

　厄密基元使用厄密多项式插值场变量。拉格朗日基元和厄密基元的主要不同点在于有效自由度（DOF）。拉格朗日基元的自由度是节点处的函数值（由节点处的所有变量构成）。而厄密基元除了节点处的函数值，也考虑了角点处变量的一阶导数。由于第 n 个和 n＋1 个节点的值 u(x n )和 u(x n+1 )，导数i/ | du dx 和1/ | i du dx未知，所以需要用四个未知数来近似 u(x)。

　2 31 2 3 3( ) u x x x x        

　 (55) 由于 x i 和 x i ＋ 1 为已知坐标，所以可以轻易写出四个方程：两个节点的函数值和两个节点的一阶导数。

　2 31222 31 31 121 41 110 1 2 310 1 2 3ii i iii iii i iii iu x x xdu x xu x x xdu x x                                 

　(56) 这里 du 表示节点处 u 对 x 的导数。通过矩阵求逆，α n 可以用节点值和导数值表示，方程变为：

　1 2 1 3 1 4( ) ( ) ( ) ( ) ( )i i i iu x u x du x u x du x        

　(57) 这里 φ i (x)是厄密插值函数（该例中为立方函数，也称为三次样条）。

　如果你仔细观察会发现，立方插值函数在拉格朗日基元中需要四个节点而在厄密基元中只需要两个节点。厄密基元通常用在求解桁架的承载和受力情况。

　对于这两种类型，COMSOL Multiphysics 进一步提供了曲线网格基元和Argyris 基元。曲线网格基元主要用于精确近似真实边界条件。Argyris 基元是五阶的厄密基元，使用节点值和高达二阶的导数。它也在边的中点使用常规算子▽u。Argyris基元的五阶多项式需要确定21个常数。除了预定义的基元，COMSOL Multiphyiscs 也允许用户自己定义基元。感兴趣的读者可以参考 COMSOL Multiphysics 手册，手册中详细介绍了如何定义一个新的基元类型。

　这里我们对初学者如何理解和使用基本基元做了详细的介绍。使用COMSOL Multiphysics，我们不需要非常精通网格划分技术和基元设定。正如前面看到的，一旦确定了偏微分方程的求解域，只需要点击一个键就可以完成网格划分工作。但是如果读者对基本概念感兴趣，可以参考[8]及其引用文献，里面详细介绍了基元和网格划分技术的发展。

　习 练习 2.6 均匀介质中的热传递。为了介绍有限元方法如何用公式表达和具体的求解方法，下面来考虑一个例子。第一章 6.1 节介绍了一个非均匀介质中的传热问题。这里我们考虑一个更简单的问题。在整个域内认为传热系数为常数，f(x)不变，域的长度 L＝1。图 11 显示了该物理系统。我们假设穿过任意垂直与中心轴线的横截面的热流均匀（如图 11 所示），那么温度只沿轴线发生改变：实现了降低维数的作用。

　 图 11 绝热棒中的轴向传热。两端温度分别为 T s 和 T e ，棒长为 L ，横截面 A＝1。棒内恒定产热速率为 Q 。

　第一章中的公式(31)完整描述了该传热问题。由于横截面 A＝1，所以常热源Q 等于(1.31)式中的 f(x)。方程变为：

　0 0 1d dTk Q xdx dx      

　 (58) 求解过程用到的数据有：

　0213.3J/ Cms10J/sm11.25J/m sxxkQTq 

　第一步：变分公式

　该偏微分方程是圆柱中热传递方程的强形式。有限元方法的第一步就是转化为方程的弱形式。用权重函数乘以方程(58)并在整个域上积分。

　100d dTw k Q dxdx dx          

　(59) 分步积分（使用一维散度定律）得：

　11 10 00dw dT dTk dx wk wQdxdk dx dx          

　(60) 根据传热理论的傅立叶定律，穿过单位横截面的热流 q=-k(dT/dx)。所以， 1 1100 0[ ]dw dTk dx wq wQdxdk dx      

　 (61) 从前几节中我们知道，必须用基底函数多项式来近似未知数 w 和 T。选择这些多项式是有限元方法的第二步。

　第二步：离散和多...