4.3.2　有限元法的原理及应用

4.3.2　有限元法的原理及应用

有限元法在早期是以变分原理为基础发展起来的，所以它广泛应用于以拉普拉斯方程和泊松方程所描述的各类物理场中。现代有限单元法的第一次成功尝试在1956年，Turner、Clough等人在分析飞机结构时，将刚架位移法推广应用于弹性力学问题，给出了用三角形单元求平面应力问题的正确答案。1960年，Clough进一步处理了平面弹性问题，并第一次提出了有限单元法，使人们认识到了它的功能。自1969年以来，有些学者在流体力学中应用加权余量法或最小二乘法等同样获得了有限元方程，这就使得有限元法可用于求解以任何微分方程所描述的物理场。

变分法和加权余量法是有限元法的重要基础。

变分法(Variational Methods)是17世纪末发展起来的一门数学分支，是处理函数的数学领域。它最终寻求的是极值函数，即变分法可使得泛函取得极大值或极小值。这里的泛函，通俗地理解就是函数的函数，也就是将y=f(x)中的自变量x用另一个函数来代替。变分法起源于一些具体的物理学问题，但最终由数学家们提出了问题的解决方案。

加权余量法(Weighted Residual Approach)，又称加权残量法、加权残余法。它是一种定义近似解(近似场函数)与真解(精确场函数)之间误差，并设法消除这种误差(或使其最小)的方法。在具体求解时，通常事先不知道真解满足什么规律，但确定真解分布可以用分段拟合的方法来逼近。这样就可以采用若干函数的线性组合来描述方程的解，也就是式(4.23)所描述的形式，其中近似解定义为

将代入要求解的微分方程和边界条件(初值条件)，与真解为φ满足的条件相比较就会出现误差，即残数或残值。若能设法构造消除残数的方程式，就可以把求解微分方程的问题转化为数值计算问题，从而得出近似解。式中，ψi称为基函数，基函数的选择要满足三个条件：①取自完备函数集；②线性独立；③要满足强制边界条件。n称为问题的自由度。一旦选定ψi，的结构就已知了，由此，式(4.23)中的待定系数Ci就成为确定的目标。

注意，在构造时要满足以下几个原则：

(1)要使得余数减小，我们可以通过寻求适当的待定系数来实现。

(2)当余数小于要求的精度时，就可以认为近似解就是偏微分方程的解。

(3)为有效表达减小余数的效果，还需选取适当的加权函数，以使余数和该加权函数的乘积积分为0。

下面，举例几个例子来讲述有限元法的原理。

例4.3　两极电容板空间分布如图4.7所示。

图4.7　平行板电容器及条件设置(书后附彩插)

已知控制方程和边界条件为

求平行板电容间的空间电势分布。

解：平行板电容器之间的空间电势分布具有均一性，因此可将本题中的三维问题简化为一维问题求解。设近似解为，定义场域Ω内的误差用RΩ表示，在边界Γ上的误差用RΓ表示，则其加权余量法误差(即余数)的定义为

设加 权 函 数：。那 么，目 标 函 数 就 变 成：∫ΩωjRΩdΩ+RΓdΓ，j=1，2，…，j表示试函数的线性组合项数。

伽辽金法是选择加权函数通常采用的方法，它选取基函数为加权函数，即

由此构建加权余量法的目标函数：

令Fj(R)=0，即余数最小时，趋于φ。通过此过程，已经将偏微分方程转化为j个代数方程组。

【解题过程】

(1)选取基函数ψi=xi(i=1，2)，构造近似解可得：

可见，中包含二次函数，这比只用一次函数的拟合效果更好。

(2)结合问题，对求解域写出余数表达式：

因此可得余数：。

(3)结合问题，对边界写出余数表达式：。那么在x=0处，，对应真解下的边界条件为。在x=d处，，对应真解下的边界条件为。

综上，可得在边界上的余量表达式为：和。

(4)写出总的加权余量的表达式：

当j=1时，得到一个代数方程：

当j=2时，又得到一个代数方程：

将式(4.31)、式(4.32)整理，得

(5)求解式(4.33)，得到待定系数，从而确定近似解。

解得：，因此可得近似解为

解毕!

事实上，求出的=φ，这是因为真解在本质上也是一个一次函数；而φ在C2=0时，自动退化为一个一次函数。这样就使得和φ具备了相同的结构特征，可使得两者误差为0。

图4.8　有源静电场的设置示意图

例4.3主要介绍了如何构造近似解，以及如何根据基函数的自由度来建立线性方程组的方法。这只是一个最基本的步骤，其中还没有涉及有限元的网格划分问题。接下来的例4.4将在例4.3的基础上进一步增加网格划分内容。

例4.4　平行平板电极如图4.8所示，求有源静电场的电势分布。其中，控制方程与边界条件为

