波前法及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列的元素非零，其余元素为零。我们知道，总刚度矩阵的列数与行数是相等的，在本列情况下，列数等于24。在第1行当中，只有6个元素非零，其余18个元素都是零。同理，第4节点只与第1和第2单元相连，其影响也就只限于第1和第2单元。因而，在总刚度矩阵第4行当中，第1、2、5、4、8、7列的元素非零，其余18个元素为零。第2节点影响4个单元，分别是第1、9、8、16单元，因而在总刚度矩阵第2行当中，非零元素最多，达到9个，即第1、2、3、4、5、6、22、23、24列的元素非零，其余15个元素为零。如此，可想而知，总刚度矩阵的每一行

的大部分元素都是零。 现在我们要考虑怎样利用这些零元素了。在以上网格中，共有16个单元，24个节点。使用高斯消去法，会生成24 x 24 的总刚度矩阵，即有24行，24列。而使用波前法，总刚度矩阵的行数虽然不变，也是 24， 但列数仅为11（至于为什么是11 ，下面要详细讨论）。最重要的，在本例情况下，是波前法根本就不需要将 24 x 11 的总刚度矩阵存在内存中，而是存在硬盘上的。内存里，波前

法只需要放一个 11 x 11 的所谓波前矩阵就行。

那么，什么是波前矩阵呢？

就是在某一时刻，彼此发生关系的节点的刚度系数组成的矩阵。这个矩阵是方

的，其中含有最多非零元素的那一行就叫波前。 什么叫某一时刻？就是某一单元！如前面MATLAB程序所示，计算有限元刚度矩阵有两重循环，最外面那层循环，是对单元循环，即从编号为第一的单元，到编号最大的单元。在波前法中，当循环到某一单元时，在计算该单元刚度矩阵以后，还要看看哪一个节点可以消去了，也就是消元。被消元的节点，对以后其它单元刚度矩阵就不再有影响了，该节点的刚度系数就可以存入硬盘指定文件中，而这些系数就可以从波前矩阵中清除掉，以便空出位置来，存储其它节点信息。因此，一个节点可以被消元的时间，是可以用单元循环的进度来度量的。 那么，一个节点，什么时候可以消元了？就是与该节点相连的所有单元都循环到了的时候。例如，若循环从第1单元开始，当循环完了第2单元（计算单元矩阵并装配到波前矩阵中）时，第4节点就可以消元了，因为第4节点所连接的2个单元都循环到了。同理，当循环完了第3单元时，第7节点也可以消元了。而第1节点的消元要等到循环到第8单元。而第2、3节点的消元时间最迟，要

循环到最后一个单元，第16单元。 据此，可以编制一个节点消元时间表，也就是循环到了哪个单元，该节点便可以被消元了。算法很简单，就是查找每一个节点所连接的最大单元编号。程序如下：

程序二：计算节点消元时间

以上程序中，第三行的ICO变量是个两维数组，18行，4列。它的每一行代表一个单元，该行的4列给出该单元的四个节点。这段程序执行的结果，是在一维数组变量NodeDeactiveTime中定义每个节点可以消元的时间（即单元号）。此

时，NodeDeactiveTime 的值是：

8 16 16 2 10 10 3 11 11 4 12 12 5 13 13

6 14 14 7 15 15 8 16 16。 第1个数“8”代表第1节点的消元时间是8（单元）；第2个数“16”代表第2节点的消元时间是16（单元）；余类推。请注意，消元最早的节点是第4节点，时间是“2”（单元）。其次是第7节点，时间是“3”（单元）。我们下面介绍在波前

矩阵里如何消元时要用到这两个信息。

知道了各节点的消元时间，就可以计算波前矩阵的最大阶数了。程序如下：

程序三：计算波前矩阵的最大阶数

在以上程序中，第1行开了一维辅助数组，NodeInFront，用于记录每一个节点是不是在波前矩阵中，“1”表示在，“0”表示不在。第2行将我们要计算的波前矩阵的最大阶数maxFrontWidth的值初始为零。第3行开始对单元循环。对编号为 i 的单元，第4行从单元总表ICO中取出该单元的节点（本例为四个节点）；第5行将这些（四个）节点在NodeInFront中的值定为“1”，代表它们进入了波前矩阵。注意，此时，有的节点可能已经在波前矩阵中了，即它们在NodeInFront中的值已经是“1”。但这没有关系，现在只是重新再植一次“1”，再一次表示该节点在波前矩阵中。第6行计算 maxFrontWidth。就是将NodeInFront中的“1”相加，再与当前maxFrontWidth比较，谁的值更大，maxFrontWidth就取谁

的值。也就是说，maxFrontWidth的值只增不减。 第8行对 i 单元的（四个）节点循环，查找其中每个节点是不是到了消元时间（由NodeDeactiveTime给出）。 如果是的，第10行的逻辑变量deactive的值为“1”，并在第12行将该节点在NodeInFront中的值改为零。这表示该节

点被清除出波前矩阵。 这段程序结束全部循环后，便得到maxFrontWidth的值为11，就是波前矩阵的最大阶数。前面说的波前法总刚度矩阵是24 行 x 11列。其中的11，就是这样得出的。它的含义，就是在整个计算过程中，某一时刻，同时存在于波前矩阵

的节点数，其值最大为11。 以上程序，实际上模拟了节点进入和离开波前矩阵的装配和消元过程。下表给出波前矩阵中的节点随着单元循环的整个变化过程。第1列是单元号码。第2列是将该单元的刚度矩阵装入波前矩阵后，波前矩阵中的节点。第3列是这些节点的数目。第4列是此时刻可以消元的节点。第5列是消元后，将消元节点清除

