COMSOL弱形式入门

COMSOL Multiphysics弱形式入门
物理问题的描述方式有三种：
1、偏微分方程
2、能量最小化形式
3、弱形式
本文希望通过比较浅显的方式来讲解弱形式，使用户更有信心通过COMSOL Multiphysics的弱形式用户界面来求解更多更复杂的问题。

COMSOL Multiphysics是唯一的直接使用弱形式来求解问题的软件，通过理解弱形式也能更进一步的理解有限元方法（FEM）以及了解COMSOL Multiphysics的实现方法。

本文假定读者没有太多的时间去研究数学细节，但是却想将弱形式快速的应用到实际工程中去。

另外，本文也会帮助理解COMSOL Multiphysics文档中常用的到一些术语和标注方法，相关理论可以参考Zienkiewicz[1]，Hughes[2]，以及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 方程
σ-∇⋅=F
其中σ是应力张量，F 是体力，比如重力等。

如果不习惯用张量的形式，你也可以将张量展开写成矩阵形式。

这个方程表示了力（或者等效力）的平衡，实际上是三个方程的合并形式——3D 中每个坐标方向有一个方程。

计算区域记为Ω，其边界记为Ω∂。

应力张量σ和应变张量ε之间的关系称为本构关系，线弹性本构一般遵循胡克HOOK 定律
εσεc =
其中εc 是弹性张量，这个关系式说明材料的行为实际上和弹簧差不多（前提是线弹性）。

最后，我们可以将应变矢量和位移的关系表述出来
ε=∇u
这里u 指的是位移矢量u =(u,v,w)，其定义就是变形体上的材料点和未变形时候的位移差。

总结以上所有的方程，我们得到了一个二阶PDE 方程（Navier 方程）， ()c -∇⋅∇=u F
需要一个边界条件来求解，
()n c ⋅∇=u P
其中n 是表面Ω∂的法矢，P 是边界上的面力或牵引力。

后面会介绍更多边界条件。