题设中，由于电势分布关于x=0对称，因此求解x∈[0，1]区域内的电势分布即可；此外，根据平行平板电容内部电势分布的空间对称性特点，可将上述问题转化为一维问题。同时注意到，本例中的电势分布将受源电势的影响，因此电势的空间梯度将不再是线性的，这时就需要采用“分块逼近”的方式考虑，需要进行网格划分。本例是一维问题，因此网格用线段表示，将场域x∈[0，1]分割成e(1)、e(2)、e(3)、e(4)四个网格单元，如图4.9所示。从图中注意到，单元e(3)的长度要长一些，说明网格尺度不相同，这在有限元方法中是允许的。事实上，绝大多数问题的网格尺度都不是均一的，求解域内的网格有疏有密。

图4.9　一维空间的有限元网格划分示意图

当求解域被离散成有限个单元后，需要对每个单元进行处理。这种处理一方面包括求解单元节点上的值，另一方面还包括根据这些值对本单元节点以外的位置进行插值。例4.3中介绍了如何采用基函数构造近似解，通常情况下，基函数阶数越高，在一个单元内就可以描述更复杂的场函数变化规律，单元适应能力就越好。这种情况下，相对于低阶次基函数，则可以采用更少的单元划分来进行求解，这意味着方程组包含的方程会减少，但建立刚度矩阵的运算会更复杂。因此，对于一个特定的问题，都有一个相对适合的基函数阶次，它能够使总的计算时间最经济，而这一般需要根据计算经验决定。如果基函数为1阶，即采用线性插值的方法来表示分布势函数，则称为一阶有限元法；如果基函数高于1阶，即采用高阶插值法表示分布势函数，则称为高阶有限元法。这里的插值函数以基函数为基础，称为单元形函数，单元形函数是有限元法中一个十分重要的概念，这将在后面试函数的构造中具体介绍。

本例采用一阶有限元方法，先求5个节点的电位。相应基函数选择为ψi，ψi为一次函数，其近似解可表示为

加权余量表达式为

由于在边界处，电势为已知，所以必然有

这样第一类边界条件的余量就自动消除了。此时Fj(R)只剩下两项，即

将式(4.39)展开，得

对于式(4.40)的第1项，由格林第一定律可得

如此，则降低了微分阶数，进而有

根据在边界上=0；若设置加权函数wj=，式(4.42)将被大大简化。

可得

理论上，这里的权函数可以为常数或任意收敛的函数。因此，只要满足收敛条件，权函数就可以根据简化计算的要求任意选取，常用的方法包括配点法、子域法、最小二乘法、力矩法和伽辽金法。这里介绍一下经常采用的伽辽金法，伽辽金法(Galerkin method)是由苏联数学家鲍里斯·格里戈里耶维奇·伽辽金发明的一种数值分析方法，其原理为通过选取基函数作为加权余量法的权函数，以得到一组易于求解的线性代数方程。其中，边界条件能够自动满足。

应用伽辽金法，用基函数代替权函数可得wj=ψj，j为近似解的线性组合数。结合式(4.32)，可得

由于是线性微分算子，故微分、求和、积分次序可调换，得

进而对一个单元采用系数代替积分结果，有

相应的代数方程为

式(4.47)可简写为KC=f+b。只要得到kij、fj和bj的具体数值，就可以求得cj，进而求得。下面，结合特定的网格单元，重点探讨式(4.46)中的kij、fj和bj的求法。

一阶有限元法的试函数构造示意如图4.10所示，选择单元e(i)作为研究对象。

图4.10　一阶有限元法的试函数构造示意图

图中的Ni称为形函数(shape function)，它是一种连续函数，满足边界点的给定值和内部连续的条件。所选取的形函数可使得各个相邻单元节点处的解具有连续性，并且任意单元内的试函数值只与该单元节点坐标有关。试函数用于构造单元内真解的近似形。为了构造图4.10中的试函数形式，需要构造如图4.11所示的形函数。

图4.11　形函数示意图

e(i)单元的形函数包含第i个节点的形函数Ni和第i+1个节点的形函数Ni+1。根据式(4.36)中，对照试函数表达式，若令Ni=ψi，即采用基函数代替形函数，那么，近似解表达式中的系数Ci也就是要求解的近似解。对于构造的形函数，不论问题的空间维度如何，需具备三个性质：其一，形函数在节点i处的值为1，在其他节点处为0；其二，在单元e(i)中任意一点的形函数之和为1；其三，形函数的值变化区间是[0，1]。根据以上要求，选取的节点形函数相应数学形式为

式中，x——单元中的位置变量；

le(i)——单元长度。

有了以上构造的节点形函数，即可开始相应系数的计算。对于每个节点，相应的系数矩阵为

式中，i——节点序号；

j——试函数的线性组合项数，j=1，2。

本例中，第二类边界条件为0，因此在式(4.46)中，h=0，故式(4.49)中的b项为0，即得到

由式(4.48)可得

代入系数矩阵，可得

注意：每个单元的积分区间都是从第i节点开始，在第i+1节点结束，所以积分下限是0，积分上限是单元长度。

系数矩阵计算结果：

因为节点和节点之间相邻，互有影响，所以进行线性方程组的叠加，即矩阵封装：

得到的线性方程组为

注意：最后一个方程位于边界条件，该方程不完备，因此不能参与方程组的求解。

