波前法及matlab实现

./有限元二维热传导波前法MATLAB程序二维热传导有限元使用高斯消去法解线性方程组的二维热传导有限元程序波前法的基本概念与算法使用波前法解线性方程组的二维热传导有限元程序消元过程波前法与高斯消去法的效率之比较小结：波前法的过去、现在和未来波前法是求解线性方程组的一种方法,广泛用于有限元程序。它最初由英国人〔？B.M.Irons于1970在"国际工程计算方法杂志"上发表。30多年来,波前法有了不少变种。本文所用算法,采于法国人PascalJOLY所著《MiseenOeuvredelaMéthodedesElémentsFinis》。这本书是我1993年在比利时一家书店买的,书中有一节"波前法",六页纸,解释了基本概念和算法,但没有程序,也没有细节讨论。我曾花了两个半天的时间,在网上寻找波前法程序,或更详细的资料,没有找到〔需要花钱才能看的文献不算。倒是看到不少中国人,也在寻找。一些人说,波前法程序太难懂了。

通过自己编写程序,我同意这些人的说法,确实难。我还真很少编如此耗费脑力的程序。完工之后,我曾对朋友老王说,有了算法,编程序还这么难,当初想出算法的人,真是了不起。

现将我对波前法的理解和编程体会解说如下,供感兴趣的网友参考,也为填补网络上波前法空白。二维热传导有限元波前法和有限元密不可分。因而,在编写波前法程序之前,必须有个有限元程序。为了简化问题,最好是能解算一个节点上只有一个自由度的问题的有限元程序,而且要尽可能地简单。我手边现有的有限元程序都不符合这个要求。就决定先开发一个解算二维热传导问题的MATLAB有限元程序。

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

对于这样的比较经典的二阶微分方程,如何导出有限元表达式？

这个问题,是有限元的首要问题！我相信,许许多多学过有限元的人,甚至每天都在用有限元的人,并不真的十分明了。

我自己曾经就是这样。在我于20XX3月到巴西教书之前,我搞过20年有限元,其中有六年,还是在比利时一个专门开发有限元程序的公司里度过的。但是,如果您那时问我,任给一个二阶偏微分方程,例如上述热传导方程,如何导出有限元表达？说老实话,不看书,我还真就不会！

直到20XX,我教了一遍有限元后,才弄明白〔惟教书才是最好的学习。

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

请看下列推导：即其中,Ωe是单元面积；Φ是任意标量函数。注意在以上积分中,温度要遭受二阶偏导,这很不好。有限元的精髓,在于通过分步积分,将其中一阶偏导转移到那个任意标量函数Φ上,这样就可降低一阶温度偏导,改善对它的苛刻待遇。这得用到您在高等数学最后一章里学过的散度定理〔Theoremofdivergence：其中,Г是面积Ω的边界；〔反Δ是梯度算子；F和G是任意两个矢量函数。

对于二维问题,上述散度定理可写为：将此散度定理应用于〔2式中的第一项积分,令得：将此积分结果带入〔3式,得到热传导偏微分方程的弱化表达〔weakform：所谓"弱化",是指对温度函数的可导阶数要求降低了。在原热传导偏微分方程〔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后,移到等号右边。这样处理后,就得到待解的线性方程组：Kx=F。

第141行,使用本人自编的高斯消去法函数,egauss,解出为未知量x,也就是个节点的温度,都等于1。这一行,也可直接用MATLAB命令,x=K\F,取代。我使用egauss的目的,是为下一步与波前法解方程比较效率。程序一：使用高斯消去法解线性方程组的二维热传导有限元程序波前法的基本概念与算法有限元形成的线性方程组的系数矩阵,即刚度矩阵,是稀疏矩阵,也就是说,矩阵里含有许许多多的零元素。有限元网格节点数目越巨大,非零元素与零元素的比值越小,刚度矩阵越稀疏。用普通高斯消去法求解这样的线性方程组,完全不考虑矩阵的稀疏性,对大量的零元素也进行加减乘除,结果浪费了大量时间。不仅如此,应用高斯消去法,因为需要将整个刚度矩阵存在计算机存里,所需计算机存量与有限元网格节点未知量总数成平方的关系。以热传导问题为例。一千个节点,光存刚度矩阵,就需存1000x1000x8/〔1024x1024

=7.8Mb。这还没问题。但若要计算一万个节点的问题,就需要10x10x7.8=780Mb来存刚度矩阵。存为1Gb的计算机已经不能计算这样的问题了,因为微软视窗等其它系统程序还需要存。您当然可以开始这样的计算。微软视窗发现存不够时,会自动启动虚拟存,也就是把硬盘上的某一块地方,当作存来使用。您的计算仍然能进行。但是,您一定看不见那最终的计算结果,除非几个月不断电,计算机不出毛病。原因在于,与存相比,虚拟存的存取时间可认为是无限长！在这种情况下,因为普通高斯消去法需要极其频繁地使用虚拟存,它的解算时间也就无限地延长了。

