波前法及matlab实现

有限元二维热传导波前法MATLAB程序•二维热传导有限元•使用高斯消去法解线性方程组的二维热传导有限元程序•波前法的基本概念与算法•使用波前法解线性方程组的二维热传导有限元程序•消元过程•波前法与高斯消去法的效率之比较•小结：波前法的过去、现在和未来波前法是求解线性方程组的一种方法，广泛用于有限元程序。

它最初由英国人（？）B.M. Irons于1970在“国际工程计算方法杂志”上发表。

30多年来，波前法有了不少变种。

本文所用算法，采于法国人Pascal JOLY所著《Mise en Oeuvre de la Méthode des Eléments Finis》。

这本书是我1993年在比利时一家书店买的，书中有一节"波前法"，六页纸，解释了基本概念和算法，但没有程序，也没有细节讨论。

我曾花了两个半天的时间，在网上寻找波前法程序，或更详细的资料，没有找到（需要花钱才能看的文献不算）。

倒是看到不少中国人，也在寻找。

一些人说，波前法程序太难懂了。

通过自己编写程序，我同意这些人的说法，确实难。

我还真很少编如此耗费脑力的程序。

完工之后，我曾对朋友老王说，有了算法，编程序还这么难，当初想出算法的人，真是了不起。

现将我对波前法的理解和编程体会解说如下，供感兴趣的网友参考，也为填补网络上波前法空白。

二维热传导有限元波前法和有限元密不可分。

因而，在编写波前法程序之前，必须有个有限元程序。

为了简化问题，最好是能解算一个节点上只有一个自由度的问题的有限元程序，而且要尽可能地简单。

我手边现有的有限元程序都不符合这个要求。

就决定先开发一个解算二维热传导问题的MATLAB有限元程序。

二维热传导问题的微分方程是其中T 是温度，Kx 和Ky 分别是x 和y 方向上的热传导系数，q 是热源。

对于这样的比较经典的二阶微分方程，如何导出有限元表达式？这个问题，是有限元的首要问题！我相信，许许多多学过有限元的人，甚至每天都在用有限元的人，并不真的十分明了。

我自己曾经就是这样。

在我于2005年3月到巴西教书之前，我搞过20年有限元，其中有六年，还是在比利时一个专门开发有限元程序的公司里度过的。

但是，如果您那时问我，任给一个二阶偏微分方程，例如上述热传导方程，如何导出有限元表达？说老实话，不看书，我还真就不会！直到2006年，我教了一遍有限元后，才弄明白（惟教书才是最好的学习）。

其实简单极了：只需将那二阶偏微分方程，乘上一个任意标量函数，然后在某个有限单元上积分。

请看下列推导：即其中，Ωe 是单元面积；Φ 是任意标量函数。

注意在以上积分中，温度要遭受二阶偏导，这很不好。

有限元的精髓，在于通过分步积分，将其中一阶偏导转移到那个任意标量函数Φ 上，这样就可降低一阶温度偏导，改善对它的苛刻待遇。

这得用到您在高等数学最后一章里学过的散度定理（Theorem of divergence）：其中，Г 是面积Ω的边界；（反）Δ 是梯度算子；F 和G 是任意两个矢量函数。

对于二维问题，上述散度定理可写为：将此散度定理应用于（2）式中的第一项积分，令得：将此积分结果带入（3）式，得到热传导偏微分方程的弱化表达（weak form）：所谓“弱化”，是指对温度函数的可导阶数要求降低了。

在原热传导偏微分方程（1）中，温度函数必须是二阶可偏导函数，而在弱化表达（6）中，温度函数只要一阶偏导就行了。

有限单元法就是以偏微分方程的弱化表达式为出发点的。

现在，到了说明那个任意标量函数，Φ，是何方神圣的时候了。

有限单元法有各种各样的变种，而最常见的，应用最广泛的，是所谓迦辽金法（Galerkin），即将这个任意标量函数，等同于，有限元的插值函数。