第5个节点上的值事实上已经由边界条件给出：。将其代入前4个方程，可求得其余4个待求节点的值，最终求得节点电势为：。

例4.4的精确解为，对照作图如图4.12所示。

图4.12　解析解与近似解对照图(书后附彩插)

思考：求x=0.68处的电势值。

首先，确定0.68位于第3个单元，涉及的节点及数值为：；然后，利用试函数公式可 得：

例4.3、例4.4主要介绍了一维情况下的有限元分析，其中例4.4相对完整地介绍了网格划分和用加权余量法求解节点方程和构造试函数等环节。通过对例4.4的学习，基本上就可以把握有限元方法的概貌和原理。下面是一个二维的例子，重要分析如何用有限元法来求解Navier-Stokes方程组。相对于前面的两个例子，例4.5增加了问题的维度和偏微分方程的复杂度，但求解思路是一致的。

例4.5　求解二维方腔流。控制方程为Navier-Stokes方程：

方腔及边界条件设置如下：

如图4.13所示，顶部边界为驱动边界，其他三个边界为静壁面。

本例参考和采用文献[4]中的求解方法，采用8节点矩形单元来设置网格。参考图4.14所示的变量设置，使用每个单元中的8个节点来描述速度分量u和v，使用4个角节点来描述压力p。

图4.13　顶盖流及边界设置示意图

图4.14　节点序号及变量设置

假设局部坐标系中的任意位置坐标为(ξ，η)，8个节点的相应坐标按图4.14中的标号顺序为(-1，-1)、(0，-1)、(1，-1)、(1，0)、(1，1)、(0，1)、(-1，1)、(-1，0)。角上4节点双线性矩形单元在局部坐标下的形函数的一般形式为

相应地，8节点二次矩形单元在局部坐标下的形函数形式为

利用形函数的Kronecker-delta性质，可以得到各节点形函数值的计算形式如下：

这样，速度分量采用二次插值函数，压力分量采用双线性插值函数，结果使得每个单元的速度和压力的未知变量是20个。因变量u、v、p表示为

式中，ui，vi——节点上的速度值；

pi——节点上的压力值。

用这些形函数来表示式(4.60)，可得

对式(4.79)采用Galerkin加权残差法，可得

由高斯定理，有

式中，Г——元素Ω的边界；

n——单元向外的单位法向量，n=(nx，ny)；

——Nj在边界Г法向方向上的方向导数，。

因此，有

对于狄利克雷边界条件，可得

对式(4.80)同理可得

利用权函数M l(·)对式(4.78)采用Galerkin加权残差法，得到

或

由于非线性的关系，所得到的代数方程组不能一次求解，需要迭代求解。在这种迭代解中，非线性项可以用许多不同的方法线性化。一种简单的做法是采用Picard线性化方法，其中非线性项被下式取代：

式中，，——速度分量的近似值。

假设单元的起始值为u01，u02，…，u08和v01，v02，…，v08，

迭代过程通过替换u0k和v0k来使计算继续，k=1，2，…，8，通过使用前两次迭代的速度分量值的平均值，直到满足允差为止。由式(4.84)、式(4.85)、式(4.87)可得矩阵形式的方程组为

其中，

这里，

然后，用局部坐标求导数和积分。设N(ξ，η)为局部坐标下的形函数，如果用x和y表示全局坐标，那么

式中，

为全局坐标系到局部坐标系的变换的雅可比矩阵，则有

接下来，获得雅可比矩阵的计算形式。设N1(ξ，η)，N2(ξ，η)，…，Nn(ξ，η)为某单元在局部坐标中的形函数。如果(x1，y1)，(x2，y2)，…，(xn，yn)为该单元节点的全局坐标，(x，y)为该单元上任一点的全局坐标，则

而且，

因此，变换的雅可比矩阵J为

离散有限元系统的形成需要计算单元上的积分。除了最简单的单元几何外，该积分不能解析求值。因此，数值积分是唯一的选择。常用的是高斯正交法。例如，使用坐标变换的微积分，一个典型的二维矩形单元的积分可以采用如下方法计算：

式中，ξi，ηj——高斯正交横坐标；

wi，wj——对应的权值；

在建立各单元的控制方程后，可进行单元方程的装配，建立适用于整个域的全局方程组。除了单元方程外，节点的全局坐标、单元连通性和节点的全局自由度也需用于建立全局系统方程。应用边界条件，则可得到修正的全局方程组。

方腔流是计算流体力学中的典型验证案例之一。本例中，设置雷诺数Re为1和10，网格包含100个正方形单元(803个节点)，对于两个速度分量采用10-6允差作为计算截止条件，相应的速度矢量场如图4.15所示。

图4.15　速度矢量场

(a)Re=1；(b)Re=10

与文献[4]中的结果对比发现，在Re=1和10时，速度矢量场的分布是基本吻合的。此外，文献[4]后附有相应的MATLAB程序，感兴趣的读者可以进一步学习和研究。

总体上，例4.5展示了一个采用二维有限元法求解流动的案例。由于其中涉及复杂的坐标变换操作，因此需要有一定数学知识基础才能更好地理解有限元法求解流动的机理。