消元过程

消元也是波前法里比较难于理解的技术细节。我在调试波前法程序过程中，所花时间最多的，就是消元。而现在，当我想解释消元过程时，发现遇到了更大的困难。因为我不知道应该从哪里说起。消元是个细活儿，必得有研究针尖上能站立多少天使的牛角尖精神才能弄懂。不是所有人都有这种精神。而为写此文，已经花了如此多工夫，几乎超过了开发程序的时间！

所以，还是算了吧。一简对三繁，我就把消去第一个可消节点（第4节点）的过

程列表于下。您若看得明白，就算彻底懂了波前法。 为方便说明问题，我将原问题的内壁边界条件换为绝热，外壁换为温度给定（程序四第63至64行须作相应修改）。这样，第4节点就不经过边界处理那一段程

序（第183至184行），而是进入消元那一段（第186至190行）。

表1，将第1单元的单元矩阵装配入波前矩阵后的波前矩阵

然后计算第2单元的单元刚度矩阵，结果如下：

注意，第2单元的单元刚度矩阵恰好与第1单元相等，这纯属巧合。巧合的原因是：1. 两单元节点关于45°线对称；2. 介质各向同性。

表3，将第2单元的单元矩阵装配入波前矩阵后的波前矩阵

表4，第 4 节点消元后的波前矩阵

红字显示因为消元而出现的非零元素。例如，第7、8两节点，本来与第1节点不共单元，对第1节点没有影响，但现在有了。在消元中，我们用第1行减去第4行（乘除两个系数后），因而第4行的元素全部出现在第1行中。

消元的结果，是第4行主元素所在的列，除主元素外，其余元素全被消成零。此时，就可以将第4行存入指定文件中，如程序第203行所示。程序第205行将此时的波前节点号写入另一指定文件。以后回代需要。程序第204行写方程Ax =

B的等号右端自由项。 有关数据存入文件后，就可以将第4行清零了（程序第206至208行）。此时

波前矩阵如下表所示。

表5，第 4 节点清除后的波前矩阵

波前法与高斯消去法的效率之比较

比较结果见下表。第1列是比较所用的6种网格。Cdiv 代表沿圆周的分割，Tdiv代表沿厚度方向的分割。第2列是节点数。第3列是波前矩阵的最大宽度。第4列和第5列分别是使用程序四（波前法）和程序一（高斯消去法）解题的时钟时

间（elapsed time）。

由上表可看出，当节点数很少（24和80）时，两种方法算题的时间差不多，但当节点数达到288时，差别开始显示出来，高斯消去法需时比波前法将近多一倍。64 x 16的网格（1088节点），时间差别达到15倍。96 x 24 （2400节电），差别到了60倍。到了最后一个网格，128 x 32（4224节点），运转程序一（高斯消去法）需内存563MB，这已经超出了我的计算机内存（512MB），因而所

需时间可认为无穷大。

小结：波前法的过去、现在和未来

在人类计算技术发展史上，有限元方法可以说是个奇迹。从来没有任何一种计算方法，能像有限元这样，深刻而广泛地影响着人们的日常生活，它使得上百的人成为大师，上千的人成为专家，上万的人靠它吃饭，上亿的人每天乘坐使用有限

元程序设计的飞机、汽车、轮船、火车。 在有限元半个世纪的发展史当中，无数人写了无数关于有限元的文献，可用汗牛充栋来形容，时至今日，已经无人能够全部读完（也没必要）。但1970年Irons发表的那篇关于波前法的文章，大概要算有限元文献里被引用次数最多的一篇了。波前法的提出，突破了制约有限元应用的瓶颈，使得当时内存量很小的计算机，可以解算规模很大的题目。可以说，没有波前法，就没有有限元的大发展，

也就没有有限元的今天。 但是，波前法的使用，也有上限。就是波前矩阵不能超过计算机内存。以我的计算机为例，516MB内存。假设有一半可用来计算，就是约256MB。这相当于256x1024x1024/8个双精度实数。将此数开平方，得到约5700。也就是说，在我的计算机上，能使用波前法解算的题目，其波前矩阵的最大阶数，约为5000。

也就是说，如果是解算三维问题，其变量总数，大约是五、六万。 所以，今天，经典波前法，例如本文程序四，您只能在大学里搞科研的那些教授自己开发的有限元恐龙化石里找得到了。今天，在像ANSYS那种商业有限元程

序里，经典波前法早已被波前法的变种所取代，就是多波前法。 简单地说，多波前法，就是将网格分成许多子区域，在每个子区域上使用波前法。这样可以将波前矩阵的阶数降低，以突破内存容许的解题上限并提高效率。每个

子区域间数据的依赖关系与传递，是用所谓消去树来管理的。 毋庸置疑，这样的多波前法，算法更加复杂，且编程量巨大，可能一万行都嫌少！

当然，商业有限元程序里，不光是多波前法，还有其它方法；不光是直接解法，还有迭代法。这些程序，光解方程的部分，据介绍，有的已经达到三万行。 但，万变不离其宗。欲理解有限元程序，欲开发有限元程序，必须理解波前法。

参考文献：

Pascal JOLY，《Mise en Oeuvre de la Méthode des Eléments Finis》。

J.N. REDDY, 《An Introduction to The Finite Element Method》, McGraw-Hill Book Company，1984

周洪伟，吴舒，陈璞，《有限元分析快速直接求解技术进展》，“力学进展”杂志，第37卷，第2期，175-188页，2007年5月25日。

本文来源：https://www.bwwdw.com/article/kqe6.html