PDE 弱形式介绍
GJ :看到一个介绍COMSOL 解决物理问题弱形式的文档,感觉很牛啊,通过COMSOL Multiphysics 的弱形式用户界面来求解更多更复杂的问题,这绝对是物理研究的利器啊!而且貌似COMSOL 是唯一可以直接使用弱形式来求解问题的软件。
为什么要理解PDE 方程的弱形式?
一般情况下,PDE 方程都已经内置在COMSOL Multiphysics 的各个模块当中,这种情况下,没有必要去了解PDE 方程和及其相关的弱形式。有时候可能问题是没有办法用COMSOL Multiphysics 内置模块来求解的,这个时候可以使用经典PDE 模版。但是,有时候可能经典PDE 模版也不包括要求解的问题,这个时候就只能使用弱形式了(虽然这种情况是极少数的)。另一个原因就是弱形式有时候描述问题比PDE 方程紧凑的多。还有,如果你是一个教授去教有限元分析方法,可以帮助学生们直接利用弱形式来更深入的了解有限元。最后,你对有限元方法了解的越多,对于COMSOL 中的一些求解器的高级设置就懂得更多。 一个重要的事实是:在所有的应用模式和PDE 模式求解的时候,COMSOL Multiphysics 都是先将方程式系统转为了弱形式,然后进行求解。
物理问题的三种描述方式
1. 偏微分方程
2. 能量最小化形式
3. 弱形式
PDE 问题常常具有最小能量问题的等效形式,这让人有一种直觉,那就是PDE 方程都可以有相应的弱形式。实际上这些PDE 方程和能量最小值问题只是同一个物理方程的两种不同表达形式罢了,同样,弱形式(几乎)是同一个物理方程的第三个等效形式。我们必须记住,这三种形式只是求解同一个问题的三种不同形式――用数学方法求解真实世界的物理现象。根据不同的需求,这三种方式又有各自不同的优点。
三种不同形式的求解
PDE 形式在各种书籍中比较常见,而且一般都提供了PDE 方程的解法。能量法一般见于结构分析的文献中,采用弹性势能最小化形式求解问题是相当自然的一件事。当我们的研究范围超出了标准有限元应用领域,比如传热和结构,这个时候弱形式是不可避免的。化工中的传质问题和流体中的N-S 方程都是没有办法用最小能量原理表述出来的。
弱形式的特点
PDE 方程是带有偏微分算子的方程,而能量方程是以积分形式表达的。积分形式的好处就是特别适合于有限元方法,而且不用担心积分变量的不连续,这在偏微分方程中比较普遍。弱形式也是积分形式,拥有和积分形式同样的优点,但是他对积分变量的连续性要求更低,可以看作是能量最小化形式的更一般形式。最重要的是,弱形式非常适合求解非线性的多物理场问题,这就是COMSOL Multiphysics的重点了。
PDE 到泛函变分
GJ :PDE 方程一般很难求出解析解,通常需要根据变分原理(数学定律)或最小能量原理(物理定律)转化为泛函变分问题,即得到积分形式,从而便于使用有限元法划分区域离散化,得到刚度矩阵,而最终求解得到PDE 的近似数值解。这基本上就是一般的工程中的有限元分析,如平面弹性力学问题、温度场分析及动力学问题等。平面弹性力学问题是通过最小势能原理或虚功原理(两者是同一问题的不同表述形式)建立积分泛函的,温度场可以通过能量法建立泛函,也可以通过变分原理裸建泛函。
下面说一说常见的PDE 问题根据最小能量原理建立泛函变分。
弹性静力学PDE 及其弹性能量方程
在静力结构分析问题中,我们需要求解的是Navier 方程
-∇⋅σ=F
其中ζ是应力张量,F 是体力,比如重力等。计算区域记为Ω,其边界记为∂Ω。
应力张量σ和应变张量ε之间的关系称为本构关系,线弹性本构一般遵循胡克HOOK 定律
σ=c εε
其中c ε是弹性张量,这个关系式说明材料的行为实际上和弹簧差不多(前提是线弹性)。 最后,我们可以将应变矢量和位移的关系表述出来
ε=∇u
这里u 指的是位移矢量u =(u,v,w),其定义就是变形体上的材料点和未变形时候的位移差。 总结以上所有的方程,我们得到了一个二阶PDE 方程(Navier 方程),
-∇⋅(c ∇u ) =F
需要一个边界条件来求解,
n ⋅(c ∇u ) =P
其中n 是表面∂Ω的法矢,P 是边界上的面力或牵引力。
可以顺便提一下,这个PDE 方程的弱形式为,
其中v=(v x , v y , v z ) 称为试函数。注意,尽管Navier 方程是一个矢量表达式,但是上面的表达式是一个标量形式。
弹性势能
在结构分析中,PDE 方程及其弱形式的表达式都不太常见,相反,能量最小化形式因为其直观的表达形式用的较多。这类问题的能量积分形式对应于总势能的最小化,即对象中存储的弹性能。
总弹性能是一个标量,可以写成:
弹性能表达式同样适用于非线性问题。在这些表达式中,我们假设体力F 为零,并忽略了边界效应。这些影响可以在以后引入。积分的意义是每个体积微元的内能总和,其中应力张量单位是Pa ,微元体上的应变d ε没有单位,dV 单位是体积,因此积分出来的单位应该是N·m 。
如果问题是线弹性的,则可以显式的写为:
联立上面的式子得到:
我们用c 代替c ε来配合COMSOL Multiphysics手册中的标记方式。
弹性能积分形式下的单位说明:
[ε]=无单位
N =Pascal 2m
[⎰d Ω]=m 3[c]=
Ω
最终给出总的积分单位是N·m ――能量。
W E 的表达式就是我们通常说的能量泛函,即位移矢量u (或实际上是u 的梯度)的泛函。这种函数的函数,而不是坐标的函数,通常被称为泛函,比单元微积分和多元微积分更加抽象。
与积分类似,我们可以说W E 就是函数u 的泛函:
我们要说明一下函数和泛函的一些区别,古典分析中的函数概念是指两个数集之间所建立的一种对应关系,现代数学的发展却是要求建立两个任意集合之间的某种对应关系。函数概念被赋予了更为一般的意义,通俗解释泛函指的就是“函数的函数”。在这里定义域为Ω,泛函可以在整个定义域内进行微分积分等操作。
泛函的变量是函数,这个函数也是有容许空间的。如果函数u 可以变化,可能会产生一些不符合物理规则的一些现象,例如结构的刚性位移等。比如一个对u 的基本约束就是材料不能穿越本身。
在有限元分析中,泛函一般是某种能量积分,比如弹性能。对于其他的物理场,可能是其他的能量积分,或者是一种等效于能量的标量也可以。至于积分区域,一般由分析对象的CAD 几何区域所确定。
静态电流传导和能量的生成
在静态导电问题中,PDE 方程由最基本的保守形式开始:
∇⋅J =0
其中J 是电流密度。
材料(或本构)模型采用欧姆Ohm 定律:
B =σr E
其中E 是电场,σr 是电导率。
另外,已知:
E =-∇V
其中V 是静电势,综合以上式子得到
-∇⋅(σr ∇V ) =0
在COMSOL Multiphysics中,这就是所谓的Conductive Media DC方程。
电阻产生的热能
稳态电流的能量问题是在电导体中的电阻热
其中J 表示电流强度,E 代表电场强度,σε是一个二阶电导张量(3×3) 。如果导体是金属,电导张量一般是一个对角矩阵,如果是晶体,情况就复杂多了。
尽量减少电阻产生的热量,也就是减少热损耗,是我们要研究的一个最小值问题。 如果问题是线性,则积分可以显式地写成:
因为E =-∇V ,其中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 的对称性,即A =A 。 T
极小值的条件也就是矩阵A 必须是一个正定矩阵,如果矩阵A 是负定矩阵(只有负特征值),则得到极大值。如果A 是不确定的(特征值有正有负),则极值可能是一个鞍点,既不是极大值,也不是极小值。如果矩阵A 是对称的,而且正定,则函数f 是超椭圆的。在2D 中,超椭圆就是椭圆。二次多项式的几何特征影响经典的PDE 方程和有限问题的分类。当利用有限元方法去离散一个椭圆的PDE 问题时候,得到一个对称矩阵(刚度矩阵)的线性代数系统。这样的问题一般等效于最小能量问题。
弹性静力学问题泛函求极值
还是以线性静态问题为例,因为这是所有有限元理论都会提到的,从而更容易进行比较。 理论概述
让我们回到线弹性问题的弹性能泛函表达式:
这里的位移矢量u 和前面讲的微积分中的点矢量x 的角色类似。
要寻找能量泛函W E 的最小值,我们首先必须得在u 上施加一个扰动δu :
上式中两个中间项实质上是一样的(因为c 的对称性),所以我们可以写成:
将上式和多元函数表达式对比,我们发现寻找极值点就是找一个使二次项为零的u :
其中δu 是任意的。
如果我们要寻找的是极小点,则还必须有:
第二项就是泛函的一阶微分:
第三项成为泛函的二级微分:
和前面一样,为了寻找极小点,我们必须保证对于任意δu 第一阶微分为零,二阶微分为正。这种寻找最小势能函数的方法也可以称作虚功原理。(这就明白理论力学里所谓虚位移的意义了,就是一个任意的扰动)
另外还有一种方法就是初始的时候将扰动写成αδu ,这时对于任意可取的δu ,其能量函数写成W (u +αδu ) 。回到微积分的基本概念,去寻找W 对于α的极值点:
如果我们将它看成是对于α的Taylor 展开,就可以找出其一阶导数(对于极值点必须为零),由于δu 是任意可取的,我们可以得到和前面相同的结果。
小结
上面的过程省略很多推导步骤,如果大家对推导有兴趣,可以试着自己推导。要说明一下的是:
1、 变量δu (而不是它的梯度)必须是很小而且是任意的。
2、 这里没有考虑边界条件和体力,比如重力等等。我们前面所讨论的问题局限于一个没有
任何约束和载荷的边界条件的区域上。
3、 一般来说δu 的限制比多元微积分中δx 宽松。在泛函中,只要δu 是在容许的范围内即可,
也就是u +δu 必须和物理位移场相对应。理解这个意思对理解有限元弱形式非常重要。
考虑边界条件和体力
如前面所讲,弹性能的泛函形式是不完整的,因为它没有加上相应的边界条件和载荷。弹性能的单位是N ⋅m ,也就是力乘上位移。在边界上,我们一般施加面力,或者指定位移,单位为m 。一般来说,我们希望附加形式是“面力乘上长度”。同样的方式可以对体力进行处理F 。
在数学上,结构场的边界条件分为两类。第一类直接定义边界上的力:
其中第一项由定义域内的方程所确定,第二项称为弹簧常数q ,等式右边是面力g 。这种边界条件就是我们通常说得流量或者Nuemann 边界条件。
第二类边界条件就是
定义一个固定的或者Dirichlet 边界条件。如果h 是矩阵的形式,r 就是定义了边界上的指定位移。固定边界条件不能直接加入泛函中去,但是可以通过反力间接加上去。当指定位移边界时,可以描述一个反力μ(N /mm ),也就是弹性体可以在固定处保持不变。反力就是我们这里用到的Lagrange 乘子,通过添加反力到力作用处的边界,可以忽略到固定边界类型。这时候我们可以形成统一的边界条件:
2
这里R 是原始的固定边界,μ是需要计算的反力。在前面的简化形式中,h 和r 都是常数,所以上式可以变化为:
记住,方程中的每一项都是矢量,表示各个方向的面力。为了得到所做的功(能量),必须点乘上位移u 。
通过合并一些系数项,将外力写成P ,可简化表达式,这时边界条件可以写成:
n ⋅(c ∇u ) =P
对于其他物理场,可能P 代表边界上的源项。 注意到上式和Navier 方程非常接近:
将能量泛函展开:
关键推导
这个时候,我们又要在u 上添加上δu ,可得:
零阶项就是泛函本身,第一阶项是:
这个方程是非常重要的一项。
从前面的讨论可知,我们应该重新组合多项式,保证带有δu 的被积函数成为一项。如果可以做到,因为δu 是任意的(事实上必须是在容许范围内),我们知道这一项必须为零。这是我们能找到极值点的唯一方法。右边第一项需要进一步处理得到我们需要的形式。 第一项我们可以根据Green 公式(有时候可能采用的是Stokes 原理)进行分部积分:
利用c 的对称性,我们可以得到:
利用Green 公式得到:
将体积项和边界项合并起来:
确定极值点,必须有:
上式应该对于任何δu 都成立。因此体积项必须有:
边界项上有:
现在我们又回到PDE 问题上了,这说明泛函的理论解就是PDE 方程的解,即通过能量最小化原理又重新推回到了PDE 形式上!这也是说明最小能量化和PDE 形式本质上是统一的一个数学证明。
GJ :要注意的是,实际有限元的求解不是从泛函又导回PDE 方程,而是通过网格划分离散化,得到数值近似解。而通过建立的泛函求解数值近似解应该算是比较完备的方法,这里说完备,是和弱形式对比来的,两者的区别后面会说。
弱形式
GJ :上面是利用能量最小化形式或者变分原理建立的泛函,即PDE 的等效积分形式,下面说的是弱形式建立。
那么,到底什么是弱形式呢?Navier 方程的弱形式实际上已经在前面的推导过程中出现过了,即一阶变分的原形式:
如果我们回到COMSOL Multiphysics的文档(或者是关于有限元和弱形式的书籍中),会发现所谓的试函数相当于扰动δu ,
v =δu
弹性静力学PDE 方程的弱形式
为了更好的理解弱形式,我们必须丢弃前面讨论的能量最小化原理,转向一种更加抽象的方法。弱形式之所以比能量最小化原理更强大,是因为它还可以应用到一些没有得到较好的能量定义的问题中。
首先我们考虑弹性静力学的PDE 方程
边界条件是:
抽象的过程如下:乘上容许范围内的试函数v ,在感兴趣的域内积分可得:
对左侧利用Green 公式进行分部积分:
应用PDE 方程的边界条件,可以得到:
整理可得:
这就是PDE 方程的弱形式。如果在积分区域内对于试函数v 都是有效的,则上式和PDE 方程是等效的。PDE 方程的解称为强解,而弱形式的解称为弱解。二者唯一的区别是弱形式对于积分参数的连续性要求比PDE 形式低。由于变形梯度和弹性张量在弱形式里面都不需要微分,所以对函数连续性要求没有那么严格,而在PDE 形式中,所有的变量都处在散度的算子下,这要求这些变量必须是可微的。在弱形式中对于可微的要求放松了(一阶)。 同时,注意到弱形式和前面的一阶变分形式保持了一致,弱形式也可以作为虚功原理的一种推广。只是虚功原理中的位移δu 换成了更加抽象的试函数v 。如果弱形式解和能量最小化原理不一致的时候,极值点变成了鞍点。也就是,在弱形式中,仍可以将试函数理解为一种推广了的虚位移。
一般性问题的弱形式
正如前面所提到的,弱形式只是PDE 方程的一种推广形式,它对变量的连续性要求比较低。那么能量方法呢?如果有一个定义好了的能量来最小化,那么能量法和弱形式是一致的。但是,在下列情形下,弱形式更具有适用性:假如PDE 方程没有相对应的能量可以进行最小化。在这种情况下,弱形式仍然是适用的。由于弱形式对解的要求较低,所以说弱形式比PDE 和能量最小化适用范围更广泛。
GJ :弱形式和最小能量形式的区别就在于虚位移δu 与试函数v 的差别,如下面两式
也就是说,泛函求极值即为泛函的变分为0,如上面的式子1,所以泛函的有限元解对任意扰动δu 成立,而从式子2可以看到,弱形式的解只是对自己设定的试函数v 成立。所以泛函求极值得到近似函数是弱形式的特殊形式,即弱形式的试函数v 可以任意取而求得的近似函数,所以从这种意义上说泛函形式求得的近似解更完备。
但很多情况下无法得到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
● ∇u 的分量标记为ux ,uy 和uz 。
● 试函数v 的标记为u_test。
● ∇v 的分量标记为ux_test,uy_test,uz_test
● 只需要输入被积函数,它将被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 中的形式非常重要,尤其是采用弱约束的时候。更多详情可以参考COMSOL Multiphysics文档以及和Lagrange 乘子相关的技术文档。
尽管弱形式是一个标量表达式,但是COMSOL Multiphysics中,弱形式有和PDE 系统一样多的未知量需要文本输入。原因在于不同的多物理场问题可能需要不同的有限元分析类型和保证其数值稳定型的积分算法。对于3D 结构分析,弱形式中有三个文本输入框。但是,在离散之前,采用了同样的有限元单元和积分类型进行合并,这样就可以选择不同的弱形式进行操作。例如你可以在第一个域内选择弱形式,而其他的域内设置为空白。
对于流动问题的Navier-stokes 方程,情况又稍微有些不同。和未知的速度场相比,未知的压力采用一个低阶有限元来离散。这种情况下,不能将所有的弱项全部在同一个弱域内输入。为了保证数值稳定型,必须依靠混阶有限元(mixed finite element) 。混阶有限元并不是COMSOL Multiphysics特别制定,而是数值算法所需要的。
有限元方法
本章说明弱形式如何利用有限元方法来进行离散。假设我们需要离散以下扩散问题:
这是一个对流-扩散问题的特殊情况,其中α=0,β=0。
有限元的基本实现是将整个计算域Ω离散为多个特别简单的形状的小单元,比如2D 中的三角形,3D 中的四面体等等。相应的网格,例如三角形,由边和节点组成。下一步就是要选择一个比较容易实现的一些近似方法,其中一种比较简单的方法就是将解表示为采用线性多项式插值的所谓基函数的和。基函数的构造方法是指定某个节点为1,而相邻的节点为0,二者之间的值就是从0到1线性变化。这里说的相邻指的是中间有一条边将其连接起来。
遍历三角形网格的所有节点(从1到N )。定义节点i 的基函数为ϕi ,也就是在节点i 处其值为1,其他点处值为0。注意ϕi 只是在节点i 及其相邻的三角形内不为零。现在假设真实值u 可以用基函数的求和u h 来近似描述:
参数U i 是u h 在节点i 的值。同样,我们可以对试函数进行类似处理:
下标h 表示离散函数属于由所有三角形边中最长边表示的具有确定的网格尺寸h 的网格。
由于我们可以任意选择试函数,因此可以将除了j 点以外的所有的V j 设置为零,
接下来我们将所有的试函数(j =1,...,N )输入到弱形式中去,每个试函数都可以得到一
个方程。这样可以生成一个线性代数系统,系统矩阵就是我们所说的刚度矩阵。为什么我们可以自由选择试函数,不妨回想一下前面提到的弱形式需要对所有可取的试函数成立。选择试函数是有限元方法的重要环节,因为他在很大程度上影响着刚度矩阵。由于刚度矩阵中很多项为零,所以一般是稀疏矩阵。
当我们使用试函数的时候,生成的有限元刚度矩阵应该是一个方阵。如果弱形式本身定义良好的话,刚度矩阵应该是非奇异的,也就是说系统有一个唯一解。
现在考虑扩散方程的弱形式:
将表达式写成离散形式:
方程重新排列:
采用矩阵标注可得:
或者:
在这里刚度矩阵K 是:
解矢量U 的单元为U i ,载荷矢量L 的单元为,
现在我们明白为什么选择基函数和试函数很关键了。如果我们关注刚度矩阵K ij ,会发现其中很多元素为零,因为前面已经提到每个ϕj 都是大部分为零,同样ϕj 的梯度也是大部分为零的。有很多有效的算法去求解这类稀疏矩阵,COMSOL Multiphysics 提供一套稀疏线性系统求解器。
有限元方法同样适用于非线性问题。非线性方法一般来说采用迭代的算法,每一次迭代就是求解一个与上面类似的线性弱形式方法。
抽象和几何解释
为什么有限元可以解决问题,它是如何解决问题的?前面的讨论中可以找到一些答案。为了有一个更清晰的答案,我们需要了解一些更多的泛函的概念。我们将发现有限元方法通过一种优化方法将解投影到一个有限维函数空间来求解。
标量积
为了得到有限元方法的几何解释,或脑海中的意象,我们需要熟悉标量积的概念。在线性代数中,我们知道两个矢量f 和g 的标量积为
这里的矢量f 和g 属于3维矢量空间。标量积可以推广到任意维N ,
标量积有时也称做内积。
如果两个矢量正交,
一个矢量f 可以通过标量积投影到另一个矢量
g
其中f p 是与g 平等的投影矢量:
矢量差
与g 正交:
投影矢量的唯一性特征是原始和投影矢量之间的差与它所投影的矢量正交。如果需要找到在方向g 上与f 最近的矢量,f p 就是我们的答案。矢量e 可以看作是关于f 到f p 之间的近似误差。换句话说,误差矢量e 与矢量g 正交。后面在对有限元进行几何解释时将用到这个结论。但首先我们得介绍一些泛函分析的概念。
与矢量不同,泛函分析讲的是函数,它们必须属于无限维的矢量空间。我们可以积分形式定义一个无限维矢量空间(函数空间)中的标量积:
如果两个函数正交,则有
进一步将标量积的概念推广,并且考虑包含函数梯度的被积函数,
下标1和2用来说明上面是两种不同的标量积。
Hilbert 空间
如果对于一个标量积,我们只考虑那些与自身进行的标量积(积分)有一个有限值的函数u ,
我们说这些函数属于一个确定的函数类,或函数空间。由标量积和定义域Ω组成的函数空间就被称为Hilbert 空间。
对应于上面的标量积1的Hilbert 空间通常标记为L 2,与标量积2对应的Hilbert 空间常称为H 1。通常空间L 2中的函数比H 1中的多,因为H 1中的函数自动地存在于L 2中,反之则不一定。我们也可以将这种现象称为H 1是L 2的一个子空间,或
有限元方法的抽象形式
现在让我们考虑PDE 问题[3]:
在域Ω中,边界条件为u =0。这是对流-扩散方程的一种特例,其中c =a =1,β=0。通过分部积分,表明当试函数选择成在边界上具有相同的边界条件v =0时,弱形式中的边界项为0。得到弱形式为:
其中包含前面提到的两种标量积,可改写为:
找到u ,使得对于所有属于相配的Hilbert 空间中的v ,有:
这里的相配的空间是H 1,且在边界上v =0。解函数u 也必须属于这个空间。注意,对于前面提到的弱形式,我们可以很自由地添加不同的标量积,因为每个积的结果是实数或虚数。 现在我们选择前面讨论的基函数的和(线性组合)来近似u 和v ,近似解被称作u h 和v h 。这些函数属于Hilbert 空间,可称为V h ,由线性基函数ϕi 扩展而来。V h 的空间维数为N (基函数的数目)。此外,V h 是H 1的子空间,
换句话说,如果函数属于V h ,则它自动地属于H 1。空间H 1是一个大得多的空间,因为它是无限维的。有限的基函数不可能扩散成属于H 1的所有函数。
弱形式的有限元现在变成了:
对于V h 中的所有v h ,在V h 中找到u h ,使得
现在我们对有限元方法在函数空间的作为一个确定的投影的几何解释有了更深入的了解。最后,在原始的弱形式中:
用v h 代替v 。这是合理的,因为v 在H 1中,而v h 在V h 中,由于V h 是H 1的子空间,因此v h 也是在H 1中,
现在从弱形式中减去有限元解,得到:
或
即离散误差:
与所有的V h 中的v h 正交。也就是说,有限元解u h 是真实解u 在属于H 1的有限维子空间V h
的投影。最终,我们得到了关于有限元方法的几何解释,
对于给定的网格,有限元解u h 是在函数空间V h 中关于标量积(,) 2最接近真实解u 的解。
PDE 弱形式介绍
GJ :看到一个介绍COMSOL 解决物理问题弱形式的文档,感觉很牛啊,通过COMSOL Multiphysics 的弱形式用户界面来求解更多更复杂的问题,这绝对是物理研究的利器啊!而且貌似COMSOL 是唯一可以直接使用弱形式来求解问题的软件。
为什么要理解PDE 方程的弱形式?
一般情况下,PDE 方程都已经内置在COMSOL Multiphysics 的各个模块当中,这种情况下,没有必要去了解PDE 方程和及其相关的弱形式。有时候可能问题是没有办法用COMSOL Multiphysics 内置模块来求解的,这个时候可以使用经典PDE 模版。但是,有时候可能经典PDE 模版也不包括要求解的问题,这个时候就只能使用弱形式了(虽然这种情况是极少数的)。另一个原因就是弱形式有时候描述问题比PDE 方程紧凑的多。还有,如果你是一个教授去教有限元分析方法,可以帮助学生们直接利用弱形式来更深入的了解有限元。最后,你对有限元方法了解的越多,对于COMSOL 中的一些求解器的高级设置就懂得更多。 一个重要的事实是:在所有的应用模式和PDE 模式求解的时候,COMSOL Multiphysics 都是先将方程式系统转为了弱形式,然后进行求解。
物理问题的三种描述方式
1. 偏微分方程
2. 能量最小化形式
3. 弱形式
PDE 问题常常具有最小能量问题的等效形式,这让人有一种直觉,那就是PDE 方程都可以有相应的弱形式。实际上这些PDE 方程和能量最小值问题只是同一个物理方程的两种不同表达形式罢了,同样,弱形式(几乎)是同一个物理方程的第三个等效形式。我们必须记住,这三种形式只是求解同一个问题的三种不同形式――用数学方法求解真实世界的物理现象。根据不同的需求,这三种方式又有各自不同的优点。
三种不同形式的求解
PDE 形式在各种书籍中比较常见,而且一般都提供了PDE 方程的解法。能量法一般见于结构分析的文献中,采用弹性势能最小化形式求解问题是相当自然的一件事。当我们的研究范围超出了标准有限元应用领域,比如传热和结构,这个时候弱形式是不可避免的。化工中的传质问题和流体中的N-S 方程都是没有办法用最小能量原理表述出来的。
弱形式的特点
PDE 方程是带有偏微分算子的方程,而能量方程是以积分形式表达的。积分形式的好处就是特别适合于有限元方法,而且不用担心积分变量的不连续,这在偏微分方程中比较普遍。弱形式也是积分形式,拥有和积分形式同样的优点,但是他对积分变量的连续性要求更低,可以看作是能量最小化形式的更一般形式。最重要的是,弱形式非常适合求解非线性的多物理场问题,这就是COMSOL Multiphysics的重点了。
PDE 到泛函变分
GJ :PDE 方程一般很难求出解析解,通常需要根据变分原理(数学定律)或最小能量原理(物理定律)转化为泛函变分问题,即得到积分形式,从而便于使用有限元法划分区域离散化,得到刚度矩阵,而最终求解得到PDE 的近似数值解。这基本上就是一般的工程中的有限元分析,如平面弹性力学问题、温度场分析及动力学问题等。平面弹性力学问题是通过最小势能原理或虚功原理(两者是同一问题的不同表述形式)建立积分泛函的,温度场可以通过能量法建立泛函,也可以通过变分原理裸建泛函。
下面说一说常见的PDE 问题根据最小能量原理建立泛函变分。
弹性静力学PDE 及其弹性能量方程
在静力结构分析问题中,我们需要求解的是Navier 方程
-∇⋅σ=F
其中ζ是应力张量,F 是体力,比如重力等。计算区域记为Ω,其边界记为∂Ω。
应力张量σ和应变张量ε之间的关系称为本构关系,线弹性本构一般遵循胡克HOOK 定律
σ=c εε
其中c ε是弹性张量,这个关系式说明材料的行为实际上和弹簧差不多(前提是线弹性)。 最后,我们可以将应变矢量和位移的关系表述出来
ε=∇u
这里u 指的是位移矢量u =(u,v,w),其定义就是变形体上的材料点和未变形时候的位移差。 总结以上所有的方程,我们得到了一个二阶PDE 方程(Navier 方程),
-∇⋅(c ∇u ) =F
需要一个边界条件来求解,
n ⋅(c ∇u ) =P
其中n 是表面∂Ω的法矢,P 是边界上的面力或牵引力。
可以顺便提一下,这个PDE 方程的弱形式为,
其中v=(v x , v y , v z ) 称为试函数。注意,尽管Navier 方程是一个矢量表达式,但是上面的表达式是一个标量形式。
弹性势能
在结构分析中,PDE 方程及其弱形式的表达式都不太常见,相反,能量最小化形式因为其直观的表达形式用的较多。这类问题的能量积分形式对应于总势能的最小化,即对象中存储的弹性能。
总弹性能是一个标量,可以写成:
弹性能表达式同样适用于非线性问题。在这些表达式中,我们假设体力F 为零,并忽略了边界效应。这些影响可以在以后引入。积分的意义是每个体积微元的内能总和,其中应力张量单位是Pa ,微元体上的应变d ε没有单位,dV 单位是体积,因此积分出来的单位应该是N·m 。
如果问题是线弹性的,则可以显式的写为:
联立上面的式子得到:
我们用c 代替c ε来配合COMSOL Multiphysics手册中的标记方式。
弹性能积分形式下的单位说明:
[ε]=无单位
N =Pascal 2m
[⎰d Ω]=m 3[c]=
Ω
最终给出总的积分单位是N·m ――能量。
W E 的表达式就是我们通常说的能量泛函,即位移矢量u (或实际上是u 的梯度)的泛函。这种函数的函数,而不是坐标的函数,通常被称为泛函,比单元微积分和多元微积分更加抽象。
与积分类似,我们可以说W E 就是函数u 的泛函:
我们要说明一下函数和泛函的一些区别,古典分析中的函数概念是指两个数集之间所建立的一种对应关系,现代数学的发展却是要求建立两个任意集合之间的某种对应关系。函数概念被赋予了更为一般的意义,通俗解释泛函指的就是“函数的函数”。在这里定义域为Ω,泛函可以在整个定义域内进行微分积分等操作。
泛函的变量是函数,这个函数也是有容许空间的。如果函数u 可以变化,可能会产生一些不符合物理规则的一些现象,例如结构的刚性位移等。比如一个对u 的基本约束就是材料不能穿越本身。
在有限元分析中,泛函一般是某种能量积分,比如弹性能。对于其他的物理场,可能是其他的能量积分,或者是一种等效于能量的标量也可以。至于积分区域,一般由分析对象的CAD 几何区域所确定。
静态电流传导和能量的生成
在静态导电问题中,PDE 方程由最基本的保守形式开始:
∇⋅J =0
其中J 是电流密度。
材料(或本构)模型采用欧姆Ohm 定律:
B =σr E
其中E 是电场,σr 是电导率。
另外,已知:
E =-∇V
其中V 是静电势,综合以上式子得到
-∇⋅(σr ∇V ) =0
在COMSOL Multiphysics中,这就是所谓的Conductive Media DC方程。
电阻产生的热能
稳态电流的能量问题是在电导体中的电阻热
其中J 表示电流强度,E 代表电场强度,σε是一个二阶电导张量(3×3) 。如果导体是金属,电导张量一般是一个对角矩阵,如果是晶体,情况就复杂多了。
尽量减少电阻产生的热量,也就是减少热损耗,是我们要研究的一个最小值问题。 如果问题是线性,则积分可以显式地写成:
因为E =-∇V ,其中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 的对称性,即A =A 。 T
极小值的条件也就是矩阵A 必须是一个正定矩阵,如果矩阵A 是负定矩阵(只有负特征值),则得到极大值。如果A 是不确定的(特征值有正有负),则极值可能是一个鞍点,既不是极大值,也不是极小值。如果矩阵A 是对称的,而且正定,则函数f 是超椭圆的。在2D 中,超椭圆就是椭圆。二次多项式的几何特征影响经典的PDE 方程和有限问题的分类。当利用有限元方法去离散一个椭圆的PDE 问题时候,得到一个对称矩阵(刚度矩阵)的线性代数系统。这样的问题一般等效于最小能量问题。
弹性静力学问题泛函求极值
还是以线性静态问题为例,因为这是所有有限元理论都会提到的,从而更容易进行比较。 理论概述
让我们回到线弹性问题的弹性能泛函表达式:
这里的位移矢量u 和前面讲的微积分中的点矢量x 的角色类似。
要寻找能量泛函W E 的最小值,我们首先必须得在u 上施加一个扰动δu :
上式中两个中间项实质上是一样的(因为c 的对称性),所以我们可以写成:
将上式和多元函数表达式对比,我们发现寻找极值点就是找一个使二次项为零的u :
其中δu 是任意的。
如果我们要寻找的是极小点,则还必须有:
第二项就是泛函的一阶微分:
第三项成为泛函的二级微分:
和前面一样,为了寻找极小点,我们必须保证对于任意δu 第一阶微分为零,二阶微分为正。这种寻找最小势能函数的方法也可以称作虚功原理。(这就明白理论力学里所谓虚位移的意义了,就是一个任意的扰动)
另外还有一种方法就是初始的时候将扰动写成αδu ,这时对于任意可取的δu ,其能量函数写成W (u +αδu ) 。回到微积分的基本概念,去寻找W 对于α的极值点:
如果我们将它看成是对于α的Taylor 展开,就可以找出其一阶导数(对于极值点必须为零),由于δu 是任意可取的,我们可以得到和前面相同的结果。
小结
上面的过程省略很多推导步骤,如果大家对推导有兴趣,可以试着自己推导。要说明一下的是:
1、 变量δu (而不是它的梯度)必须是很小而且是任意的。
2、 这里没有考虑边界条件和体力,比如重力等等。我们前面所讨论的问题局限于一个没有
任何约束和载荷的边界条件的区域上。
3、 一般来说δu 的限制比多元微积分中δx 宽松。在泛函中,只要δu 是在容许的范围内即可,
也就是u +δu 必须和物理位移场相对应。理解这个意思对理解有限元弱形式非常重要。
考虑边界条件和体力
如前面所讲,弹性能的泛函形式是不完整的,因为它没有加上相应的边界条件和载荷。弹性能的单位是N ⋅m ,也就是力乘上位移。在边界上,我们一般施加面力,或者指定位移,单位为m 。一般来说,我们希望附加形式是“面力乘上长度”。同样的方式可以对体力进行处理F 。
在数学上,结构场的边界条件分为两类。第一类直接定义边界上的力:
其中第一项由定义域内的方程所确定,第二项称为弹簧常数q ,等式右边是面力g 。这种边界条件就是我们通常说得流量或者Nuemann 边界条件。
第二类边界条件就是
定义一个固定的或者Dirichlet 边界条件。如果h 是矩阵的形式,r 就是定义了边界上的指定位移。固定边界条件不能直接加入泛函中去,但是可以通过反力间接加上去。当指定位移边界时,可以描述一个反力μ(N /mm ),也就是弹性体可以在固定处保持不变。反力就是我们这里用到的Lagrange 乘子,通过添加反力到力作用处的边界,可以忽略到固定边界类型。这时候我们可以形成统一的边界条件:
2
这里R 是原始的固定边界,μ是需要计算的反力。在前面的简化形式中,h 和r 都是常数,所以上式可以变化为:
记住,方程中的每一项都是矢量,表示各个方向的面力。为了得到所做的功(能量),必须点乘上位移u 。
通过合并一些系数项,将外力写成P ,可简化表达式,这时边界条件可以写成:
n ⋅(c ∇u ) =P
对于其他物理场,可能P 代表边界上的源项。 注意到上式和Navier 方程非常接近:
将能量泛函展开:
关键推导
这个时候,我们又要在u 上添加上δu ,可得:
零阶项就是泛函本身,第一阶项是:
这个方程是非常重要的一项。
从前面的讨论可知,我们应该重新组合多项式,保证带有δu 的被积函数成为一项。如果可以做到,因为δu 是任意的(事实上必须是在容许范围内),我们知道这一项必须为零。这是我们能找到极值点的唯一方法。右边第一项需要进一步处理得到我们需要的形式。 第一项我们可以根据Green 公式(有时候可能采用的是Stokes 原理)进行分部积分:
利用c 的对称性,我们可以得到:
利用Green 公式得到:
将体积项和边界项合并起来:
确定极值点,必须有:
上式应该对于任何δu 都成立。因此体积项必须有:
边界项上有:
现在我们又回到PDE 问题上了,这说明泛函的理论解就是PDE 方程的解,即通过能量最小化原理又重新推回到了PDE 形式上!这也是说明最小能量化和PDE 形式本质上是统一的一个数学证明。
GJ :要注意的是,实际有限元的求解不是从泛函又导回PDE 方程,而是通过网格划分离散化,得到数值近似解。而通过建立的泛函求解数值近似解应该算是比较完备的方法,这里说完备,是和弱形式对比来的,两者的区别后面会说。
弱形式
GJ :上面是利用能量最小化形式或者变分原理建立的泛函,即PDE 的等效积分形式,下面说的是弱形式建立。
那么,到底什么是弱形式呢?Navier 方程的弱形式实际上已经在前面的推导过程中出现过了,即一阶变分的原形式:
如果我们回到COMSOL Multiphysics的文档(或者是关于有限元和弱形式的书籍中),会发现所谓的试函数相当于扰动δu ,
v =δu
弹性静力学PDE 方程的弱形式
为了更好的理解弱形式,我们必须丢弃前面讨论的能量最小化原理,转向一种更加抽象的方法。弱形式之所以比能量最小化原理更强大,是因为它还可以应用到一些没有得到较好的能量定义的问题中。
首先我们考虑弹性静力学的PDE 方程
边界条件是:
抽象的过程如下:乘上容许范围内的试函数v ,在感兴趣的域内积分可得:
对左侧利用Green 公式进行分部积分:
应用PDE 方程的边界条件,可以得到:
整理可得:
这就是PDE 方程的弱形式。如果在积分区域内对于试函数v 都是有效的,则上式和PDE 方程是等效的。PDE 方程的解称为强解,而弱形式的解称为弱解。二者唯一的区别是弱形式对于积分参数的连续性要求比PDE 形式低。由于变形梯度和弹性张量在弱形式里面都不需要微分,所以对函数连续性要求没有那么严格,而在PDE 形式中,所有的变量都处在散度的算子下,这要求这些变量必须是可微的。在弱形式中对于可微的要求放松了(一阶)。 同时,注意到弱形式和前面的一阶变分形式保持了一致,弱形式也可以作为虚功原理的一种推广。只是虚功原理中的位移δu 换成了更加抽象的试函数v 。如果弱形式解和能量最小化原理不一致的时候,极值点变成了鞍点。也就是,在弱形式中,仍可以将试函数理解为一种推广了的虚位移。
一般性问题的弱形式
正如前面所提到的,弱形式只是PDE 方程的一种推广形式,它对变量的连续性要求比较低。那么能量方法呢?如果有一个定义好了的能量来最小化,那么能量法和弱形式是一致的。但是,在下列情形下,弱形式更具有适用性:假如PDE 方程没有相对应的能量可以进行最小化。在这种情况下,弱形式仍然是适用的。由于弱形式对解的要求较低,所以说弱形式比PDE 和能量最小化适用范围更广泛。
GJ :弱形式和最小能量形式的区别就在于虚位移δu 与试函数v 的差别,如下面两式
也就是说,泛函求极值即为泛函的变分为0,如上面的式子1,所以泛函的有限元解对任意扰动δu 成立,而从式子2可以看到,弱形式的解只是对自己设定的试函数v 成立。所以泛函求极值得到近似函数是弱形式的特殊形式,即弱形式的试函数v 可以任意取而求得的近似函数,所以从这种意义上说泛函形式求得的近似解更完备。
但很多情况下无法得到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
● ∇u 的分量标记为ux ,uy 和uz 。
● 试函数v 的标记为u_test。
● ∇v 的分量标记为ux_test,uy_test,uz_test
● 只需要输入被积函数,它将被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 中的形式非常重要,尤其是采用弱约束的时候。更多详情可以参考COMSOL Multiphysics文档以及和Lagrange 乘子相关的技术文档。
尽管弱形式是一个标量表达式,但是COMSOL Multiphysics中,弱形式有和PDE 系统一样多的未知量需要文本输入。原因在于不同的多物理场问题可能需要不同的有限元分析类型和保证其数值稳定型的积分算法。对于3D 结构分析,弱形式中有三个文本输入框。但是,在离散之前,采用了同样的有限元单元和积分类型进行合并,这样就可以选择不同的弱形式进行操作。例如你可以在第一个域内选择弱形式,而其他的域内设置为空白。
对于流动问题的Navier-stokes 方程,情况又稍微有些不同。和未知的速度场相比,未知的压力采用一个低阶有限元来离散。这种情况下,不能将所有的弱项全部在同一个弱域内输入。为了保证数值稳定型,必须依靠混阶有限元(mixed finite element) 。混阶有限元并不是COMSOL Multiphysics特别制定,而是数值算法所需要的。
有限元方法
本章说明弱形式如何利用有限元方法来进行离散。假设我们需要离散以下扩散问题:
这是一个对流-扩散问题的特殊情况,其中α=0,β=0。
有限元的基本实现是将整个计算域Ω离散为多个特别简单的形状的小单元,比如2D 中的三角形,3D 中的四面体等等。相应的网格,例如三角形,由边和节点组成。下一步就是要选择一个比较容易实现的一些近似方法,其中一种比较简单的方法就是将解表示为采用线性多项式插值的所谓基函数的和。基函数的构造方法是指定某个节点为1,而相邻的节点为0,二者之间的值就是从0到1线性变化。这里说的相邻指的是中间有一条边将其连接起来。
遍历三角形网格的所有节点(从1到N )。定义节点i 的基函数为ϕi ,也就是在节点i 处其值为1,其他点处值为0。注意ϕi 只是在节点i 及其相邻的三角形内不为零。现在假设真实值u 可以用基函数的求和u h 来近似描述:
参数U i 是u h 在节点i 的值。同样,我们可以对试函数进行类似处理:
下标h 表示离散函数属于由所有三角形边中最长边表示的具有确定的网格尺寸h 的网格。
由于我们可以任意选择试函数,因此可以将除了j 点以外的所有的V j 设置为零,
接下来我们将所有的试函数(j =1,...,N )输入到弱形式中去,每个试函数都可以得到一
个方程。这样可以生成一个线性代数系统,系统矩阵就是我们所说的刚度矩阵。为什么我们可以自由选择试函数,不妨回想一下前面提到的弱形式需要对所有可取的试函数成立。选择试函数是有限元方法的重要环节,因为他在很大程度上影响着刚度矩阵。由于刚度矩阵中很多项为零,所以一般是稀疏矩阵。
当我们使用试函数的时候,生成的有限元刚度矩阵应该是一个方阵。如果弱形式本身定义良好的话,刚度矩阵应该是非奇异的,也就是说系统有一个唯一解。
现在考虑扩散方程的弱形式:
将表达式写成离散形式:
方程重新排列:
采用矩阵标注可得:
或者:
在这里刚度矩阵K 是:
解矢量U 的单元为U i ,载荷矢量L 的单元为,
现在我们明白为什么选择基函数和试函数很关键了。如果我们关注刚度矩阵K ij ,会发现其中很多元素为零,因为前面已经提到每个ϕj 都是大部分为零,同样ϕj 的梯度也是大部分为零的。有很多有效的算法去求解这类稀疏矩阵,COMSOL Multiphysics 提供一套稀疏线性系统求解器。
有限元方法同样适用于非线性问题。非线性方法一般来说采用迭代的算法,每一次迭代就是求解一个与上面类似的线性弱形式方法。
抽象和几何解释
为什么有限元可以解决问题,它是如何解决问题的?前面的讨论中可以找到一些答案。为了有一个更清晰的答案,我们需要了解一些更多的泛函的概念。我们将发现有限元方法通过一种优化方法将解投影到一个有限维函数空间来求解。
标量积
为了得到有限元方法的几何解释,或脑海中的意象,我们需要熟悉标量积的概念。在线性代数中,我们知道两个矢量f 和g 的标量积为
这里的矢量f 和g 属于3维矢量空间。标量积可以推广到任意维N ,
标量积有时也称做内积。
如果两个矢量正交,
一个矢量f 可以通过标量积投影到另一个矢量
g
其中f p 是与g 平等的投影矢量:
矢量差
与g 正交:
投影矢量的唯一性特征是原始和投影矢量之间的差与它所投影的矢量正交。如果需要找到在方向g 上与f 最近的矢量,f p 就是我们的答案。矢量e 可以看作是关于f 到f p 之间的近似误差。换句话说,误差矢量e 与矢量g 正交。后面在对有限元进行几何解释时将用到这个结论。但首先我们得介绍一些泛函分析的概念。
与矢量不同,泛函分析讲的是函数,它们必须属于无限维的矢量空间。我们可以积分形式定义一个无限维矢量空间(函数空间)中的标量积:
如果两个函数正交,则有
进一步将标量积的概念推广,并且考虑包含函数梯度的被积函数,
下标1和2用来说明上面是两种不同的标量积。
Hilbert 空间
如果对于一个标量积,我们只考虑那些与自身进行的标量积(积分)有一个有限值的函数u ,
我们说这些函数属于一个确定的函数类,或函数空间。由标量积和定义域Ω组成的函数空间就被称为Hilbert 空间。
对应于上面的标量积1的Hilbert 空间通常标记为L 2,与标量积2对应的Hilbert 空间常称为H 1。通常空间L 2中的函数比H 1中的多,因为H 1中的函数自动地存在于L 2中,反之则不一定。我们也可以将这种现象称为H 1是L 2的一个子空间,或
有限元方法的抽象形式
现在让我们考虑PDE 问题[3]:
在域Ω中,边界条件为u =0。这是对流-扩散方程的一种特例,其中c =a =1,β=0。通过分部积分,表明当试函数选择成在边界上具有相同的边界条件v =0时,弱形式中的边界项为0。得到弱形式为:
其中包含前面提到的两种标量积,可改写为:
找到u ,使得对于所有属于相配的Hilbert 空间中的v ,有:
这里的相配的空间是H 1,且在边界上v =0。解函数u 也必须属于这个空间。注意,对于前面提到的弱形式,我们可以很自由地添加不同的标量积,因为每个积的结果是实数或虚数。 现在我们选择前面讨论的基函数的和(线性组合)来近似u 和v ,近似解被称作u h 和v h 。这些函数属于Hilbert 空间,可称为V h ,由线性基函数ϕi 扩展而来。V h 的空间维数为N (基函数的数目)。此外,V h 是H 1的子空间,
换句话说,如果函数属于V h ,则它自动地属于H 1。空间H 1是一个大得多的空间,因为它是无限维的。有限的基函数不可能扩散成属于H 1的所有函数。
弱形式的有限元现在变成了:
对于V h 中的所有v h ,在V h 中找到u h ,使得
现在我们对有限元方法在函数空间的作为一个确定的投影的几何解释有了更深入的了解。最后,在原始的弱形式中:
用v h 代替v 。这是合理的,因为v 在H 1中,而v h 在V h 中,由于V h 是H 1的子空间,因此v h 也是在H 1中,
现在从弱形式中减去有限元解,得到:
或
即离散误差:
与所有的V h 中的v h 正交。也就是说,有限元解u h 是真实解u 在属于H 1的有限维子空间V h
的投影。最终,我们得到了关于有限元方法的几何解释,
对于给定的网格,有限元解u h 是在函数空间V h 中关于标量积(,) 2最接近真实解u 的解。