2 发展方程的有限元分析.doc

第二章 发展方程的有限元分析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科学研究和工程应用中的偏微分方程（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 偏微分方程偏微分方程可以根据阶数、边界条件类型和线性度（线性、非线性和准线性）进行分类。令人惊讶的是，大多数科学研究和工程应用中碰到的偏微分方程都是二阶的，即最高导数项是二阶偏导数。这是偶然的么？这一巧合往往是在各自领域分别解释的。但是，最近Frieden2证明所有已知物理定律都可以由最小Fisher信息原理推导出，而该原律通常将一个二阶项作为最高项引入相应的物理定律中包括从Schrodinger波函数方程到经典电动力学等各种物理定律。所以二维和三维的二阶“空间时间”系统求解和分类在科学研究和工程中有着广泛的应用，这非常重要！出于这个原因，且有限元方法本质上非常适合处理二阶系统，所以有限元是有着广泛应用的技术。本章我们将集中介绍二维和三维空间的二阶系统。这里有三个典型的例子，在课本中经常见到，它们是：下面用COMSOL Multiphysics对这几种偏微分方程进行建模，同时在每个例子中我们将引入一个新的或特殊的COMSOL Multiphysics功能。这些椭圆、抛物线和双曲线方程可以归结为一种通用的线性二阶偏微分方程，具有两个独立变量和一个非独立变量： (1)方程中的系数只是独立变量x和y的函数，或者是常数。根据以下原则分为三种经典形式：这种分类方式大致反应了域中的信息流。例如，对于椭圆方程，边界点的信息会快速地传到所有内部点。所以椭圆方程是“非本地的”，意味着很远处的信息会对给定位置有影响，而“本地的”则表示只有附近信息才对其有影响。抛物线系统“扩散”信息，即信息会向所有方向传播。双曲系统“传播”信息，也就是说已经接收到信息的区域和将要接受到信息的区域有一个明显的划分。如果系统是线性的或者准线性的（即一些系数依赖于非独立变量，或者乘以一个更低阶的偏导数），这种分类方式和信息传播方式可以作为对二阶系统的一个引导。但是在非线性系统中，非线性破坏了信息传递结构，信息可能发生“异常”，也就是说，不再传递、超出给定范围、来自噪声或者在给定窗口内因丢失初始条件而消失。2.1 泊松方程：椭圆偏微分方程拉普拉斯方程经过适当变化就成为泊松方程： (2)我们在第一章描述非均匀介质中具有分布热源的热传导问题时见过该方程的1一维形式，但是这里的热传导率是均匀的。下面给出(2)式的另一个侧面，应该注意到，该方程满足给定涡量曲线的流函数方程： (3)有两种较为容易的常见涡量类型兰金涡（某一区域内的涡量是常数）和点源涡（涡量快速降低，可以理想为一个点涡）。你可能会对这两种类型涡产生的流线感兴趣。下面我们来研究一下这些流线。打开COMSOL Multiphysics模型导航栏。按照表1的步骤求解具有单位分布源项（常涡量）的泊松方程。该模型具有一个独立变量u，二维空间坐标x和y。我们绘制一个实心圆作为研究域，COMSOL Multiphysics将边界分为四部分，默认网格划分762个。等流函数线如图1所示。表1 单位分布源项的泊松方程Model Navigator选择2D空间维数，COMSOL Multiphysics| PDE Modes| Classical PDEs| Poissons Equation设定Element：Lagrange-Quadratic，完成Draw菜单Specify objects：Circle默认设置（单位半径；圆形在原点），完成Physics菜单：Boundary settingsDirichlet默认模式，h=1，r=0(u=0)完成Physics菜单：Subdomain settings选择域1接受默认选项，c1；f1。完成Meshing点击三角形符号绘制网格Solver点击求解按钮（）求解Post-processing：Plot parameters在General选项卡中取消选中Surface选项Contour选项卡：选中contour框，完成图1 光滑圆形域中的常涡量流线流线用等高线的形式显示，图1通过“Menu”菜单中“Export：Image”选项导出。显然流线是一组同心圆。整个边界也是条流线（0），最大流函数出现在圆心（0.250），所以由常涡量引入的体积流量为0.25。不改变几何形状，将网格数加密到3045个单元。由于流函数等高线相减的差值就是流线间的体积流率，流线的切线方向就是流动方向，所以图1的含义就很清楚了常涡量导致圆形域中的环形流动。涡量从边界扩散到每个地方，导致边界上产生最大流速，圆心速度最小。下面来看域中的点源情况。我们将在域中画一个非常小的圆形作为第二个域，来近似点源情况。在第二个域中设定f=1/R2，这里R是小圆的半经，小圆外侧f0。这样f积分后为单位1，当R0时，f趋近于Dirac delta函数。但是极限将通过细化网格和使用弱形式来逼近。打开“Draw”下拉菜单，在域中定义一个点。在域中明确绘制一个容易区分的点，会出现两种情况可能（但不是必须）会将点贡献或点约束施加在方程系统上，也可能点本身成为产生网格的节点，约束网格。在网格节点容易（可能计算的更准确）求解独立变量的值，因为不需要进行值的网格内插。在该点附近可能需要更精细（或粗糙）的网格划分。按照表2的步骤来实现。表2 位于圆形域圆点处的点涡量产生的流线Draw菜单Specify objects：Point设置点x=0, y=0，完成Physics：Point settingsPoint selection：设置为3（原点）在week term项输入u_test完成Physics菜单：Subdomain settings选择域1设置c=1；f1完成Meshing点击三角形按钮绘制网格Solve菜单：Solver parametersAdvanced选项卡：改变解的形式为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.2 扩散方程：抛物线偏微分方程一维非稳态扩散方程的形式能够很好地适用于以下三种常见的传递性现象：这里c，T和分别是浓度、温度和二维流动中z方向的涡量，它们相应的扩散系数为D，和。该方程已经在本科课程中详细介绍过。它可以通过傅立叶和拉普拉斯变换求解，初始和边界条件下的相似解是否失效取决于变量。这里没有有限元方法求解的机会只有求解这类老问题的老办法，对吗？不对。COMSOL Multiphysics仍然可以改进、计算这类问题，这一点常被忽略。COMSOL Multiphysics非常适用于非常数系数情况，即传递性能取决于场变量。例如，对于低压和高温情况，气体必须满足理想气体定律： (4)这里R是气体常数，M是组分相对原子质量。在这种情况下，很少有气体热容能够保持常数。例如，在一定温度范围内，CO2气体的热容可以用温度的二次函数近似(MJ/kg mol )： (5)需要满足： (6)对于时间和位置满足： (7)带入热传导方程得到简化形式为： (8)下面来考虑在CO2气体层中的热传导，设气体层长度为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条件，设定h1；r500(T=500)选择域2选择Dirichlet条件，设定h1；r400(T=400)完成Options菜单：Constants名称 表达式a1 1.172E-3a2 7.995E-7F400 421.5 完成Options菜单：Scalar Expressions名称 表达式FofT (u+273)/(1+a1*u+a2*u2/F400)完成Physics菜单：Subdomain settings选择域1设定f=0; c=1; =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(t0)400，即t0时左侧边界条件跃升至u500。打开“Solver”菜单选择“Parameters”选项，会弹出求解器变量设置对话框。按表4内容进行瞬态分析设置。经过后处理生成图4和图5。图4给出了不断逼近稳态的温度曲线（几乎线性），图5给出了整个区域中间点的温度变化。表4 非线性介质中的热量传递扩散系数随温度变化。瞬态分析Physics菜单：Subdomain settingsInit选项卡，设定u(t0)=400完成Solve菜单：Solver parametersGeneral选项卡，选择time-dependent solver设定output times：0:0.01:0.2Advanced选项卡，设定解形式：General。完成求解（）Post-processing：Domain Plot接受默认设定，完成Post-processing：Cross-section plotparametersPoint选项卡，设定x：0.5，完成图4 从t0到t0.2的温度曲线图5 x0.5处的温度变化图4和图5给出了向前扩散的速率，问题是非线性介质在其中起什么作用。特别是由于扩散率随温度升高而增大，我们发现与常数扩散率相比，曲线能够更快达到稳态。由于表4中没有出现自相似项，所以高温比低温更快达到稳态线性曲线。图5给出了域中心点处温度升至稳态值的过程，曲线呈S形，比预计要快。练习2.2设常数扩散率f(T)=1求解图4和图5，这些曲线有什么区别？2.3 波动方程：双曲偏微分方程一维正则波动方程已经研究的非常透彻，尤其在化工实例中经常见到。最典型的就是化工课程中容易忽视的声波。这是否意味着化学和过程工业中波动方程不重要呢？实际上总是会遇到很多复杂的波动方程或者非线性发展方程。例如，前沿研究经常关注的化学波反应器，冷凝器表面的波，旋转喷油嘴，分裂蒸馏塔质量传递影响，声学，超声波和声化学等等。但是碰巧化学工程师很少学习波的知识，市面上也很难找到经典的课本，分析声学在化工单元中的应用。我们将一维正则波动方程留作练习，这里以Korteweg-de Vries波传播方程为例加以介绍： (9)这是薄层表面波的波方程常见的用石头打水飘就是这样一个过程在物理学中广泛存在，因为它可以在任何波传播介质中产生，不会消散，从某种意义上讲是非常“微弱的”。作者最近在地球物理3和表面涂层应用4方面对这类波方程加以推导。在这两方面工作中，都使用特殊的线性化方法(MOL)模拟作为数值分析的指导。线性化方法使用有限差分方法离散偏微分项，所以在网格点产生了一套可以同时求解u值的常微分方程。COMSOL Multiphysics用有限元方法求解时间依赖的偏微分方程系统，本章后续章节还会遇到类似情况。与使用常微分方程计算网格点的u值不同，从COMSOL Multiphysics求解常微分方程中产生的“自由度”，更容易得到任何位置独立变量u的值。出于技术原因，我们不考虑用FEMLAB建模计算从34中推导出的类似与(9)的波方程我们希望严格控制数值积分过程，使用刚性系统或者全隐式进行时间积分，因为这类方程很容易导致数值非稳定性。但是COMSOL Multiphysics具有强有效的时间积分求解器，会自动选择刚性求解器。所以仅仅出于好奇，我们决定用COMSOL Multiphysics求解孤波传递的Korteweg-de Vries方程。我们选择COMSOL Multiphysics而不用FEMLAB的一个重要原因在于对高阶偏导数项的处理。据我所知，FEMLAB不能处理类似的项。在COMSOL Multiphysics中也没有定义uxxx，但是定义了uxxt，使得我们可以通过一些办法来构造uxxx。我的办法是引入一个辅助变量v，生成一个与(9)等价的系统： (10)第一个方程是通用守恒形式，第二个方程是一个关于v的椭圆方程。所以整个系统是一个“微分几何”方程，不是所有的变量都需要“发展”。但是初始系统具有双曲特性。由于假设波传播没有方向性，理想情况下应当在无限介质中模拟。因为没有哪台计算机有无限的内存，所以提供足够大介质使得边界效应可以忽略的情况就显得更有意义，于是引入周期性边界条件。这里我们将根据方程(9)模拟一个孤波的传播。如果这些孤立或者本地的结构足够紧凑，域足够长，波就不会和它的周期性反射波相互作用。通过结果分析我们将看到这一点。方程(9)的孤波解析解为：这里a是振幅。注意到传播速率（系数x与t的比值）依赖于a2，因此波振幅越大传播速度越快。这是由Korteweg-de Vries方程的非线性引起的。有必要知道初始条件并进行模拟，用二阶偏微分项定义v：打开COMSOL Multiphysics模型导航栏。表5中包含建立孤波传递模型的具体步骤。表5 非线性波方程举例Korteweg-de Vries孤波传递Model Navigator选择1D空间维数，COMSOL Multiphysics| PDE Modes| General Mode。设定两个因变量：u v设定Element：Lagrange-Cubic，完成Draw菜单Specify objects: Linex:-10 10，名称：interval，完成Options菜单：Constants定义常数a1。我们希望通过调节a研究振幅的影响Options菜单：Scalar Expressions定义以下变量：initialu=2*a2*(sech(a*x)2initialv=4*a4*(-2+cosh(2*a*x)*(sech(a*x)4Physics菜单：Boundary settings设定两个端点都为Neumann边界条件，G1G20Physics菜单：periodic conditions选择边界1，输入表达式：u。点击约束名框自动生成约束名pconstr1。Destination选项卡，选中边界2，选中“Use selected boundaries as destination”选项，输入表达式u。Source vertices选项卡，选择边界1，添加。Destination选项卡，选择边界2，添加。返回Source选项卡，对于pconstr2，将u替换为v，1和2调换，重复所有步骤。Physics菜单：Subdomain settings选择域1设定13*u2+v；2=ux；F1=0；F2=v；(v2 entry)=0Init选项卡，U(t0)=initialu; v(t0)=initialv，完成Mesh菜单：Mesh parametersGlobal选项卡，设定最大基元尺寸0.02重绘网格，完成Solve菜单：Solver parametersGeneral选项卡，设定output times：0:0.01:3Advanced选项卡，设定约束保持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正则波动方程，具有相应的经典偏微分方程应用模式。该应用模式具有二阶空间和时间偏导数。改变边界条件（Neumann边界条件），对于正则波动方程建立表5中的周期性边界条件。尝试初始条件u(t0)=sin(20*Pi*x)和u_t(t0)=-10*Pi*cos(10*Pi*x)。你预计u(t)和u_t(t)的结果如何？计算结果是否与你预料的不同？注意COMSOL Multiphysics有预定义的pi。练习2.4当a2时，求解Korteweg-de Vries模型。理论上讲，振幅因子会导致传播速度增大a2倍。模拟结果如何？对克服这个问题有什么建议？3有限元方法到目前为止，你一定很想知道COMSOL Multiphysics是如何求解偏微分方程系统的。有限元分析已经存在了几十年，从20世纪80年代开始就出现了相应的商业软件，Reddy5的书对此给出了很好的介绍。这里不是想详细介绍有限元方法，也不是仔细讲解COMSOL Multiphysics求解是如何实现的，而是希望对有限元方法计算类型有一个大致印象，理解为什么说COMSOL Multiphysics是实现有限元方法的一个很方便的工具。有限元方法的本质是将任何场变量约束以弱形式来表示。理解什么是弱形式（以及为什么数学上要求必须是弱形式），就必须要理解系统约束的强形式是偏微分方程系统和适当的边界条件。它为什么强？由于场变量需要连续并且要有与方程阶数相同的连续偏导数。这就是强约束。而弱形式对函数满足的约束条件限制要相对弱一些不连续但是可积。来看一个偏微分方程和它的等价弱形式。考虑在三维空间域中的具有一个独立变量u的稳态偏微分方程，通用形式为： (11)我们假设测试函数v是域中定义的一个任意函数，并且满足函数vV。方程(11)两侧同乘以v，并在整个域内积分有： (12)这里dx表示体积微元。根据散度定律，我们得到： (13)当偏微分方程被Neumann边界条件约束时，方程(13)中的边界项消失。这也是有限元方法有Neumann自然边界条件的一个原因。回顾第一章，我们发现有限微分法有自然的Dirichlet边界条件。这导致体积积分变为： (14)这必须对所有的vV都满足。现在来看有限元和基底函数。假设u可以分解为一系列基底函数： (15)例如，如果是正弦、余弦等基本渐进谐波函数，那么(15)就是由一系列傅立叶函数构成的。相反，在有限元方法中，要选择那些只在某一基元中作用的基底函数，即除了某一基元，在其它基元中都为零。图8一维空间临近基元的两段线性基底函数图8给出了一维空间中两个拉格朗日线性基底函数的例子。显然，任意函数u(x)都可以由分段线性基底函数和足够小的基元进行任意精确度的近似。基底函数可以使用更高阶数，但是相应的每个基元也需要更多的未知变量ui。所以未知变量个数随着基底函数阶数的增大而增大，基底函数的个数也随着阶数的增大而增大。对于拉格朗日线性基底函数，在任何单元中都是由两个基底函数组成的： (16)这里是基元中的本地坐标。所以对于N个单元，就有2N个基底函数i。二次拉格朗日单元有三个基底函数： (17)所以对于拉格朗日二次单元，N个单元有3N个基底函数。回顾前面讲的，方程(14)必须对所有的vV满足，下面我们考虑所有由基底函数线性组合构成的函数空间，即： (18)但是由于v线性加入(14)，这就说明，如果以任意基底函数i作为v，方程(14)都成立，那么对于所有基底函数的线性组合(18)也成立，即所有vV都成立。所以(18)等价于一个(k+1)N个方程（每个i都有方程(14)）、(k+1)N个未知变量（ui）的系统，这里k是基元的阶数（k1线性，k2二阶）。下面解释为什么COMSOL Multiphysics有限元方法具有这样的作用。COMSOL Multiphysics建立(k+1)N个方程(14)的集合。首先，注意到(u)和F(u)都是u的通用非线性函数。所以一般没有闭合解。在第一章中，我们介绍了COMSOL Multiphysics中预置的零维问题非线性求解器，即f(u)=0，u是唯一未知变量。非线性求解器基于牛顿方法。牛顿方法的N维矢量方程为： (19)这里U是未知数ui的矢量，L(U)是取代(14)中v的基底函数i的方程系统，即： (20)这里的K(U0)是刚性矩阵刚度矩阵，L(U0)是荷载向量。刚性矩阵刚度矩阵是L的负雅克比矩阵： (21)所以(20)是给定近似解U0情况下U的线性方程。因子此，如果U0能够足够接近真实解，线性方程(20)会找出更加精确的解U，该过程可以不断重复，直到最后解足够精确。显然，牛顿法非线性求解器是COMSOL Multiphysics偏微分求解器的核心。COMSOL Multiphyiscs可以实现偏微分方程有限元分析的所有步骤。如果可以的话，它用符号组成非线性运算符L(U)的雅克比矩阵，如果不行的话，采用数值方法建立雅克比矩阵。如果偏微分方程本身是线性的，就不会很复杂。但是建立刚性矩阵刚度矩阵是一个非常庞大的工作这在有限元分析中很常见。不久以前，有限元分析，包括网格划分和建立刚性矩阵刚度矩阵还经常是很多博士的研究课题。对于新的偏微分方程组合，甚至变系数问题（例如1.2节介绍的准线性系统），建立刚性矩阵刚度矩阵的工作量大的吓人。进一步讲，方程(20)的稀疏求解方法和时间积分需要通过编程工作来调整使其适合于某一问题。好在COMSOL Multiphysics将这个工作作为一个子程序（MATLAB函数）来完成，实现与复杂偏微分方程系统的无缝衔接。在介绍边界条件实现方法之前，我们先通过一个常微分方程练习来强调一下刚才讨论的概念。通过该练习了解一下弱形式，找出给定二阶常微分方程的近似数值解。练习2.5有限元详细计算实例。使用有限元方法核心变分原理求解一个简单的常微分方程，详细介绍计算过程中出现的概念。问题比较简单，二阶常微分方程如下： (22)边界条件为：使用了弱方程式。该二阶常微分方程解析解如下（试证明之！）： (23)因为已经知道了解析解，就可以比较近似解的误差。建立弱方程式的第一步是假设一个权重函数和试探函数。将U(x)作为试探函数，(x)作为权重函数。我们一会儿再讨论这些函数的准确形式。试探函数U(x)构成了常微分方程的解。所以如果将U(x)带入(22)式，可以得到残差为： (24)接下来的工作就是最小化残差。要最小化残差首先要计算权重残差，将方程(24)两侧同乘以(x)，然后在整个域中进行积分（即0x/4） (25)使用分步积分法简化上式为： (26)为后续计算方便，我们需要做一些关键假设。因为只要满足边界条件，我们可以任意设定U(x)和(x)，我们令U(x)(x)。这被称作Galerkin方法。如果U(x)(x)，给出Rayleigh-Ritz公式。我们需要选择一个x的代数函数以满足边界条件u(0)=u(/4)=0。 (27)假设N个方程为： (28)所以 (29)这种选择满足边界条件，且不需要考虑方程包含的项数。由于U(x)(x)，权重残差为： (30)将(29)式带入(30)并积分，得到一个与x无关的R表达式。但是该表达式包含N个未知变量ui。我们需要计算ui的值使得权重残差最小。由(29)式得： (31)所以我们有 (32)然后通过将R对uj取导数，最小化残差。对于预定义的N，产生了N个需要同时求解的代数方程。 (33)当N2时只有两个未知变量，u1和u2，产生两个线性方程： (34)矩阵形式表示为： (35)通用形式为： (36)这里K是雅克比矩阵（刚性矩阵刚度矩阵），u是未知向量，L是限制向量（负载向量）。方程(35)的解为：所以方程(22)的解为： (37)如果进一步假设N3，得到含三个未知变量u1，u2和u3的代数方程，矩阵形式为： (38)方程(38)的解为：所以方程(22)的解变为： (39)图9给出了(37)(39)式和解析解的比较。可以清楚看到，只要取的项数足够多，近似代数解能够与精确解很好的吻合。图9 (39)式曲线。三项几何方程能够与解析解很好的符合通过本例，你应该对解的弱形式有了清晰的概念。下一节我们将讨论如何在COMSOL Multiphysics中实现边界条件。3.1 边界条件正如前面例子中提到的，读者应当注意到刚性矩阵刚度矩阵K等价于Neumann边界条件。在第一章中学过，纯Neumann条件将会导致刚性矩阵刚度矩阵奇异，这样COMSOL Multiphysics就不能直接求解，因为它会对用投影法求出的刚性矩阵刚度矩阵特征系统的解加上任意常数。有限元方法的一个特别之处就在于边界条件的处理。我们打算像有限微分法那样处理边界条件。矩阵方程中的某些行将会被对未知变量ui的直接约束所取代，以保证矩阵的阶数。但是这会破坏刚性矩阵刚度矩阵的稀疏性，对一些边界条件有不利影响，所以需要人为的全矩阵求解器。对于相同矩阵方程，全矩阵求解器与稀疏求解器相比计算精度和效率都会有所降低。通过拉格朗日因子，有限元方法可以得到一个很好的解，这是处理等式约束优化问题非常著名的方法。假设除了偏微分方程约束，我们还有一系列边界条件，对所有的vV都满足弱形式。通过将每个基底函数展开并写出边界积分，类比偏微分方程约束(19)，可以得到边界约束的矢量方程： (40)该约束的剩余方程并不需要N个方程。通常它只是其中包含N个未知变量的少数几个方程，因为不是所有的基底向量都会对边界约束起到测试函数v的作用。方程(40)的线性形式类似与(20)： (41)这里N是M的负雅克比矩阵： (42)下面需要一定的处理技巧。刚性方程乘以未知向量，是拉格朗日因子乘以NT，这里T表示转置： (43)这里有什么技巧呢？如果约束(40)成立，将会有唯一的拉格朗日因子满足(43)式。但是(41)和(43)式不只适用与边界条件，任何约束内部逐点，域积分，边或边界等所有能够以弱形式表示的约束都可以用拉格朗日因子法处理。拉格朗日因子那么拉格朗日因子如何满足M(U)=0呢？根据变分原理。当0时，(43)式等价于对下式应用最小化原理： (44)如果想确认约束(40)是否也同时满足，我们对(44)额外增加一个权重项，确保不满足M(U)=0。权重项就称作拉格朗日因子。 (45)下面用线性化M(U)=M(U0)+N(U0)(U-U0)来简化(45)。注意到常数项对最小化不起作用。 (46)最小化(46)等价于(43)，即(43)式的解能够最小化(46)。在有限元方法中，(46)式相当于“能量最小化”，即用刚性矩阵刚度矩阵K作为权重。注意不要与最小二乘法的最小化混淆，它与(45)式类似： (47)线性化L和M，并经过整理并忽略常数项为： (48)所以根据线性代数理论，(48)式的解U就是正则方程6的解： (49)也是下式最小二乘法的解 (50)所以从最小二乘误差意义上讲，方程(50)的解满足约束(40)，M(U)=0；同时从最小能量意义上讲，约束(43)也满足(40)。为使求解简单、计算精确，COMSOL Multiphysics用(43)代替(50)。两者的差别非常明显，最小二乘误差方程(47)适用于通用非线性算子L(U)，而刚性能量最小化方程(45)只适用于对称正定矩阵K。否则，(43)式中的拉格朗日算子就只是出于方便，并不能保证约束(40)满足近似条件。(50)是强更强的条件，但是是以额外矩阵操作为代价的。(47)表明M(U)不是惩罚性的约束，所以更强的条件必须明确将每个约束分别看作是一个惩罚7 (51)该方法没有将计算结果简介表达为矩阵方程的形式，因为约束项含有N，M，K和。弱条件下面来看如何处理2.1节中涡量的点源。它满足弱条件，只计算如下积分， (52)并对刚性矩阵刚度矩阵和(44)中的荷载向量有适当贡献。3.2 基元有限元方法的一个基本思想是将每个域看作由很多个各种形状的更小子域组成。这些子域就被称作有限元。基元的角点被称作节点，场变量在节点上进行计算。在角点之间也可以有节点，称作边节点。在COMSOL Mulitiphyiscs中划分网格时，将整个计算域分解为一个个可选基元并且构建相应的节点。在实际应用中可以找到一百多种类型的基元。作为初学者，可能会对使用哪种类型基元和基元的个数感到困惑。离散化过程设定了基元的类型和个数。基元的个数与解的准确性直接相关。基元个数越多，误差越小。但是大的基元个数也会增大计算成本，需要更大的内存空间和更长的计算时间。经常见到适当往多的设定基元个数，因为没有优化基元个数的公式。一个域中选择的基元个数完全靠经验来决定。虽然准确性随着基元个数N的增大而增大，但是存在某个数Nc，当N大于Nc时，精确度敏感性可以忽略。图10 归一化误差随基元个数变化曲线。随着N的增大准确性提高，但是当超过某特定值Nc时，改进效果可以忽略。图10给出了归一化误差随基元个数N变化的曲线。在每次迭代时基元个数翻倍。可以看到，最后三个点在计算精度上并没有明显的改进。可以通过适当增加一些运算来得到使用的基元个数。在很多情况下，人们对某个区域更感兴趣，而不是整个计算空间。例如考虑液体流过圆筒的过程，希望研究圆筒壁面的边界层现象。在这种情况下就可以、并且应该增大圆筒壁面处的基元密度，降低远离壁面处的基元密度。这样就可以在不增大基元个数的情况下提高准确度。COMSOL Multiphysics允许使用这种网格拉伸的方法。在这点上还需要一提的是，拉伸网格中基元发生偏斜。偏斜基元使得雅克比矩阵不可用。所以在使用拉伸网格时必须非常小心。基元的类型取决于需要求解的问题。而域的维数决定了基元的维数。最简单的一维基元就是在两个端点间划分的一段直线。最基础的二维基元是三角形基元，在三维中就是四面体基元。表6给出了实际应用中对应维数下的一些基本基元。虽然我们在节点间一般使用直线段，但是曲线段能够提供更为通用的基元形式。表6 基本基元维数形状1-D2-D三角形矩形菱形3-D四面体正六面体不规则六面体虽然基元的几何形状是重要的区分方式，但是基元常根据对应的插值多项式进行分类。根据这种分类方法可以将其分为三种类型的基元：单纯型，复合型和多元型。如果多项式在边角节点处只有线性项和常数项，就称为单纯型基元。而符合型基元在角节点、边节点和内部节点处使用更高阶的多项式（平方、立方、五次方等等）。而多元型基元使其边界与坐标轴平行，并使用更高阶的多项式。正如前面提到的，复合型基元使用更高阶的多项式。而多项式的组成和节点的结构则由所使用的Pascal三角形，Pascal四面体和Pascal超立方体来决定。表7给出了05阶多项式的Pascal三角形。而采用的多项式必须要完整，即包含从最低阶到最高阶的每一项。例如c1+c2x+c3x2是完整的，而c1+c2x2就不完整，因为缺少一次项。对于任意二维基元，只要包含了表7中对应阶数的所有项，就可以获得所需的完整多项式。通常对于多项式的每一项都应该有一个节点。例如立方形基元就应该在每个边上设定四个节点，但是也有更复杂的结合情况。表8给出了二维三角形的线性、平方和立方型基元。表7 二维基元的Pascal三角形1常数x y线性x2 xy y2平方x3 x2y xy2 y3立方x4 xy3 x2y2 x3y y4四次方x5 x4y x3y2 x2y3 xy4 y5五次方表8 1-3阶基底函数复杂基元的类型线性平方立方COMSOL Multiphysics提供了两种忽略几何形状的主要基元类型，称作拉格朗日基元和厄密基元。这种命名方式主要来自它们所使用的差值多项式。拉格朗日基元使用拉格朗日多项式作为基底函数。假设在一维基元中，使用拉格朗日多项式Ln表示场变量u(x)，则有： (53)这里un是未知系数。Ln(x)为： (54)展开式给出了对应阶数的多项式。拉格朗日基元通常出现在计算流体力学中。它们可以提供节点处的变量值。厄密基元使用厄密多项式插值场变量。拉格朗日基元和厄密基元的主要不同点在于有效自由度（DOF）。拉格朗日基元的自由度是节点处的函数值（由节点处的所有变量构成）。而厄密基元除了节点处的函数值，也考虑了角点处变量的一阶导数。由于第n个和n1个节点的值u(xn)和u(xn+1)，导数和未知，所以需要用四个未知数来近似u(x)。 (55)由于xi和xi1为已知坐标，所以可以轻易写出四个方程：两个节点的函数值和两个节点的一阶导数。 (56)这里du表示节点处u对x的导数。通过矩阵求逆，n可以用节点值和导数值表示，方程变为： (57)这里i(x)是厄密插值函数（该例中为立方函数，也称为三次样条）。如果你仔细观察会发现，立方插值函数在拉格朗日基元中需要四个节点而在厄密基元中只需要两个节点。厄密基元通常用在求解桁架的承载和受力情况。对于这两种类型，COMSOL Multiphysics进一步提供了曲线网格基元和Argyris基元。曲线网格基元主要用于精确近似真实边界条件。Argyris基元是五阶的厄密基元，使用节点值和高达二阶的导数。它也在边的中点使用常规算子u。Argyris基元的五阶多项式需要确定21个常数。除了预定义的基元，COMSOL Multiphyiscs也允许用户自己定义基元。感兴趣的读者可以参考COMSOL Multiphysics手册，手册中详细介绍了如何定义一个新的基元类型。这里我们对初学者如何理解和使用基本基元做了详细的介绍。使用COMSOL Multiphysics，我们不需要非常精通网格划分技术和基元设定。正如前面看到的，一旦确定了偏微分方程的求解域，只需要点击一个键就可以完成网格划分工作。但是如果读者对基本概念感兴趣，可以参考8及其引用文献，里面详细介绍了基元和网格划分技术的发展。练习2.6均匀介质中的热传递。为了介绍有限元方法如何用公式表达和具体的求解方法，下面来考虑一个例子。第一章6.1节介绍了一个非均匀介质中的传热问题。这里我们考虑一个更简单的问题。在整个域内认为传热系数为常数，f(x)不变，域的长度L1。图11显示了该物理系统。我们假设穿过任意垂直与中心轴线的横截面的热流均匀（如图11所示），那么温