这个PDE 方程的弱形式为，
其中v=),,(z y x v v v 称为试函数。

注意，尽管Navier 方程是一个矢量表达式，但是上面的表达式是一个标量形式。

下面介绍如何去推导以及理解弱形式。

弹性势能
在结构分析中，PDE 方程及其弱形式的表达式都不太常见，相反，能量最小化形式因为其直观的表达形式用的较多。

这类问题的能量积分形式对应于总势能的最小化，即对象中存储的弹性能。

总弹性能是一个标量，可以写成：
弹性能表达式同样适用于非线性问题。

在这些表达式中，我们假设体力F 为零，并忽略了边界效应。

这些影响可以在以后引入。

积分的意义是每个体积微元的内能总和，其中应力张量单位是Pa ，微元体上的应变εd 没有单位，dV 单位是体积，因此积分出来的单位应该是N·m 。

如果问题是线弹性的，则可以显式的写为：
利用下面的通用公式：
用应变张量ε替换上式中的标量变量x ，弹性张量εc 替换上标量常量a 。

联立上面的式子得到：
我们用c 代替εc 来配合COMSOL Multiphysics 手册中的标记方式。

再提醒一次，如果你不习惯用张量，可以将张量看成是一个3×3的矩阵，点乘是一种张量的运算符号，弹性张量εc 是一个4阶张量（看上去就像4维矩阵）。

更多的标记方法可以参考COMSOL Multiphysics 的Anisotropic Structural Analysis 中的Matrix Notation 。

弹性能积分形式下的单位说明：
3
2
][m N [c]无][m d Pascal =Ω===⎰Ω
单位
ε 最终给出总的积分单位是N·m ――能量。

E W 的表达式就是我们通常说的能量泛函，即位移矢量u （或实际上是u 的梯度）的泛函。

这种函数的函数，而不是坐标的函数，通常被称为泛函，比单元微积分和多元微积分更加抽象。

与积分类似，我们可以说E W 就是函数u 的泛函：
这好比是一个2D 的变量x ，y 的二元函数：
其中x ＝),(y x ，[]a
b
A c d =，[]d b e =。

采用这样的类比是因为在后面我们会看到矩阵A 与有限元的刚度矩阵比较类似。

我们要说明一下函数和泛函的一些区别，古典分析中的函数概念是指两个数集之间所建立的一种对应关系，现代数学的发展却是要求建立两个任意集合之间的某种对应关系。

函数概念被赋予了更为一般的意义，通俗解释泛函指的就是“函数的函数”。

在这里定义域为Ω，泛函可以在整个定义域内进行微分积分等操作。

泛函的变量是函数，这个函数也是有容许空间的。

如果函数u 可以变化，可能会产生一些不符合物理规则的一些现象，例如结构的刚性位移等。

比如一个对u 的基本约束就是材料不能穿越本身。

在有限元分析中，泛函一般是某种能量积分，比如弹性能。

对于其他的物理场，可能是其他的能量积分，或者是一种等效于能量的标量也可以。

至于积分区域，一般由分析对象的CAD 几何区域所确定。

静态电流传导和能量的生成
在静态导电问题中，PDE 方程由最基本的保守形式开始：
0∇⋅=J
其中J 是电流密度。

材料（或本构）模型采用欧姆Ohm 定律：
r σ=B E
其中E 是电场，r σ是电导率。

另外，已知：
V =-∇E
其中V 是静电势，综合以上式子得到
()0r V σ-∇⋅∇=
在COMSOL Multiphysics 中，这就是所谓的Conductive Media DC 方程。

电阻产生的热能
稳态电流的能量问题是在电导体中的电阻热
其中J 表示电流强度，E 代表电场强度，εσ是一个二阶电导张量(3×3)。

如果导体是金属，电导张量一般是一个对角矩阵，如果是晶体，情况就复杂多了。

尽量减少电阻产生的热量，也就是减少热损耗，是我们要研究的一个最小值问题。

如果问题是线性，则积分可以显式地写成：
因为V =-∇E ，其中V 是电势，可以得到：
将这个式子与结构力学中的式子进行对比，发现他们非常相似。

V 的梯度对应于位移梯度，电导率张量εσ对应于弹性张量c 。

在稳态电流和结构力学的计算过程中，张量形式都可以改写为矩阵形式。

传热PDE 方程和能量形式
对于稳态传热问题，PDE 形式为：
其中T 是温度，k 是热传导系数，Q 是空间分布的热源。

热能
基于传热方程的典型泛函为：
其中T 是温度，k 是热传导系数张量（3×3）。

泛函极小值
泛函极值的概念借用了微积分中的不少方法。

本节首先会介绍函数微积分的求极值方法，接下来，我们会借用有限元中常用的术语和标注方法来推导我们熟悉的结果。

这个过程
可以被看作是微积分方法的一种推广。

考虑一个多元微积分函数f，我们要求最小值：
寻找x使得f(x) 最小化
这里x是一个矢量，或者点的坐标。

通过微积分我们知道，这个时候首先必须求函数f的梯度。

将梯度的设置为0，我们可得到一个非线性方程组。

求解方程，我们可以得到一系列的坐标点x，如果在其中某点处的二阶倒数(一般称为Hessian矩阵)为正(或者说有正的特征值)，就说这点就是我们要求的极小点，就好像该点是整个函数的一个谷底一样。

利用Taylor展开的观点，假设已知一个最小值x，我们可以在上面施加一个小的扰动，由Taylor展开可得：
这里H就是前面所说的Hessian矩阵。

现在我们用其他的方法来说明函数f在x最小。

首先，假设x是一个极值点，当添加了一个xδ后，f对于其一阶值不改变。

换句话说，如果我们在x上添加一个xδ来扰动f，其一阶Taylor级数应该为0。

这个条件应该对每个方向都是成立的，否则该点就不是极值点了。

如果上式第二项为0：
对于任意小的xδ都成立，也就是：
我们这里只是用一个稍微有点不同的方法得到了一个同样的结果。

但是，这只是给了我们一个极值点的信息，如果要确定其是最小极值点，必须保证第三项（二阶项）对于任意xδ都为正：
只有当H的特征值都为正时，上式成立（参考线性代数）。

有可能会遇到二阶项也总为0，这个时候我们必须借助更高阶项来判断极值点。

下面是函数f的一个特例：
二次多项式：
其中A是对称矩阵。

如果我们应用Taylor展开，可得到：
或者
这里零阶，一阶和二级项都在独立的中括号内。

为了得到一阶变分，矩阵A必须是对称的。

极值的条件成了：
对于任意小xδ都必须成立，则上式成为：
A T=。

这里我们对矩阵进行了转置，而且利用了矩阵A的对称性，即A
极小值的条件也就是矩阵A必须是一个正定矩阵，如果矩阵A是负定矩阵（只有负特征值），则得到极大值。

如果A是不确定的（特征值有正有负），则极值可能是一个鞍点，既不是极大值，也不是极小值。

如果矩阵A是对称的，而且正定，则函数f是超椭圆的。

在2D中，超椭圆就是椭圆。

二次多项式的几何特征影响经典的PDE方程和有限问题的分类。

当利用有限元方法去离散一个椭圆的PDE问题时候，得到一个对称矩阵（刚度矩阵）的线性代数系统。

这样的问题一般等效于最小能量问题。

弹性静力学问题变分
我们将通过两个步骤来介绍最小能量法理论。

首先粗略说明，让大家熟悉基本概念；接下来考虑细节。

还是以线性静态问题为例，因为这是所有有限元理论都会提到的，从而更容易进行比较。

理论概述
让我们回到线弹性问题的弹性能泛函表达式：
这里的位移矢量u和前面讲的微积分中的点矢量x的角色类似。

W的最小值，我们首先必须得在u上施加一个扰动uδ：要寻找能量泛函
E
上式中两个中间项实质上是一样的（因为c的对称性），所以我们可以写成：
将上式和多元函数表达式对比，我们发现寻找极值点就是找一个使二次项为零的u ：
其中u δ是任意的。

如果我们要寻找的是极小点，则还必须有：
第二项就是泛函的一阶微分：
第三项成为泛函的二级微分：
和前面一样，为了寻找极小点，我们必须保证对于任意δu 第一阶微分为零，二阶微分为正。

这种寻找最小势能函数的方法也可以称作虚功原理。

另外还有一种方法就是初始的时候将扰动写成αδu ，这时对于任意可取的δu ，其能量函数写成()W αδ+u u 。

回到微积分的基本概念，去寻找W 对于α的极值点：
如果我们将它看成是对于α的Taylor 展开，就可以找出其一阶导数（对于极值点必须为零），由于u δ是任意可取的，我们可以得到和前面相同的结果。

小结
上面的过程省略很多推导步骤，如果大家对推导有兴趣，可以试着自己推导。

我们要说
明一下的是：
1、 变量u δ（而不是它的梯度）必须是很小而且是任意的。

2、 这里没有考虑边界条件和体力，比如重力等等。

我们前面所讨论的问题局限于一个没有
任何约束和载荷的边界条件的区域上。

3、 一般来说u δ的限制比多元微积分中x δ宽松。

在泛函中，只要u δ是在容许的范围内即可，
也就是u u δ+必须和物理位移场相对应。

理解这个意思对理解有限元弱形式非常重要。

考虑边界条件和体力
如前面所讲，弹性能的泛函形式是不完整的，因为它没有加上相应的边界条件和载荷。

弹性能的单位是m N ⋅，也就是力乘上位移。

在边界上，我们一般施加面力，或者指定位移，单位为m 。

一般来说，我们希望附加形式是“面力乘上长度”。

同样的方式可以对体力进行处理F 。

在数学上，结构场的边界条件分为两类。

第一类直接定义边界上的力：
其中第一项由定义域内的方程所确定，第二项称为弹簧常数q ，等式右边是面力g 。

这种边界条件就是我们通常说得流量或者Nuemann 边界条件。

第二那边界条件就是
定义一个固定的或者Dirichlet 边界条件。

如果h 是矩阵的形式，r 就是定义了边界上的指定位移。

固定边界条件不能直接加入泛函中去，但是可以通过反力间接加上去。

当指定位移边界时，可以描述一个反力μ（2
/mm N ），也就是弹性体可以在固定处保持不变。

反力就是我们这里用到的Lagrange 乘子，通过添加反力到力作用处的边界，可以忽略到固定边界类型。

这时候我们可以形成统一的边界条件：
这里R 是原始的固定边界，μ是需要计算的反力。

在前面的简化形式中，h 和r 都是常数，所以上式可以变化为：
记住，方程中的每一项都是矢量，表示各个方向的面力。

为了得到所做的功（能量），必须点乘上位移u 。

通过合并一些系数项，将外力写成P ，可简化表达式，这时边界条件可以写成：
P u c n =∇⋅)(
对于其他物理场，可能P 代表边界上的源项。

注意到上式和Navier 方程非常接近：
将能量泛函展开：
关键推导
这个时候，我们又要在u上添加上uδ，可得：
零阶项就是泛函本身，第一阶项是：
这个方程是非常重要的一项。

接下来我们详细介绍。

从前面的讨论可知，我们应该重新组合多项式，保证带有uδ的被积函数成为一项。

如果可以做到，因为uδ是任意的（事实上必须是在容许范围内），我们知道这一项必须为零。

这是我们能找到极值点的唯一方法。

右边第一项需要进一步处理得到我们需要的形式。

第一项我们可以根据Green公式（有时候可能采用的是Stokes原理）进行分部积分：
利用c的对称性，我们可以得到：
利用Green公式得到：
将体积项和边界项合并起来：
确定极值点，必须有：
上式应该对于任何uδ都成立。

因此体积项必须有：
边界项上有：
现在我们又回到PDE问题上了，通过能量最小化原理又重新推回到了PDE形式上！这也是说明最小能量化和PDE形式本质上是统一的一个数学证明。

弱形式
那么，到底什么是弱形式呢？Navier方程的弱形式实际上已经在前面的推导过程中出现过了，即一阶变分的原形式：
如果我们回到COMSOL Multiphysics的文档（或者是关于有限元和弱形式的书籍中），会发现所谓的试函数相当于扰动uδ，
vδ
＝
u
弹性静力学PDE方程的弱形式
为了更好的理解弱形式，我们必须丢弃前面讨论的能量最小化原理，转向一种更加抽象的方法。

弱形式之所以比能量最小化原理更强大，是因为它还可以应用到一些没有得到较好的能量定义的问题中。

首先我们考虑弹性静力学的PDE方程
边界条件是：
抽象的过程如下：乘上容许范围内的试函数v，在感兴趣的域内积分可得：
对左侧利用Green公式进行分部积分：
应用PDE方程的边界条件，可以得到：
整理可得：
这就是PDE方程的弱形式。

如果在积分区域内对于试函数v都是有效的，则上式和PDE 方程是等效的。

PDE方程的解称为强解，而弱形式的解称为弱解。

二者唯一的区别是弱形式对于积分参数的连续性要求比PDE形式低。

由于变形梯度和弹性张量在弱形式里面都不需要微分，所以对函数连续性要求没有那么严格，而在PDE形式中，所有的变量都处在散度的算子下，这要求这些变量必须是可微的。

在弱形式中对于可微的要求放松了（一阶）。

同时，注意到弱形式和前面的一阶变分形式保持了一致，弱形式也可以作为虚功原理的
δ换成了更加抽象的试函数v。

如果弱形式解和能量最一种推广。

只是虚功原理中的位移u
小化原理不一致的时候，极值点变成了鞍点。

也就是，在弱形式中，仍可以将试函数理解为一种推广了的虚位移。

一般性问题的弱形式
正如前面所提到的，弱形式只是PDE方程的一种推广形式，它对变量的连续性要求比较低。

那么能量方法呢？如果有一个定义好了的能量来最小化，那么能量法和弱形式是一致的。

但是，在下列情形下，弱形式更具有适用性：假如PDE方程没有相对应的能量可以进行最小化。

在这种情况下，弱形式仍然是适用的。

由于弱形式对解的要求较低，所以说弱形式比PDE和能量最小化适用范围更广泛。

我们将给出一个没有对应能量最小化的PDE的例子。

对流－扩散PDE问题
对流－扩散PDE问题没有与之相对应的可最小化的能量：
这里c是扩散系数，β是对流系数，α是反应/吸收系数，f是源项。

变量u是标量函数，代表浓度（在COMSOL Multiphysics手册中的Convection-Diffusion模块中，浓度是用变量c表示，扩散系数用D表示）。

在这里我们考虑Neumann边界：
所有困难将集中在刚度K的提取上，主要是对u和Lagrange乘子的线性表达式的集成。

为了得到弱形式，将PDE方程乘以一个试函数v，积分：
这里的试函数v是一个标量函数。

将第一项分部积分，并将所有的项都移到左边，可得到：
加上边界条件，得到：
这就是对流－扩散PDE方程的弱形式。

这个弱形式不能像前面一阶变分那样进行重排。

因为他的对流项，
使得整个系数无法重排。

具体说来，解函数u和试函数v必须在弱形式中的形式保持一致才能和能量泛函的形式保持一致。

但是，在对流项中，u前面带有梯度乘子，而试函数前面却没有任何微分算子的。

没有什么分部积分可以改变这种形式了。

当然，我们也可以看到，实际上弱形式的解和PDE形式的解是保持一致的。

对流项非对称的行为通过数值离散扩展到有限元刚度矩阵上：和能量最小化保持一致的弱形式可以推导出一个对称的刚度矩阵，但是对流－扩散方程推导出来的却是一个非对称的矩阵。

在COMSOL Multiphysics中应用弱形式用户界面的时候，可以输入任意的表达式，包括未知函数u和试函数v的零阶和一阶导数。

你所键入的是弱形式积分中的微分项。

COMSOL Multiphysics的弱形式用法
本章介绍如何在COMSOL Multiphysics中输入弱形式表达式。

对流－扩散PDE问题
假设我们要在COMSOL Multiphysics的用户界面下输入表达式：
约定：COMSOL Multiphysics将所有的项要放在等号右边。

可得到：
区域积分和边界积分可分别在Subdomain Setting和Boundary Setting对话框下设置。

另外，假设我们已经将系数定义为常数或者表达式：
●系数c，P，a和f分别由c，P，a和f表示。

●矢量 的分量由bx，by和bz表示。

在COMSOL Multiphysics中未知函数（因变量）u和试函数v标记如下：
●未知函数u的标记为u
∇的分量标记为ux，uy和uz。

●u
●试函数v的标记为u_test。

∇的分量标记为ux_test，uy_test，uz_test
●v
●只需要输入被积函数，它将被COMSOL Multiphysics自动积分处理。

每一个子域的弱
形式可以有不同的表达式，COMSOL Multiphysics会将各个子域的弱形式整合起来。

输入对流－扩散问题的弱形式：
选择PDE mode下的Weak Form, Subdomain。

在Physics->Subdomain Setting，在Weak Term编辑框中输入：
边界设定，Physics->Boundary Setting，Weak Term编辑框中输入：
COMSOL Multiphysics将边界设置和子域设置分开，因为子域和边界上可以设置不同的数值积分算法。

弱项
如果想要扩展内建的经典PDE模板或者物理应用模式（比如传热），也可以在Physics->Equation System中对应的对话框中输入相同的表达式。

弱形式方程会自动添加在控制方程中。

（通过设置所有的PDE或材料参数为0，选择齐次Neumann边界（流量=0），可以去掉应用模式自动创建的弱形式。

）
Dirichlet或者固定边界，在Boundary setting对话框中的constr编辑框输入弱形式，COMSOL Multiphysics会添加相应的Lagrange乘子（参见用户手册中的边界条件章节）。

结构力学PDE问题
静态结构力学的基本方程是Navier方程：
边界条件：
对流－扩散方程中的标量项现在全部成了矢量和张量，Navier方程的弱形式为：
约定标记如下：
●矢量u的分量：u，v和w。

●位移矢量梯度∇u的分量：ux，uy，uz，vx，vy，vz，wx，wy，wz。

●试位移矢量v的分量：u_test，v_test，w_test。

●试位移矢量梯度∇v的分量：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，ey，...）和应力分量（sx，sy，...）。

后面带有_test后缀的，COMSOL Multiphysics都会和上式一样建立相应的试函数和试函数梯度的表达式。

比如，exy_test等效于0.5*(uy_test+vx_test)。

另一种方式是test()，其中test(xy)表示0.5*(test(uy)+test(vx))，也就相当于0.5*(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和泊松比ν的简单函数。

详情参考COMSOL Multiphysics文档。

在边界上，对于载荷类边界条件，弱形式可以在weak编辑框中写成标量的形式：
Px*u_test+Py*v_test+Pz*w_test
如果采用固定边界，我们必须在其中一个constr编辑框中输入相应的表达式。

对于多物理场仿真，约束和载荷在weak和constr中的形式非常重要，尤其是采用弱约。