但如果您在这样的计算机上运行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个节点。使用高斯消去法,会生成24x24的总刚度矩阵,即有24行,24列。而使用波前法,总刚度矩阵的行数虽然不变,也是24,但列数仅为11〔至于为什么是11,下面要详细讨论。最重要的,在本例情况下,是波前法根本就不需要将24x11的总刚度矩阵存在存中,而是存在硬盘上的。存里,波前法只需要放一个11x11的所谓波前矩阵就行。那么,什么是波前矩阵呢？就是在某一时刻,彼此发生关系的节点的刚度系数组成的矩阵。这个矩阵是方的,其中含有最多非零元素的那一行就叫波前。什么叫某一时刻？就是某一单元！如前面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行x11列。其中的11,就是这样得出的。它的含义,就是在整个计算过程中,某一时刻,同时存在于波前矩阵的节点数,其值最大为11。以上程序,实际上模拟了节点进入和离开波前矩阵的装配和消元过程。下表给出波前矩阵中的节点随着单元循环的整个变化过程。第1列是单元。第2列是将该单元的刚度矩阵装入波前矩阵后,波前矩阵中的节点。第3列是这些节点的数目。第4列是此时刻可以消元的节点。第5列是消元后,将消元节点清除之后,波前矩阵里剩下的节点数,即第3列减去第4列。可以看出,两次达到最大波前,分别当循环到第7单元和第10单元时。读到这里,对照前面那个有限元网格,您要是能够理解这个表中所有数字的来龙去脉,您就理解了波前法！剩下的,是一些技术上的细节,可以通过阅读下列全部程序来最后彻底"抠"懂。这些技术细节里,比较难懂的,是如何将单元刚度矩阵装配入波前矩阵。单元刚度矩阵是个4阶矩阵,即4x4矩阵,代表了单元4个节点之间的相互影响。例如,第1行里的4列元素,是单元4个节点对单元第1节点的影响；第2行里的4列元素,是那4个节点对单元第2节点的影响；余类推。

每个单元节点的局部编号都是1、2、3、4,但它们的总节点编号是不同的,而且是唯一的。单元节点的局部编号与其总节点编号的对应关系,在本文的程序中,由二维数组ICO给出。在用普通高斯消去法解方程时,任给一个单元刚度矩阵,我们可以通过ICO表,查到该单元4个节点的总编号。而节点总编号与总刚度矩阵K的行与列之间具有一一对应的关系。第1节点占据第1行第1列,第2节点占据第2行第2列,……,第n节点占据第n行第n列。所以,使用普通高斯消去法解方程时,把单元刚度矩阵装配进总刚度矩阵的算法很简单,如本文第一个程序的第122行所示,只一行程序便可实现。打个比方,这种情况下的总刚度矩阵K,就像占地广阔的大观园,有许多房子,《红楼梦》里的所有人等,例如金陵十二钗,各有各的住所,还有许多空地〔相当于K里的零元素,可供黛玉葬花,姐妹们赏月吟诗。波前法的装配算法则要复杂得多。因为波前矩阵小,不能同时装下所有节点的所有元素。这种情况,就如同旅店。世界上所有旅店是住不下装得下世界上所有居民的,旅客必须有进有出才行。节点进出波前矩阵,就如同旅客进出旅馆。有的旅客住的时间长,有的短。节点在波前矩阵呆的时间也一样,有长有短。随着节点的进进出出,波前矩阵的每一行和列都可能先后被不同节点占用,就像旅馆里的每个房间都会被不同旅客先后住过。所以,波前矩阵的行与列,与总节点编号之间,不再有像使用普通高斯消去法解方程时那样的一一对应的固定关系。单元矩阵的某一行、某一列的元素,应该放到波前矩阵的哪一行、哪一列,是动态分配的。有两种可能：如果某节点已经存在于波前矩阵中,那么就把该节点在单元矩阵中的元素加到波前矩阵的相应元素上〔因而需要知道它在哪里；反之,如果某节点还没进入过波前矩阵,那么就在波前矩阵给它分配一个自由的行与列。在下面的程序中,我们使用一维数组freeLines来记录波前矩阵每一行〔同时也是每一列。我们假定行数等于列数,也就是说,一个节点占据了第n行,它也就占据了第n列的状态：0表示自由；大于

0表示已占用〔即占用该行的节点>。我们还使用一维数组GlobalID2FrontId来记录每一节点在波前矩阵占据的行数〔=列数。也就是说,给出一个节点的总编号,就能在GlobalID2FrontId中查到它在波前矩阵中的位置。这个位置如果是零,就表示该节点还没进入过波前矩阵。这些算法的具体实施可见下列程序中的第148-168行。程序四：使用波前法解线性方程组的二维热传导有限元程序消元过程消元也是波前法里比较难于理解的技术细节。我在调试波前法程序过程中,所花时间最多的,就是消元。而现在,当我想解释消元过程时,发现遇到了更大的困难。因为我不知道应该从哪里说起。消元是个细活儿,必得有研究针尖上能站立多少天使的牛角尖精神才能弄懂。不是所有人都有这种精神。而为写此文,已经花了如此多工夫,几乎超过了开发程序的时间！

所以,还是算了吧。一简对三繁,我就把消去第一个可消节点〔第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列分别是使用程序四〔波前法和程序一〔高斯消去法解题的时钟时间〔elapsedtime。由上表可看出,当节点数很少〔24和80时,两种方法算题的时间差不多,但当节点数达到288时,差别开始显示出来,高斯消去法需时比波前法将近多一倍。64x16的网格〔1088节点,时间差别达到15倍。96x24〔2400节电,差别到了60倍。到了最后一个网格,128x32〔4224节点,运转程序一〔高斯消去法需存563MB,这已经超出了我的计算机存〔512MB,因而所需时间可认为无穷大。小结：波前法的过去、现在和未来在人类计算技术发展史上,有限元方法可以说是个奇迹。从来没有任何一种计算方法,