迦辽金法的优点是可以最终形成对称劲度矩阵，从而具有较好的数值稳定性。

我们知道，在一个有限单元中，任意一点的值，例如温度，是用节点上的温度表达的。

对于一个四节点双线性单元来说，设四个节点的温度分别是T1,T2,T3,T4，则单元内任意一点的温度T 可表达为其中Ψ1，Ψ2，Ψ3，Ψ4，为插值函数，也称为形函数，定义如下：其中ξ 和η 称为形坐标，取值区间为[-1,1]。

基于式（7），对温度的偏导数可写为：将此二偏导数代入弱化表达式（6），该式就转化为以节点温度为变量的代数方程：到此，为得到原偏微分方程的有限元表达，我们只需将任意标量函数，Φ，依次取为四个插值函数，Ψ1-4，就得到四个代数方程：注意到式(9)-(12)中下标的规律，我们可将这四个方程合并写成矩阵的形式：使用下标表达，式（13）可进一步缩写为：式中对应于下标i = 1…4 的每一个值，下标j 取遍1…4。

式（14）就是热传导偏微分方程（1）的四节点双线性有限元表达，也是我们接下来编写有限元程序的出发点。

使用高斯消去法解线性方程组的二维热传导有限元程序这个程序是专为解一个特殊的热传导问题而设计的。

所解问题是：已知一个无限长圆筒，内径100毫米，外径200毫米，筒内壁表面温度1°C，外壁绝热，求圆筒截面上的温度分布。

圆筒材料各向均质，热传导系数为 1 （单位还得查一下，但也无所谓）。

问题的解答很简单：均布，截面各点温度都是1°C。

为什么？因为外部绝热，没有热量损失。

温度只能是均布。

而内壁温度为1°C，所以到处都是1°C。

因为问题的几何图形简单，有限元网格便容易自动生成。

在以下程序中，第12行至第51行，生成四节点单元的有限元网格。

控制变量有两个，Cdiv 和Tdiv。

前者定义沿圆周分成多少单元；后者定义沿径向即筒壁厚度方向分成多少单元。

如果Cdiv = 8，Tdiv = 2，则所生成有限元网格如下（由第52行子程序DrawFEM画出）：第64行使用MATLAB命令syms 定义了两个符号变量ksi(即前面公式中的ξ),eta(η)。

在MATLAB中，可对符号变量进行代数运算，例如定义公式，求导，积分等。

第72行就利用代数运算定义了本文前面给出的四节点单元的形函数。

第76和77行分别对形函数关于ksi 和eta 求导。

第78至第99行，计算这些导数在高斯积分点的数值。

本程序中，每个单元有四个高斯积分点，也就是说，在ksi 和eta 两个方向上，各有二个高斯积分点。

第101至124行，根据式（14）计算单元劲度矩阵，Kelem，并将其装配入总劲度矩阵K。

本问题没有热源，所以在单元循环水平上，不需计算式（14）中的热源积分项。

第127至139行，施加边界条件。

圆筒外壁是绝热条件，即法向热流等于零。

在有限元中，这是自然满足的，所以式（14）中的边界热流积分项为零，不必考虑。

唯一需要考虑的，是圆筒内壁温度等于1°C 。

这种温度给定的边界条件，在数学上叫第一类边界条件。

在有限元技术中，有各种各样的方法施加第一类边界条件。

主要是考虑提高内存效率。

鉴于本程序的目的是进一步开发波前法，不需要仔细考虑如何更有效地施加这种温度给定的边界条件，因而所用的方法是最简单的一种：即将内壁边界节点的各行方程，全部换为T = Tin (Tin = 1°C)。

相应地，将这些行的主元素所占据的列，乘以Tin后，移到等号右边。

这样处理后，就得到待解的线性方程组：K x = F。

第141行，使用本人自编的高斯消去法函数，egauss，解出为未知量x，也就是个节点的温度，都等于1。

