有限元的弱形式

PDE弱形式介绍GJ：看到一个介绍COMSOL解决物理问题弱形式的文档，感觉很牛啊，通过COMSOLMultiphysics的弱形式用户界面来求解更多更复杂的问题，这绝对是物理研究的利器啊！而且貌似COMSOL是唯一可以直接使用弱形式来求解问题的软件。为什么要理解PDE方程的弱形式？一般情况下，PDE方程都已经内置在COMSOLMultiphysics的各个模块当中，这种情况下，没有必要去了解PDE方程和及其相关的弱形式。有时候可能问题是没有办法用COMSOLMultiphysics内置模块来求解的，这个时候可以使用经典PDE模版。但是，有时候可能经典PDE模版也不包括要求解的问题，这个时候就只能使用弱形式了（虽然这种情况是极少数的）。另一个原因就是弱形式有时候描述问题比PDE方程紧凑的多。还有，如果你是一个教授去教有限元分析方法，可以帮助学生们直接利用弱形式来更深入的了解有限元。最后，你对有限元方法了解的越多，对于COMSOL中的一些求解器的高级设置就懂得更多。一个重要的事实是：在所有的应用模式和PDE模式求解的时候，COMSOLMultiphysics都是先将方程式系统转为了弱形式，然后进行求解。物理问题的三种描述方式1.偏微分方程2.能量最小化形式3.弱形式PDE问题常常具有最小能量问题的等效形式，这让人有一种直觉，那就是PDE方程都可以有相应的弱形式。实际上这些PDE方程和能量最小值问题只是同一个物理方程的两种不同表达形式罢了，同样，弱形式（几乎）是同一个物理方程的第三个等效形式。我们必须记住，这三种形式只是求解同一个问题的三种不同形式――用数学方法求解真实世界的物理现象。根据不同的需求，这三种方式又有各自不同的优点。三种不同形式的求解PDE形式在各种书籍中比较常见，而且一般都提供了PDE方程的解法。能量法一般见于结构分析的文献中，采用弹性势能最小化形式求解问题是相当自然的一件事。当我们的研究范围超出了标准有限元应用领域，比如传热和结构，这个时候弱形式是不可避免的。化工中的传质问题和流体中的N-S方程都是没有办法用最小能量原理表述出来的。弱形式的特点PDE方程是带有偏微分算子的方程，而能量方程是以积分形式表达的。积分形式的好处就是特别适合于有限元方法，而且不用担心积分变量的不连续，这在偏微分方程中比较普遍。弱形式也是积分形式，拥有和积分形式同样的优点，但是他对积分变量的连续性要求更低，可以看作是能量最小化形式的更一般形式。最重要的是，弱形式非常适合求解非线性的多物理场问题，这就是COMSOLMultiphysics的重点了。PDE到泛函变分GJ：PDE方程一般很难求出解析解，通常需要根据变分原理（数学定律）或最小能量原理（物理定律）转化为泛函变分问题，即得到积分形式，从而便于使用有限元法划分区域离散化，得到刚度矩阵，而最终求解得到PDE的近似数值解。这基本上就是一般的工程中的有限元分析，如平面弹性力学问题、温度场分析及动力学问题等。平面弹性力学问题是通过最小势能原理或虚功原理（两者是同一问题的不同表述形式）建立积分泛函的，温度场可以通过能量法建立泛函，也可以通过变分原理裸建泛函。下面说一说常见的PDE问题根据最小能量原理建立泛函变分。弹性静力学PDE及其弹性能量方程在静力结构分析问题中，我们需要求解的是Navier方程F其中σ是应力张量，F是体力，比如重力等。计算区域记为，其边界记为。应力张量和应变张量之间的关系称为本构关系，线弹性本构一般遵循胡克HOOK定律c其中c是弹性张量，这个关系式说明材料的行为实际上和弹簧差不多（前提是线弹性）。最后，我们可以将应变矢量和位移的关系表述出来u这里u指的是位移矢量u=(u,v,w)，其定义就是变形体上的材料点和未变形时候的位移差。总结以上所有的方程，我们得到了一个二阶PDE方程（Navier方程），()cuF需要一个边界条件来求解，()ncuP其中n是表面的法矢，P是边界上的面力或牵引力。可以顺便提一下，这个PDE方程的弱形式为，其中v=),,(zyxvvv称为试函数。注意，尽管Navier方程是一个矢量表达式，但是上面的表达式是一个标量形式。弹性势能在结构分析中，PDE方程及其弱形式的表达式都不太常见，相反，能量最小化形式因为其直观的表达形式用的较多。这类问题的能量积分形式对应于总势能的最小化，即对象中存储的弹性能。总弹性能是一个标量，可以写成：弹性能表达式同样适用于非线性问题。在这些表达式中，我们假设体力F为零，并忽略了边界效应。这些影响可以在以后引入。积分的意义是每个体积微元的内能总和，其中应力张量单位是Pa，微元体上的应变dε没有单位，dV单位是体积，因此积分出来的单位应该是N·m。如果问题是线弹性的，则可以显式的写为：联立上面的式子得到：我们用c代替c来配合COMSOLMultiphysics手册中的标记方式。弹性能积分形式下的单位说明：32][mN[c]无][mdPascal单位最终给出总的积分单位是N·m――能量。EW的表达式就是我们通常说的能量泛函，即位移矢量u（或实际上是u的梯度）的泛函。这种函数的函数，而不是坐标的函数，通常被称为泛函，比单元微积分和多元微积分更加抽象。与积分类似，我们可以说EW就是函数u的泛函：我们要说明一下函数和泛函的一些区别，古典分析中的函数概念是指两个数集之间所建立的一种对应关系，现代数学的发展却是要求建立两个任意集合之间的某种对应关系。函数概念被赋予了更为一般的意义，通俗解释泛函指的就是“函数的函数”。在这里定义域为，泛函可以在整个定义域内进行微分积分等操作。泛函的变量是函数，这个函数也是有容许空间的。如果函数u可以变化，可能会产生一些不符合物理规则的一些现象，例如结构的刚性位移等。比如一个对u的基本约束就是材料不能穿越本身。在有限元分析中，泛函一般是某种能量积分，比如弹性能。对于其他的物理场，可能是其他的能量积分，或者是一种等效于能量的标量也可以。至于积分区域，一般由分析对象的CAD几何区域所确定。静态电流传导和能量的生成在静态导电问题中，PDE方程由最基本的保守形式开始：0J其中J是电流密度。材料（或本构）模型采用欧姆Ohm定律：rBE其中E是电场，r是电导率。另外，已知：VE其中V是静电势，综合以上式子得到()0rV在COMSOLMultiphysics中，这就是所谓的ConductiveMediaDC方程。电阻产生的热能稳态电流的能量问题是在电导体中的电阻热其中J表示电流强度，E代表电场强度，是一个二阶电导张量(3×3)。如果导体是金属，电导张量一般是一个对角矩阵，如果是晶体，情况就复杂多了。尽量减少电阻产生的热量，也就是减少热损耗，是我们要研究的一个最小值问题。如果问题是线性，则积分可以显式地写成：因为VE，其中V是电势，可以得到：将这个式子与结构力学中的式子进行对比，发现他们非常相似。V的梯度对应于位移梯度，电导率张量对应于弹性张量c。传热PDE方程和能量形式对于稳态传热问题，PDE形式为：其中T是温度，k是热传导系数，Q是空间分布的热源。热能基于传热方程的典型泛函为：其中T是温度，k是热传导系数张量（3×3）。泛函求极值GJ：泛函求极值，即泛函变分，之前写过博客说过它的具体思想，下面的介绍可以说是从另一个角度解释。通过推导会发现，通过能量最小化原理会重新回到了PDE形式上，从而说明能量最小化形式和PDE是同一问题的不同表述。函数求极值考虑一个多元微积分函数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的对称性，即AAT。极小值的条件也就是矩阵A必须是一个正定矩阵，如果矩阵A是负定矩阵（只有负特征值），则得到极大值。如果A是不确定的（特征值有正有负），则极值可能是一个鞍点，既不是极大值，也不是极小值。如果矩阵A是对称的，而且正定，则函数f是超椭圆的。在2D中，超椭圆就是椭圆。二次多项式的几何特征影响经典的PDE方程和有限问题的分类。当利用有限元方法去离散一个椭圆的PDE问题时候，得到一个对称矩阵（刚度矩阵）的线性代数系统。这样的问题一般等效于最小能量问题。弹性静力学问题泛函求极值还是以线性静态问题为例，因为这是所有有限元理论都会提到的，从而更容易进行比较。理论概述让我们回到线弹性问题的弹性能泛函表达式：这里的位移矢量u和前面讲的微积分中的点矢量x的角色类似。要寻找能量泛函EW的最小值，我们首先必须得在u上施加一个扰动u：上式中两个中间项实质上是一样的（因为c的对称性），所以我们可以写成：将上式和多元函数表达式对比，我们发现寻找极值点就是找一个使二次项为零的u：其中u是任意的。如果我们要寻找的是极小点，则还必须有：第二项就是泛函的一阶微分：第三项成为泛函的二级微分：和前面一样，为了寻找极小点，我们必须保证对于任意u第一阶微分为零，二阶微分为正。这种寻找最小势能函数的方法也可以称作虚功原理。（这就明白理论力学里所谓虚位移的意义了，就是一个任意的扰动）另外还有一种方法就是初始的时候将扰动写成u，这时对于任意可取的u，其能量函数写成()Wuu。回到微积分的基本概念，去寻找W对于的极值点：如果我们将它看成是对于的Taylor展开，就可以找出其一阶导数（对于极值点必须为零），由于u是任意可取的，我们可以得到和前面相同的结果。小结上面的过程省略很多推导步骤，如果大家对推导有兴趣，可以试着自己推导。要说明一下的是：1、变量u（而不是它的梯度）必须是很小而且是任意的。2、这里没有考虑边界条件和体力，比如重力等等。我们前面所讨论的问题局限于一个没有任何约束和载荷的边界条件的区域上。3、一般来说u的限制比多元微积分中x宽松。在泛函中，只要u是在容许的范围内即可，也就是uu必须和物理位移场相对应。理解这个意思对理解有限元弱形式非常重要。考虑边界条件和体力如前面所讲，弹性能的泛函形式是不完整的，因为它没有加上相应的边界条件和载荷。弹性能的单位是mN，也就是力乘上位移。在边界上，我们一般施加面力，或者指定位移，单位为m。一般来说，我们希望附加形式是“面力乘上长度”。同样的方式可以对体力进行处理F。在数学上，结构场的边界条件分为两类。第一类直接定义边界上的力：其中第一项由定义域内的方程所确定，第二项称为弹簧常数q，等式右边是面力g。这种边界条件就是我们通常说得流量或者Nuemann