COMSOL_弱形式详解

COMSOL Multiphysics 弱形式入门物理问题的描述方式有三种：1、偏微分方程2、能量最小化形式3、弱形式本文希望通过比较浅显的方式来讲解弱形式，使用户更有信心通过COMSOL Multiphysics的弱形式用户界面来求解更多更复杂的问题。COMSOL Multiphysics是唯一的 直接使用弱形式来求解问题的软件，通过理解弱形式也能更进一步的理解有限元方法（FEM）以及了解COMSOL Multiphysics的实现方法。本文假定读者没有太多的时间去研 究数学细节，但是却想将弱形式快速的应用到实际工程中去。另外，本文也会帮助理解 COMSOL Multiphysics文档中常用的到一些术语和标注方法，相关理论可以参考 Zienkiewicz1，Hughes2，以及 Johnson 3等。为什么必须要理解PDE方程的弱形式？ 一般情况下，PDE方程都已经内置在COMSOL Multiphysics的各个模块当中，这种情况下，没有必要去了解PDE方程和及其相关的弱形式。 有时候可能问题是没有办法用COMSOL Multiphysics内置模块来求解的，这个时候可以使 用经典PDE模版。但是，有时候可能经典PDE模版也不包括要求解的问题，这个时候就只 能使用弱形式了（虽然这种情况是极少数的）。掌握弱形式可以使你的水平超过一般的 COMSOL Multiphysics用户，让你更容易去理解模型库中利用弱形式做的算例。另一个原因 就是弱形式有时候描述问题比PDE方程紧凑的多。还有，如果你是一个教授去教有限元分 析方法，可以帮助学生们直接利用弱形式来更深入的了解有限元。最后，你对有限元方法了 解的越多，对于COMSOL Multiphysics中的一些求解器的高级设置就懂得更多。一个重要的事实是:在所有的应用模式和PDE模式求解的时候，COMSOL Multiphysics 都是先将方程式系统转为了弱形式，然后进行求解。PDE问题常常具有最小能量问题的等效形式，这让人有一种直觉，那就是PDE方程都 可以有相应的弱形式。实际上这些PDE方程和能量最小值问题只是同一个物理方程的两种 不同表达形式罢了，同样，弱形式（几乎）是同一个物理方程的第三个等效形式。这三种形式的区别虽然不大，但绝对是很关键的。我们必须记住，这三种形式只是求 解同一个问题的三种不同形式一一用数学方法求解真实世界的物理现象。根据不同的需求， 这三种方式又有各自不同的优点。PDE形式在各种书籍中比较常见，而且一般都提供了 PDE方程的解法。能量法一般见 于结构分析的文献中，采用弹性势能最小化形式求解问题是相当自然的一件事。当我们的研 究范围超出了标准有限元应用领域，比如传热和结构，这个时候弱形式是不可避免的。化工 中的传质问题和流体中的N-S方程都是没有办法用最小能量原理表述出来的。本文后面还 有很多这样的例子。PDE方程是带有偏微分算子的方程，而能量方程是以积分形式表达的。积分形式的好 处就是特别适合于有限元方法，而且不用担心积分变量的不连续，这在偏微分方程中比较普 遍。弱形式也是积分形式，拥有和积分形式同样的优点，但是他对积分变量的连续性要求更 低，可以看作是能量最小化形式的更一般形式。最重要的是，弱形式非常适合求解非线性的 多物理场问题，这就是COMSOL Multiphysics的重点了。小结：为了理解PDE方程的弱形式，我们必须跳开常规的偏微分形式，对于积分形式要好好研究。由于最小于能原理对比弱形式来说好理解的多，所以我们将从线弹性开始学 习，依次到热传导，电流传导等问题。这几种物理问题都有相关的能量和功率可以进行最 小化。我们将只涉及到静态问题，重点是在结构分析和更特殊的线弹性分析。弹性静力学PDE及其弹性能量方程在静力结构分析问题中，我们需要求解的是Navier方程-v-c = F其中。是应力张量F是体力，比如重力等。如果不习惯用张量的形式，你也可以将张量展 开写成矩阵形式。这个方程表示了力（或者等效力）的平衡，实际上是三个方程的合并形式 3D中每个坐标方向有一个方程。计算区域记为Q，其边界记为0Q。应力张量c和应变张量之间的关系称为本构关系，线弹性本构一般遵循胡克HOOK 定律。=CsG其中c是弹性张量，这个关系式说明材料的行为实际上和弹簧差不多（前提是线弹性）。最后，我们可以将应变矢量和位移的关系表述出来 =Vu这里u指的是位移矢量u=（u,v,w），其定义就是变形体上的材料点和未变形时候的位移差。总结以上所有的方程，我们得到了一个二阶PDE方程（Navier方程），-V- （cVu） = F需要一个边界条件来求解，n - （cVu） = P其中n是表面ae的法矢，p是边界上的面力或牵引力。后面会介绍更多边界条件。这个PDE方程的弱形式为，j VuVVVfl-Q-j* E&Q- j PTxdS = 0其中v=（七,七,七）称为试函数。注意，尽管Navier方程是一个矢量表达式，但是上面的表 达式是一个标量形式。下面介绍如何去推导以及理解弱形式。弹性势能在结构分析中，PDE方程及其弱形式的表达式都不太常见，相反，能量最小化形式因 为其直观的表达形式用的较多。这类问题的能量积分形式对应于总势能的最小化，即对象中 存储的弹性能。总弹性能是一个标量，可以写成：dV单位是体积，因此积分出来的单位应弹性能表达式同样适用于非线性问题。在这些表达式中，我们假设体力F为零，并忽 略了边界效应。这些影响可以在以后引入。积分的意义是每个体积微元的内能总和，其中应 力张量单位是Pa，微元体上的应变d没有单位 该是Nm。如果问题是线弹性的，则可以显式的写为：WE = j j(re)r/er/Q（C萨）8如利用下面的通用公式:用应变张量替换上式中的标量变量X，弹性张量匕替换上标量常量a。联立上面的式子得到:我们用c代替c来配合COMSOL Multiphysics手册中的标记方式。再提醒一次，如果你不习惯用张量，可以将张量看成是一个3X3的矩阵，点乘是一种张量的运算符号，弹性张量c 是一个4阶张量（看上去就像4维矩阵）。更多的标记方法可以参考COMSOL Multiphysics的 Anisotropic Structural Analysis 中的 Matrix Notationo弹性能积分形式下的单位说明： =无单位N 八 ,c = Pascalm2j dQ = m 3Q最终给出总的积分单位是N-m能量。气的表达式就是我们通常说的能量泛函，即位移矢量u（或实际上是u的梯度）的泛 函。这种函数的函数，而不是坐标的函数，通常被称为乏函，比单元微积分和多元微积分更 加抽象。与积分类似，我们可以说气就是函数u的泛函：u = r VuJ cVurfQQ这好比是一个2D的变量x, y的二元函数：/（X） =Lg其中 x= 3 y）, A = b, b = d。 c d e采用这样的类比是因为在后面我们会看到矩阵A与有限元的刚度矩阵比较类似。我们要说明一下函数和泛函的一些区别，古典分析中的函数概念是指两个数集之间所 建立的一种对应关系，现代数学的发展却是要求建立两个任意集合之间的某种对应关系。 函数概念被赋予了更为一般的意义，通俗解释泛函指的就是“函数的函数”。在这里定义域 为Q，泛函可以在整个定义域内进行微分积分等操作。泛函的变量是函数，这个函数也是有容许空间的。如果函数u可以变化，可能会产生一 些不符合物理规则的一些现象，例如结构的刚性位移等。比如一个对u的基本约束就是材料 不能穿越本身。在有限元分析中，泛函一般是某种能量积分，比如弹性能。对于其他的物理场，可能是 其他的能量积分，或者是一种等效于能量的标量也可以。至于积分区域，一般由分析对象的 CAD几何区域所确定。静态电流传导和能量的生成在静态导电问题中，PDE方程由最基本的保守形式开始：V- J = 0其中J是电流密度。材料(或本构)模型采用欧姆Ohm定律：B =。E其中E是电场，。尸是电导率。另外，已知：E = -W其中V是静电势，综合以上式子得到-V- (q VV) = 0在 COMSOL Multiphysics 中，这就是所谓的 Conductive Media DC 方程。电阻产生的热能稳态电流的能量问题是在电导体中的电阻热屿=jjdEd。=诉虹 dQJJC1EQE其中J表示电流强度，E代表电场强度，。是一个二阶电导张量(3X3)。如果导体是金属， g电导张量一般是一个对角矩阵，如果是晶体，情况就复杂多了。尽量减少电阻产生的热量，也就是减少热损耗，是我们要研究的一个最小值问题。 如果问题是线性，则积分可以显式地写成：=打QEQ因为E = -W，其中V是电势，可以得到：q将这个式子与结构力学中的式子进行对比，发现他们非常相似。V的梯度对应于位移梯度，电导率张量b对应于弹性张量C。在稳态电流和结构力学的计算过程中，张量形式都可以 g改写为矩阵形式。传热PDE方程和能量形式对于稳态传热问题，PDE形式为：-V -( kT = Q其中T是温度，k是热传导系数，Q是空间分布的热源。热能基于传热方程的典型泛函为：其中T是温度，k是热传导系数张量(3X3)。泛函极小值泛函极值的概念借用了微积分中的不少方法。本节首先会介绍函数微积分的求极值方 法，接下来，我们会借用有限元中常用的术语和标注方法来推导我们熟悉的结果。这个过程 可以被看作是微积分方法的一种推广。考虑一个多元微积分函数f,我们要求最小值:寻找x使得f（x）最小化这里x是一个矢量，或者点的坐标。通过微积分我们知道，这个时候首先必须求函麴的梯 度。将梯度的设置为0，我们可得到一个非线性方程组。求解方程，我们可以得到一系列的 坐标点x，如果在其中某点处的二阶倒数（一般称为Hessian矩阵）为正（或者说有正的特征值）， 就说这点就是我们要求的极小点，就好像该点是整个函数的一个谷底一样。利用Taylor展开的观点，假设已知一个最小值x，我们可以在上面施加一个小的扰动， 由Taylor展开可得：fix-dx） =y（x）-V/l）6x-3xrH（x）3x-,这里H就是前面所说的Hessian矩阵。现在我们用其他的方法来说明函数f在x最小。首先， 假设x是一个极值点，当添加了一个8X后，f对于其一阶值不改变。换句话说，如果我们在 x上添加一个&来扰动f，其一阶Taylor级数应该为0。这个条件应该对每个方向都是成立 的，否则该点就不是极值点了。如果上式第二项为0：y/（ X = 0对于任意小的&都成立，也就是：以 X）= 0我们这里只是用一个稍微有点不同的方法得到了一个同样的结果。但是，这只是给了我们一个极值点的信息，如果要确定其是最小极值点，必须保证第三 项（二阶项）对于任意&都为正：3xJ，H（x）3（x）0只有当H的特征值都为正时，上式成立（参考线性代数）。有可能会遇到二阶项也总为0， 这个时候我们必须借助更高阶项来判断极值点。下面是函数f的一个特例：二次多项式：必）=* Ax-b-g其中A是对称矩阵。如果我们应用Taylor展开，可得到：或者/（x + 5x） =+ 宰/&一这里零阶，一阶和二级项都在独立的中括号内。为了得到一阶变分，矩阵A必须是对称的。 极值的条件成了：对于任意小5x都必须成立，则上式成为:这里我们对矩阵进行了转置，而且利用了矩阵A的对称性，即AT = A。极小值的条件也就是矩阵A必须是一个正定矩阵，如果矩阵A是负定矩阵(只有负特 征值)，则得到极大值。如果A是不确定的(特征值有正有负)，则极值可能是一个鞍点， 既不是极大值，也不是极小值。如果矩阵A是对称的，而且正定，则函数f是超椭圆的。 在2D中，超椭圆就是椭圆。二次多项式的几何特征影响经典的PDE方程和有限问题的分 类。当利用有限元方法去离散一个椭圆的PDE问题时候，得到一个对称矩阵(刚度矩阵) 的线性代数系统。这样的问题一般等效于最小能量问题。弹性静力学问题变分我们将通过两个步骤来介绍最小能量法理论。首先粗略说明，让大家熟悉基本概念；接 下来考虑细节。还是以线性静态问题为例，因为这是所有有限元理论都会提到的，从而更容 易进行比较。理论概述让我们回到线弹性问题的弹性能泛函表达式：吨(11) = | |(Vu)rcVurfQ这里的位移矢量u和前面讲的微积分中的点矢量x的角色类似。要寻找能量泛函w的最小值，我们首先必须得在u上施加一个扰动5u :W u - 6 u)= | u - 6 H)/ rV( u + 6 n J dQ.=ij(Vu)JVut/Q- j(V6u)JcVuQ- j(Vu)JcV3uc7Q +。口(V5u)JcV3ut/Q上式中两个中间项实质上是一样的(因为c的对称性)，所以我们可以写成：W(u - 6u ) = | (+ J( Vu)J +1 rr口-j(V6n) cVSuO将上式和多元函数表达式对比，我们发现寻找极值点就是找一个使二次项为零的u：j (Vu)cV6iifi/Q = 0其中5u是任意的。如果我们要寻找的是极小点，则还必须有：(V6n)ZrV3uQ0Q第二项就是泛函的一阶微分：3 Wv u ) = ( Vu J rVSut/Q第三项成为泛函的二级微分:和前面一样，为了寻找极小点，我们必须保证对于任意5 U第一阶微分为零，二阶微分 为正。这种寻找最小势能函数的方法也可以称作虚功原理。另外还有一种方法就是初始的时候将扰动写成a5 U，这时对于任意可取的5 U，其能量函数写成W(U + a5 U)。回到微积分的基本概念，去寻找W对于a的极值点：- c(6u) =|( Vu cVudQ. cz J (Vu)J cVSut/O -QQ:c?j (宁如果我们将它看成是对于a的Taylor展开，就可以找出其一阶导数(对于极值点必须 为零)，由于5u是任意可取的，我们可以得到和前面相同的结果。小结上面的过程省略很多推导步骤，如果大家对推导有兴趣，可以试着自己推导。我们要说明一下的是：1、变量8(而不是它的梯度)必须是很小而且是任意的。2、这里没有考虑边界条件和体力，比如重力等等。我们前面所讨论的问题局限于一个没有 任何约束和载荷的边界条件的区域上。3、一般来说诵的限制比多元微积分中8x宽松。在泛函中，只要8u是在容许的范围内即可， 也就是u +8 u必须和物理位移场相对应。理解这个意思对理解有限元弱形式非常重要。考虑边界条件和体力如前面所讲，弹性能的泛函形式是不完整的，因为它没有加上相应的边界条件和载荷。 弹性能的单位是N - m，也就是力乘上位移。在边界上，我们一般施加面力，或者指定位移， 单位为m。一般来说，我们希望附加形式是“面力乘上长度，。同样的方式可以对体力进行 处理F。在数学上，结构场的边界条件分为两类。第一类直接定义边界上的力：n - i/Vu) += g其中第一项由定义域内的方程所确定，第二项称为弹簧常数q，等式右边是面力g。这种边 界条件就是我们通常说得流量或者Nuemann边界条件。第二那边界条件就是R = hvi-j- = 0定义一个固定的或者Dirichlet边界条件。如果h是矩阵的形式，r就是定义了边界上的指定 位移。固定边界条件不能直接加入泛函中去，但是可以通过反力间接加上去。当指定位移边界时，可以描述一个反力P( N/mm2)，也就是弹性体可以在固定处保持不变。反力就是我们这里用到的Lagrange乘子，通过添加反力到力作用处的边界，可以忽略到固定边界类 型。这时候我们可以形成统一的边界条件：n - ( l Vu .1 + * u = g 4 RCU这里R是原始的固定边界，P是需要计算的反力。在前面的简化形式中，h和尸都是常数， 所以上式可以变化为：ti i rVu) + qu = g -记住，方程中的每一项都是矢量，表示各个方向的面力。为了得到所做的功(能量)， 必须点乘上位移u。通过合并一些系数项，将外力写成P，可简化表达式，这时边界条件可以写成：n - (cVu) = P对于其他物理场，可能P代表边界上的源项。注意到上式和Navier方程非常接近：-V (rVu) = F将能量泛函展开:关键推导这个时候，我们又要在u上添加上8u，可得：Rz(u-6u) = J u + V(u + 6ujr/Q-+- j Ku+ 跖)肉8*u=JF零阶项就是泛函本身，第一阶项是：汕0- j以口nco这个方程是非常重要的一项。接下来我们详细介绍。从前面的讨论可知，我们应该重新组合多项式，保证带有8的被积函数成为一项。如 果可以做到，因为8u是任意的(事实上必须是在容许范围内)，我们知道这一项必须为零。 这是我们能找到极值点的唯一方法。右边第一项需要进一步处理得到我们需要的形式。第一项我们可以根据Green公式(有时候可能采用的是Stokes原理)进行分部积分：vuVV(6ur?Q = J (n (rVui)ur?1.?- J(V (cVu)7GuQ利用c的对称性，我们可以得到：J VTcV(S)rfQ =利用Green公式得到：15无3)= J (d lcVu)：/6uS-仔丁如仙j 而SF 口QQcQ将体积项和边界项合并起来：%何)=J(*EUF代何- j (ng-Fg 口沁确定极值点，必须有：1次=j(-V icVu)-F)T5iiQ- f iu-Vu-P：j6uf?5 = 0上式应该对于任何8u都成立。因此体积项必须有:V icVu) = F边界项上有：u (Eu) = P现在我们又回到PDE问题上了，通过能量最小化原理又重新推回到了 PDE形式上！这也是 说明最小能量化和PDE形式本质上是统一的一个数学证明。弱形式那么，到底什么是弱形式呢？ Navier方程的弱形式实际上已经在前面的推导过程中出现 过了，即一阶变分的原形式：I6TF-111) = RuLw.Bu /Q-冬丁如忒1- J PTducS 口nco如果我们回到COMSOL Multiphysics的文档(或者是关于有限元和弱形式的书籍中)， 会发现所谓的试函数相当于扰动8U，v = 6 u弹性静力学PDE方程的弱形式为了更好的理解弱形式，我们必须丢弃前面讨论的能量最小化原理，转向一种更加抽象 的方法。弱形式之所以比能量最小化原理更强大，是因为它还可以应用到一些没有得到较好 的能量定义的问题中。首先我们考虑弹性静力学的PDE方程V ieVu) = F边界条件是：u (i?Vu = P抽象的过程如下：乘上容许范围内的试函数V，在感兴趣的域内积分可得：p：V (eVu )TF = JFTvdV对左侧利用Green公式进行分部积分：| (q= 住匠应用PDE方程的边界条件，可以得到：5口口口整理可得:JVuVVvr,n-j PTvdS = 0 n典w这就是PDE方程的弱形式。如果在积分区域内对于试函数v都是有效的，则上式和PDE 方程是等效的。PDE方程的解称为强解，而弱形式的解称为弱解。二者唯一的区别是弱形 式对于积分参数的连续性要求比PDE形式低。由于变形梯度和弹性张量在弱形式里面都不 需要微分，所以对函数连续性要求没有那么严格，而在PDE形式中，所有的变量都处在散 度的算子下，这要求这些变量必须是可微的。在弱形式中对于可微的要求放松了（一阶）。同时，注意到弱形式和前面的一阶变分形式保持了一致，弱形式也可以作为虚功原理的 一种推广。只是虚功原理中的位移5u换成了更加抽象的试函数V。如果弱形式解和能量最 小化原理不一致的时候，极值点变成了鞍点。也就是，在弱形式中，仍可以将试函数理解为 一种推广了的虚位移。一般性问题的弱形式正如前面所提到的，弱形式只是PDE方程的一种推广形式，它对变量的连续性要求比 较低。那么能量方法呢？如果有一个定义好了的能量来最小化，那么能量法和弱形式是一致 的。但是，在下列情形下，弱形式更具有适用性：假如PDE方程没有相对应的能量可以进 行最小化。在这种情况下，弱形式仍然是适用的。由于弱形式对解的要求较低，所以说弱形 式比PDE和能量最小化适用范围更广泛。我们将给出一个没有对应能量最小化的PDE的例子。对流一扩散PDE问题对流一扩散PDE问题没有与之相对应的可最小化的能量：这里c是扩散系数，P是对流系数，a是反应/吸收系数，f是源项。变量u是标量函数，代表浓度（在COMSOL Multiphysics手册中的Convection-Diffusion模块中，浓度是用变量 c表示，扩散系数用D表示）。在这里我们考虑Neumann边界：所有困难将集中在刚度K的提取上，主要是对u和Lagrange乘子的线性表达式的集成。 为了得到弱形式，将PDE方程乘以一个试函数v，积分：j一 * Wm） + & 研一网 E萍=月沽F口Q这里的试函数v是一个标量函数。将第一项分部积分，并将所有的项都移到左边，可得到：J (r 研。. n (+pp-研小wF f audV- fMF = 0Q口口口口加上边界条件，得到：jcSubdomain Setting，在 Weak Term编辑框中输入：-(cKux ux_test+tuy uy_test+cuzuz_t&st+ f bx ux+byruy+bzuz+au -f) =u_t日st)边界设定，Physics-Boundary Setting, Weak Term编辑框中输入：Pu testCOMSOL Multiphysics将边界设置和子域设置分开，因为子域和边界上可以设置不同的 数值积分算法。弱项，也可以在如果想要扩展内建的经典PDE模板或者物理应用模式（比如传热）Physics-Equation System中对应的对话框中输入相同的表达式。弱形式方程会自动添加在控制方程中。（通过设置所有的PDE或材料参数为0，选择齐 次Neumann边界（流量=0），可以去掉应用模式自动创建的弱形式。）Dirichlet或者固定边界，在Boundary setting对话框中的constr编辑框输入弱形式， COMSOL Multiphysics会添加相应的Lagrange乘子（参见用户手册中的边界条件章节）。结构力学PDE问题静态结构力学的基本方程是Navier方程：V-kVu） = F边界条件:对流一扩散方程中的标量项现在全部成了矢量和张量，Navier方程的弱形式为：10 = - |VurcVvr,Q + jFr7Q + j FE5约定标记如下： 矢量u的分量：u，v和w。 位移矢量梯度Vu的分量： ux，uy，uz， vx， vy， vz， wx， wy， wz。 试位移矢量v的分量：u_test， v_test， w_test。 试位移矢量梯度Vv的分量: ux_test，uy_test，uz_test，vx_test，vy_test，vz_test，wx_test， wy_test，wz_test。弹性张量的分量：c11，c12，c13，c14，c15，c16，c22，c23，c24，c25，c26，c33, c34， c35， c36， c44， c45， c46， c55， c56， c66 体力矢量F的分量：Fx，Fy，Fz。 边界面力矢量P的分量：Px，Py，Pz。在子域内，弱形式输入为：-ex_testsx-ey_testsy-ez_testsz-2exy_testsxy-2eyz_testsyz-2exz_test “sxz+Fxu_test+Fyv_test+Fzw_test其中ey=vyez=wzexy=0.5T(uy+vxeyz=O.5(QXZ=O.6(U Z+WX)1&x+c1ez+2o14xy+2d 5&yz+2d 6ex2&y=c12 ex+c22ey+c3 Kez+2c24?,exy+2 KG26eyz+2az=c13ex+c23Koy+c33 Tez+2c34exy+2 Tc36eyz+2T c36kdxzsxy=c14Tex+c：24xQy+c：34Texy+2T c45eyz+2 tc46xqxz&yz=Gl4Tex+c24ey+c34T e2+2G44c45eyz+2 TG46exzsxz=c16QX+c26ey+c3 6 ez+2C46 exy+2： c66eyz+2 kc66qxz这些表达式定义了应变分量(ex，ey，.)和应力分量(sx，sy，.)。后面带有test后缀的，COMSOL Multiphysics都会和上式一样建立相应的试函数和试 函数梯度的表达式。比如，exy_test等效于05*(uy_test+vx_test)。另一种方式是test()，其 中 test(xy)表示 0.5*(test(uy)+test(vx)，也就相当于 05*(uy_test+vx_test)。对于其他一些张量表述如有疑问，可以参考COMSOL Multiphysics中的Anisotropic Structural Analysis 的 Matrix Notation。如果想更直观的表述弱形式，我们可以用原始定义代替变量，最后变成： -ux_test*(c11*ux+c12*vy+c13*wz+c14*(uy+vx)+c15*(vz+wy)+c16*(uz+wx)-vy_test*(c12*ux +.对于各向同性体，其实cij就是杨式模量E和泊松比v的简单函数。详情参考COMSOL Multiphysics 文档。在边界上，对于载荷类边界条件，弱形式可以在weak编辑框中写成标量的形式：Px*u_test+Py*v_test+Pz*w_test如果采用固定边界，我们必须在其中一个constr编辑框中输入相应的表达式。对于多物理场仿真，约束和载荷在weak和constr中的形式非常重要，尤其是采用弱约 束的时候。更多详情可以参考COMSOL Multiphysics文档以及和Lagrange乘子相关的技术 文档。尽管弱形式是一个标量表达式，但是COMSOL Multiphysics中，弱形式有和PDE系统 一样多的未知量需要文本输入。原因在于不同的多物理场问题可能需要不同的有限元分析类 型和保证其数值稳定型的积分算法。对于3D结构分析，弱形式中有三个文本输入框。但是， 在离散之前，采用了同样的有限元单元和积分类型进行合并，这样就可以选择不同的弱形式 进行操作。例如你可以在第一个域内选择弱形式，而其他的域内设置为空白。对于流动问题的Navier-stokes方程，情况又稍微有些不同。和未知的速度场相比，未知 的压力采用一个低阶有限元来离散。这种情况下，不能将所有的弱项全部在同一个弱域内输 入。为了保证数值稳定型，必须依靠混阶有限元(mixed finite element)。混阶有限元并不是 COMSOL Multiphysics特别制定，而是数值算法所需要的。有限元方法本章说明弱形式如何利用有限元方法来进行离散。假设我们需要离散以下扩散问题：10 = - f(icVu) Vvf PvdSasn这是一个对流一扩散问题的特殊情况，其中a =0，P =0。有限元的基本实现是将整个计算域Q离散为多个特别简单的形状的小单元，比如2D中 的三角形，3D中的四面体等等。相应的网格，例如三角形，由边和节点组成。下一步就是 要选择一个比较容易实现的一些近似方法，其中一种比较简单的方法就是将解表示为采用线 性多项式插值的所谓基函数的和。基函数的构造方法是指定某个节点为1，而相邻的节点为 0，二者之间的值就是从0到1线性变化。这里说的相邻指的是中间有一条边将其连接起来。遍历三角形网格的所有节点(从1到N)。定义节点i的基函数为甲注也就是在节点i 处其值为1，其他点处值为0。注意g只是在节点i及其相邻的三角形内不为零。现在假设 真实值u可以用基函数的求和匕来近似描述：yi = 1参数气是u在节点i的值。同样，我们可以对试函数进行类似处理:v = E 啊；= 1下标h表示离散函数属于由所有三角形边中最长边表示的具有确定的网格尺寸h的网格。由于我们可以任意选择试函数，因此可以将除了 j点以外的所有的七设置为零，接下来我们将所有的试函数(j = 1,.,N)输入到弱形式中去，每个试函数都可以得到一 个方程。这样可以生成一个线性代数系统，系统矩阵就是我们所说的刚度矩阵。为什么我们 可以自由选择试函数，不妨回想一下前面提到的弱形式需要对所有可取的试函数成立。选择 试函数是有限元方法的重要环节，因为他在很大程度上影响着刚度矩阵。由于刚度矩阵中很 多项为零，所以一般是稀疏矩阵。当我们使用试函数的时候，生成的有限元刚度矩阵应该是一个方阵。如果弱形式本身定 义良好的话，刚度矩阵应该是非奇异的，也就是说系统有一个唯一解。现在考虑扩散方程的弱形式：10 = -I ES口M将表达式写成离散形式:方程重新排列：0 =-云|J相计+ 场妣=1 nq sei采用矩阵标注可得：0 =-KU-L或者：KU = L在这里刚度矩阵K是：白=两/勺：出解矢量U的单元为U载荷矢量L的单元为，I上=J *膏 口口现在我们明白为什么选择基函数和试函数很关键了。如果我们关注刚度矩阵K，会发 现其中很多元素为零，因为前面已经提到每个9j都是大部分为零，同样9 j的梯度也是大部 分为零的。有很多有效的算法去求解这类稀疏矩阵，COMSOL MultiphysicS提供一套稀疏线 性系统求解器。有限元方法同样适用于非线性问题。非线性方法一般来说采用迭代的算法，每一次迭代 就是求解一个与上面类似的线性弱形式方法。抽象和几何解释为什么有限元可以解决问题，它是如何解决问题的？前面的讨论中可以找到一些答 案。为了有一个更清晰的答案，我们需要了解一些更多的泛函的概念。我们将发现有限元方 法通过一种优化方法将解投影到一个有限维函数空间来求解。标量积为了得到有限元方法的几何解释，或脑海中的意象，我们需要熟悉标量积的概念。在 线性代数中，我们知道两个矢量f和g的标量积为f g = fXSX fyg. f-S- = Ng)这里的矢量f和g属于3维矢量空间。标量积可以推广到任意维N，(fg) =-f2S2+A4 + 为 = 22= 1标量积有时也称做内积。如果两个矢量正交，(f-g)二 o一个矢量f可以通过标量积投影到另一个矢量g矢量差与g正交:其中fp是与g平等的投影矢量:投影矢量的唯一性特征是原始和投影矢量之间的差与它所投影的矢量正交。如果需要 找到在方向g上与f最近的矢量，fp就是我们的答案。矢量e可以看作是关于f到之间的 近似误差。换句话说，误差矢量e与矢量g正交。后面在对有限元进行几何解释时将用到这 个结论。但首先我们得介绍一些泛函分析的概念。与矢量不同，泛函分析讲的是函数，它们必须属于无限维的矢量空间。我们可以积分 形式定义一个无限维矢量空间（函数空间）中的标量积：（心）1 = j itvdV如果两个函数正交，则有进一步将标量积的概念推广，并且考虑包含函数梯度的被积函数，下标1和2用来说明上面是两种不同的标量积。Hilbert 空间如果对于一个标量积，我们只考虑那些与自身进行的标量积（积分）有一个有限值的 函数u，（iLli） z我们说这些函数属于一个确定的函数类，或函数空间。由标量积和定义域Q组成的函数空间 就被称为Hilbert空间。对应于上面的标量积1的Hilbert空间通常标记为L2，与标量积2对应的Hilbert空间 常称为H1。通常空间L2中的函数比H1中的多，因为H1中的函数自动地存在于L2中，反之 则不一定。我们也可以将这种现象称为H1是L2的一个子空间，或有限元方法的抽象形式现在让我们考虑PDE问题3：-V （研+ w =/在域Q中，边界条件为u = 0。这是对流一扩散方程的一种特例，其中c = a = 1，P = 0。通过分部积分，表明当试函数选择成在边界上具有相同的边界条件v = 0时，弱形式中的边 界项为0。得到弱形式为：0 =-n其中包含前面提到的两种标量积，可改写为：找到u，使得对于所有属于相配的Hilbert空间中的V，有:0 = （. + （源1这里的相配的空间是Hi,且在边界上v = 0。解函数u也必须属于这个空间。注意，对于前 面提到的弱形式，我们可以很自由地添加不同的标量积，因为每个积的结果是实数或虚数。现在我们选择前面讨论的基函数的和（线性组合）来近似u和v，近似解被称作/叫。这些函数属于Hilbert空间，可称为匕,由线性基函数气扩展而来。匕的空间维数为N （基 函数的数目）。此外，匕是Hi的子空间，换句话说，如果函数属于匕，则它自动地属于Hi。空间Hi是一个大得多的空间，因 为它是无限维的。有限的基函数不可能扩散成属于Hi的所有函数。弱形式的有限元现在变成了：对于匕中的所有匕，在匕中找到正，使得。=-（人）2 + 印现在我们对有限元方法在函数空间的作为一个确定的投影的几何解释有了更深入的 了解。最后，在原始的弱形式中：0 = -（）/（31用V代替v。这是合理的，因为V在Hi中，而V在V中，由于V是Hi的子空间，因此Vhh hhh也是在Hi中，o =-（以 + 加）现在从弱形式中减去有限元解，得到：0 = _（ -瞄仇或 * W2 = 0即离散误差：e = ii - tth与所有的Vh中的Vh正交。也就是说，有限元解uh是真实解u在属于Hi的有限维子空间V的投影。最终，我们得到了关于有限元方法的几何解释，对于给定的网格，有限元解u是在函数空间匕中关于标量积（，）2最接近真实解u的解。参考资料:1 Zienkiewicz, Taylor, The Finite Element Method: Volume 1-3, Butterworth-Heinemann; 5th edition, 20052 T. J. R. Hughes, The Finite Element Method : Linear Static and Dynamic Finite Element Analysis, Dover Publications, 20003 C. Johnson, Numerical solution of partial differential equations by the finite element method, Studentliteratur, 1987, ISBN 91-44-25241-2