这一行，也可直接用MATLAB命令，x = K \ F, 取代。

我使用egauss的目的，是为下一步与波前法解方程比较效率。

程序一：使用高斯消去法解线性方程组的二维热传导有限元程序波前法的基本概念与算法有限元形成的线性方程组的系数矩阵，即刚度矩阵，是稀疏矩阵，也就是说，矩阵里含有许许多多的零元素。

有限元网格节点数目越巨大，非零元素与零元素的比值越小，刚度矩阵越稀疏。

用普通高斯消去法求解这样的线性方程组，完全不考虑矩阵的稀疏性，对大量的零元素也进行加减乘除，结果浪费了大量时间。

不仅如此，应用高斯消去法，因为需要将整个刚度矩阵存在计算机内存里，所需计算机内存量与有限元网格节点未知量总数成平方的关系。

以热传导问题为例。

一千个节点，光存刚度矩阵，就需内存1000 x 1000 x 8 / （1024 x 1024） = 7.8 Mb。

这还没问题。

但若要计算一万个节点的问题，就需要10 x 10 x 7.8 = 780 Mb 来存刚度矩阵。

内存为1 Gb的计算机已经不能计算这样的问题了，因为微软视窗等其它系统程序还需要内存。

您当然可以开始这样的计算。

微软视窗发现内存不够时，会自动启动虚拟内存，也就是把硬盘上的某一块地方，当作内存来使用。

您的计算仍然能进行。

但是，您一定看不见那最终的计算结果，除非几个月内不断电，计算机不出毛病。

原因在于，与内存相比，虚拟内存的存取时间可认为是无限长！在这种情况下，因为普通高斯消去法需要极其频繁地使用虚拟内存，它的解算时间也就无限地延长了。

但如果您在这样的计算机上运行ANSYS，或任何需要花钱买的有限元程序，就会发现，解算同样的问题，只需几分钟。

您甚至可以毫无困难地解算十万个节点的热传导问题。

秘密就在于，这些商业有限元软件，在求解线性方程组时，都尽可能地利用刚度矩阵的稀疏性。

波前法就是这样一种充分考虑了刚度矩阵的稀疏性求解方程的方法。

刚度矩阵为什么会稀疏？因为在有限元中，一个节点的影响，只限于它所连接的那些单元。

为什么？就是因为在前面，我们推导有限元时，在式（2）中，将热传导偏微分方程乘以的那个神奇函数Φ。

我们说过，Φ是任意标量函数。

既然是任意的，当然可以任意选取。

然而我们没有“任意”地、胡乱地选取，而是居心叵测地，让它恰恰等于有限元的插值函数！而这些插值函数，恰恰只在给定单元内部非零，在单元外边一律为零！换句话说，插值函数只是些局部函数。

我们让Φ等于插值函数，它也就具有了这种局部性。

正是Φ的这种局部性，使得一个节点的影响，只限于它所连接的单元。

有限元方法，之所以能够在计算力学领域里令人眼花缭乱的各种各样的计算方法中，独领风骚，傲视群雄，鹤立鸡群，至今几达50年之久，其全部奥妙，皆出于此！为进一步考察这些影响到底有多少，我们来看下面的例子。

使用前面的MATLAB 有限元程序，给Cdiv 的值输入8, Tdiv 输入 2 ，就会生成以下网格。

它将圆周分成8份，厚度分成2份。

图中括弧内是单元号码，其余数字为节点号码。

可以看出，第1节点只与第1和第8单元相连，其影响也就只限于这两个单元。

这里所说的影响，就是在刚度矩阵中，第1和第8单元的所有节点，即第1、2、5、4、22、23节点，要发生关系。

也就是说，在总刚度矩阵（高斯消去法中的K矩阵）中，相应的行与列上的元素非零。

例如在第1行当中，第1、2、5、4、22、23列的元素非零，其余元素为零。