波前法及matlab實(shí)現(xiàn)

./有限元二維熱傳導(dǎo)波前法MATLAB程序二維熱傳導(dǎo)有限元使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序波前法的基本概念與算法使用波前法解線性方程組的二維熱傳導(dǎo)有限元程序消元過程波前法與高斯消去法的效率之比較小結(jié)：波前法的過去、現(xiàn)在和未來波前法是求解線性方程組的一種方法,廣泛用于有限元程序。它最初由英國人〔？B.M.Irons于1970在"國際工程計(jì)算方法雜志"上發(fā)表。30多年來,波前法有了不少變種。本文所用算法,采于法國人PascalJOLY所著《MiseenOeuvredelaMéthodedesElémentsFinis》。這本書是我1993年在比利時一家書店買的,書中有一節(jié)"波前法",六頁紙,解釋了基本概念和算法,但沒有程序,也沒有細(xì)節(jié)討論。我曾花了兩個半天的時間,在網(wǎng)上尋找波前法程序,或更詳細(xì)的資料,沒有找到〔需要花錢才能看的文獻(xiàn)不算。倒是看到不少中國人,也在尋找。一些人說,波前法程序太難懂了。

通過自己編寫程序,我同意這些人的說法,確實(shí)難。我還真很少編如此耗費(fèi)腦力的程序。完工之后,我曾對朋友老王說,有了算法,編程序還這么難,當(dāng)初想出算法的人,真是了不起。

現(xiàn)將我對波前法的理解和編程體會解說如下,供感興趣的網(wǎng)友參考,也為填補(bǔ)網(wǎng)絡(luò)上波前法空白。二維熱傳導(dǎo)有限元波前法和有限元密不可分。因而,在編寫波前法程序之前,必須有個有限元程序。為了簡化問題,最好是能解算一個節(jié)點(diǎn)上只有一個自由度的問題的有限元程序,而且要盡可能地簡單。我手邊現(xiàn)有的有限元程序都不符合這個要求。就決定先開發(fā)一個解算二維熱傳導(dǎo)問題的MATLAB有限元程序。

二維熱傳導(dǎo)問題的微分方程是其中T是溫度,Kx和Ky分別是x和y方向上的熱傳導(dǎo)系數(shù),q是熱源。

對于這樣的比較經(jīng)典的二階微分方程,如何導(dǎo)出有限元表達(dá)式？

這個問題,是有限元的首要問題！我相信,許許多多學(xué)過有限元的人,甚至每天都在用有限元的人,并不真的十分明了。

我自己曾經(jīng)就是這樣。在我于20XX3月到巴西教書之前,我搞過20年有限元,其中有六年,還是在比利時一個專門開發(fā)有限元程序的公司里度過的。但是,如果您那時問我,任給一個二階偏微分方程,例如上述熱傳導(dǎo)方程,如何導(dǎo)出有限元表達(dá)？說老實(shí)話,不看書,我還真就不會！

直到20XX,我教了一遍有限元后,才弄明白〔惟教書才是最好的學(xué)習(xí)。

其實(shí)簡單極了：只需將那二階偏微分方程,乘上一個任意標(biāo)量函數(shù),然后在某個有限單元上積分。

請看下列推導(dǎo)：即其中,Ωe是單元面積；Φ是任意標(biāo)量函數(shù)。注意在以上積分中,溫度要遭受二階偏導(dǎo),這很不好。有限元的精髓,在于通過分步積分,將其中一階偏導(dǎo)轉(zhuǎn)移到那個任意標(biāo)量函數(shù)Φ上,這樣就可降低一階溫度偏導(dǎo),改善對它的苛刻待遇。這得用到您在高等數(shù)學(xué)最后一章里學(xué)過的散度定理〔Theoremofdivergence：其中,Г是面積Ω的邊界；〔反Δ是梯度算子；F和G是任意兩個矢量函數(shù)。

對于二維問題,上述散度定理可寫為：將此散度定理應(yīng)用于〔2式中的第一項(xiàng)積分,令得：將此積分結(jié)果帶入〔3式,得到熱傳導(dǎo)偏微分方程的弱化表達(dá)〔weakform：所謂"弱化",是指對溫度函數(shù)的可導(dǎo)階數(shù)要求降低了。在原熱傳導(dǎo)偏微分方程〔1中,溫度函數(shù)必須是二階可偏導(dǎo)函數(shù),而在弱化表達(dá)〔6中,溫度函數(shù)只要一階偏導(dǎo)就行了。

有限單元法就是以偏微分方程的弱化表達(dá)式為出發(fā)點(diǎn)的。

現(xiàn)在,到了說明那個任意標(biāo)量函數(shù),Φ,是何方神圣的時候了。

有限單元法有各種各樣的變種,而最常見的,應(yīng)用最廣泛的,是所謂迦遼金法〔Galerkin,即將這個任意標(biāo)量函數(shù),等同于,有限元的插值函數(shù)。迦遼金法的優(yōu)點(diǎn)是可以最終形成對稱勁度矩陣,從而具有較好的數(shù)值穩(wěn)定性。

我們知道,在一個有限單元中,任意一點(diǎn)的值,例如溫度,是用節(jié)點(diǎn)上的溫度表達(dá)的。對于一個四節(jié)點(diǎn)雙線性單元來說,設(shè)四個節(jié)點(diǎn)的溫度分別是T1,T2,T3,T4,則單元任意一點(diǎn)的溫度T可表達(dá)為其中Ψ1,Ψ2,Ψ3,Ψ4,為插值函數(shù),也稱為形函數(shù),定義如下：其中ξ和η稱為形坐標(biāo),取值區(qū)間為[-1,1]。

基于式〔7,對溫度的偏導(dǎo)數(shù)可寫為：將此二偏導(dǎo)數(shù)代入弱化表達(dá)式〔6,該式就轉(zhuǎn)化為以節(jié)點(diǎn)溫度為變量的代數(shù)方程：到此,為得到原偏微分方程的有限元表達(dá),我們只需將任意標(biāo)量函數(shù),Φ,依次取為四個插值函數(shù),Ψ1-4,就得到四個代數(shù)方程：注意到式<9>-<12>中下標(biāo)的規(guī)律,我們可將這四個方程合并寫成矩陣的形式：使用下標(biāo)表達(dá),式〔13可進(jìn)一步縮寫為：式中對應(yīng)于下標(biāo)i=1…4的每一個值,下標(biāo)j取遍1…4。

式〔14就是熱傳導(dǎo)偏微分方程〔1的四節(jié)點(diǎn)雙線性有限元表達(dá),也是我們接下來編寫有限元程序的出發(fā)點(diǎn)。使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序這個程序是專為解一個特殊的熱傳導(dǎo)問題而設(shè)計(jì)的。所解問題是：已知一個無限長圓筒,徑100毫米,外徑200毫米,筒壁表面溫度1°C,外壁絕熱,求圓筒截面上的溫度分布。圓筒材料各向均質(zhì),熱傳導(dǎo)系數(shù)為1〔單位還得查一下,但也無所謂。問題的解答很簡單：均布,截面各點(diǎn)溫度都是1°C。為什么？因?yàn)橥獠拷^熱,沒有熱量損失。溫度只能是均布。而壁溫度為1°C,所以到處都是1°C。

因?yàn)閱栴}的幾何圖形簡單,有限元網(wǎng)格便容易自動生成。在以下程序中,第12行至第51行,生成四節(jié)點(diǎn)單元的有限元網(wǎng)格?？刂谱兞坑袃蓚€,Cdiv和Tdiv。前者定義沿圓周分成多少單元；后者定義沿徑向即筒壁厚度方向分成多少單元。如果Cdiv=8,Tdiv=2,則所生成有限元網(wǎng)格如下〔由第52行子程序DrawFEM畫出：第64行使用MATLAB命令syms定義了兩個符號變量ksi<即前面公式中的ξ>,eta<η>。在MATLAB中,可對符號變量進(jìn)行代數(shù)運(yùn)算,例如定義公式,求導(dǎo),積分等。第72行就利用代數(shù)運(yùn)算定義了本文前面給出的四節(jié)點(diǎn)單元的形函數(shù)。第76和77行分別對形函數(shù)關(guān)于ksi和eta求導(dǎo)。第78至第99行,計(jì)算這些導(dǎo)數(shù)在高斯積分點(diǎn)的數(shù)值。本程序中,每個單元有四個高斯積分點(diǎn),也就是說,在ksi和eta兩個方向上,各有二個高斯積分點(diǎn)。

第101至124行,根據(jù)式〔14計(jì)算單元勁度矩陣,Kelem,并將其裝配入總勁度矩陣K。本問題沒有熱源,所以在單元循環(huán)水平上,不需計(jì)算式〔14中的熱源積分項(xiàng)。

第127至139行,施加邊界條件。圓筒外壁是絕熱條件,即法向熱流等于零。在有限元中,這是自然滿足的,所以式〔14中的邊界熱流積分項(xiàng)為零,不必考慮。唯一需要考慮的,是圓筒壁溫度等于1°C

。這種溫度給定的邊界條件,在數(shù)學(xué)上叫第一類邊界條件。在有限元技術(shù)中,有各種各樣的方法施加第一類邊界條件。主要是考慮提高存效率。鑒于本程序的目的是進(jìn)一步開發(fā)波前法,不需要仔細(xì)考慮如何更有效地施加這種溫度給定的邊界條件,因而所用的方法是最簡單的一種：即將壁邊界節(jié)點(diǎn)的各行方程,全部換為T=Tin<Tin=1°C>。相應(yīng)地,將這些行的主元素所占據(jù)的列,乘以Tin后,移到等號右邊。這樣處理后,就得到待解的線性方程組：Kx=F。

第141行,使用本人自編的高斯消去法函數(shù),egauss,解出為未知量x,也就是個節(jié)點(diǎn)的溫度,都等于1。這一行,也可直接用MATLAB命令,x=K\F,取代。我使用egauss的目的,是為下一步與波前法解方程比較效率。程序一：使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序波前法的基本概念與算法有限元形成的線性方程組的系數(shù)矩陣,即剛度矩陣,是稀疏矩陣,也就是說,矩陣?yán)锖性S許多多的零元素。有限元網(wǎng)格節(jié)點(diǎn)數(shù)目越巨大,非零元素與零元素的比值越小,剛度矩陣越稀疏。用普通高斯消去法求解這樣的線性方程組,完全不考慮矩陣的稀疏性,對大量的零元素也進(jìn)行加減乘除,結(jié)果浪費(fèi)了大量時間。不僅如此,應(yīng)用高斯消去法,因?yàn)樾枰獙⒄麄€剛度矩陣存在計(jì)算機(jī)存里,所需計(jì)算機(jī)存量與有限元網(wǎng)格節(jié)點(diǎn)未知量總數(shù)成平方的關(guān)系。以熱傳導(dǎo)問題為例。一千個節(jié)點(diǎn),光存剛度矩陣,就需存1000x1000x8/〔1024x1024

=7.8Mb。這還沒問題。但若要計(jì)算一萬個節(jié)點(diǎn)的問題,就需要10x10x7.8=780Mb來存剛度矩陣。存為1Gb的計(jì)算機(jī)已經(jīng)不能計(jì)算這樣的問題了,因?yàn)槲④浺暣暗绕渌到y(tǒng)程序還需要存。您當(dāng)然可以開始這樣的計(jì)算。微軟視窗發(fā)現(xiàn)存不夠時,會自動啟動虛擬存,也就是把硬盤上的某一塊地方,當(dāng)作存來使用。您的計(jì)算仍然能進(jìn)行。但是,您一定看不見那最終的計(jì)算結(jié)果,除非幾個月不斷電,計(jì)算機(jī)不出毛病。原因在于,與存相比,虛擬存的存取時間可認(rèn)為是無限長！在這種情況下,因?yàn)槠胀ǜ咚瓜シㄐ枰獦O其頻繁地使用虛擬存,它的解算時間也就無限地延長了。

但如果您在這樣的計(jì)算機(jī)上運(yùn)行ANSYS,或任何需要花錢買的有限元程序,就會發(fā)現(xiàn),解算同樣的問題,只需幾分鐘。您甚至可以毫無困難地解算十萬個節(jié)點(diǎn)的熱傳導(dǎo)問題。

秘密就在于,這些商業(yè)有限元軟件,在求解線性方程組時,都盡可能地利用剛度矩陣的稀疏性。波前法就是這樣一種充分考慮了剛度矩陣的稀疏性求解方程的方法。

剛度矩陣為什么會稀疏？因?yàn)樵谟邢拊?一個節(jié)點(diǎn)的影響,只限于它所連接的那些單元。為什么？就是因?yàn)樵谇懊?我們推導(dǎo)有限元時,在式〔2中,將熱傳導(dǎo)偏微分方程乘以的那個神奇函數(shù)Φ。我們說過,Φ是任意標(biāo)量函數(shù)。既然是任意的,當(dāng)然可以任意選取。然而我們沒有"任意"地、胡亂地選取,而是居心叵測地,讓它恰恰等于有限元的插值函數(shù)！而這些插值函數(shù),恰恰只在給定單元部非零,在單元外邊一律為零！換句話說,插值函數(shù)只是些局部函數(shù)。我們讓Φ等于插值函數(shù),它也就具有了這種局部性。正是Φ的這種局部性,使得一個節(jié)點(diǎn)的影響,只限于它所連接的單元。有限元方法,之所以能夠在計(jì)算力學(xué)領(lǐng)域里令人眼花繚亂的各種各樣的計(jì)算方法中,獨(dú)領(lǐng)風(fēng)騷,傲視群雄,鶴立雞群,至今幾達(dá)50年之久,其全部奧妙,皆出于此！

為進(jìn)一步考察這些影響到底有多少,我們來看下面的例子。使用前面的MATLAB有限元程序,給Cdiv的值輸入8,Tdiv輸入2,就會生成以下網(wǎng)格。它將圓周分成8份,厚度分成2份。圖中括弧是單元,其余數(shù)字為節(jié)點(diǎn)?？梢钥闯?第1節(jié)點(diǎn)只與第1和第8單元相連,其影響也就只限于這兩個單元。這里所說的影響,就是在剛度矩陣中,第1和第8單元的所有節(jié)點(diǎn),即第1、2、5、4、22、23節(jié)點(diǎn),要發(fā)生關(guān)系。也就是說,在總剛度矩陣〔高斯消去法中的K矩陣中,相應(yīng)的行與列上的元素非零。例如在第1行當(dāng)中,第1、2、5、4、22、23列的元素非零,其余元素為零。我們知道,總剛度矩陣的列數(shù)與行數(shù)是相等的,在本列情況下,列數(shù)等于24。在第1行當(dāng)中,只有6個元素非零,其余18個元素都是零。同理,第4節(jié)點(diǎn)只與第1和第2單元相連,其影響也就只限于第1和第2單元。因而,在總剛度矩陣第4行當(dāng)中,第1、2、5、4、8、7列的元素非零,其余18個元素為零。第2節(jié)點(diǎn)影響4個單元,分別是第1、9、8、16單元,因而在總剛度矩陣第2行當(dāng)中,非零元素最多,達(dá)到9個,即第1、2、3、4、5、6、22、23、24列的元素非零,其余15個元素為零。如此,可想而知,總剛度矩陣的每一行的大部分元素都是零。

現(xiàn)在我們要考慮怎樣利用這些零元素了。在以上網(wǎng)格中,共有16個單元,24個節(jié)點(diǎn)。使用高斯消去法,會生成24x24的總剛度矩陣,即有24行,24列。而使用波前法,總剛度矩陣的行數(shù)雖然不變,也是24,但列數(shù)僅為11〔至于為什么是11,下面要詳細(xì)討論。最重要的,在本例情況下,是波前法根本就不需要將24x11的總剛度矩陣存在存中,而是存在硬盤上的。存里,波前法只需要放一個11x11的所謂波前矩陣就行。那么,什么是波前矩陣呢？就是在某一時刻,彼此發(fā)生關(guān)系的節(jié)點(diǎn)的剛度系數(shù)組成的矩陣。這個矩陣是方的,其中含有最多非零元素的那一行就叫波前。什么叫某一時刻？就是某一單元！如前面MATLAB程序所示,計(jì)算有限元剛度矩陣有兩重循環(huán),最外面那層循環(huán),是對單元循環(huán),即從編號為第一的單元,到編號最大的單元。在波前法中,當(dāng)循環(huán)到某一單元時,在計(jì)算該單元剛度矩陣以后,還要看看哪一個節(jié)點(diǎn)可以消去了,也就是消元。被消元的節(jié)點(diǎn),對以后其它單元剛度矩陣就不再有影響了,該節(jié)點(diǎn)的剛度系數(shù)就可以存入硬盤指定文件中,而這些系數(shù)就可以從波前矩陣中清除掉,以便空出位置來,存儲其它節(jié)點(diǎn)信息。因此,一個節(jié)點(diǎn)可以被消元的時間,是可以用單元循環(huán)的進(jìn)度來度量的。

那么,一個節(jié)點(diǎn),什么時候可以消元了？就是與該節(jié)點(diǎn)相連的所有單元都循環(huán)到了的時候。例如,若循環(huán)從第1單元開始,當(dāng)循環(huán)完了第2單元〔計(jì)算單元矩陣并裝配到波前矩陣中時,第4節(jié)點(diǎn)就可以消元了,因?yàn)榈?節(jié)點(diǎn)所連接的2個單元都循環(huán)到了。同理,當(dāng)循環(huán)完了第3單元時,第7節(jié)點(diǎn)也可以消元了。而第1節(jié)點(diǎn)的消元要等到循環(huán)到第8單元。而第2、3節(jié)點(diǎn)的消元時間最遲,要循環(huán)到最后一個單元,第16單元。

據(jù)此,可以編制一個節(jié)點(diǎn)消元時間表,也就是循環(huán)到了哪個單元,該節(jié)點(diǎn)便可以被消元了。算法很簡單,就是查找每一個節(jié)點(diǎn)所連接的最大單元編號。程序如下：程序二：計(jì)算節(jié)點(diǎn)消元時間以上程序中,第三行的ICO變量是個兩維數(shù)組,18行,4列。它的每一行代表一個單元,該行的4列給出該單元的四個節(jié)點(diǎn)。這段程序執(zhí)行的結(jié)果,是在一維數(shù)組變量NodeDeactiveTime中定義每個節(jié)點(diǎn)可以消元的時間〔即單元號。此時,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個數(shù)"8”代表第1節(jié)點(diǎn)的消元時間是8〔單元；第2個數(shù)"16”代表第2節(jié)點(diǎn)的消元時間是16〔單元；余類推。請注意,消元最早的節(jié)點(diǎn)是第4節(jié)點(diǎn),時間是"2”〔單元。其次是第7節(jié)點(diǎn),時間是"3”〔單元。我們下面介紹在波前矩陣?yán)锶绾蜗獣r要用到這兩個信息。

知道了各節(jié)點(diǎn)的消元時間,就可以計(jì)算波前矩陣的最大階數(shù)了。程序如下：程序三：計(jì)算波前矩陣的最大階數(shù)在以上程序中,第1行開了一維輔助數(shù)組,NodeInFront,用于記錄每一個節(jié)點(diǎn)是不是在波前矩陣中,"1”表示在,"0”表示不在。第2行將我們要計(jì)算的波前矩陣的最大階數(shù)maxFrontWidth的值初始為零。第3行開始對單元循環(huán)。對編號為i的單元,第4行從單元總表ICO中取出該單元的節(jié)點(diǎn)〔本例為四個節(jié)點(diǎn)；第5行將這些〔四個節(jié)點(diǎn)在NodeInFront中的值定為"1”,代表它們進(jìn)入了波前矩陣。注意,此時,有的節(jié)點(diǎn)可能已經(jīng)在波前矩陣中了,即它們在NodeInFront中的值已經(jīng)是"1”。但這沒有關(guān)系,現(xiàn)在只是重新再植一次"1”,再一次表示該節(jié)點(diǎn)在波前矩陣中。第6行計(jì)算maxFrontWidth。就是將NodeInFront中的"1”相加,再與當(dāng)前maxFrontWidth比較,誰的值更大,maxFrontWidth就取誰的值。也就是說,maxFrontWidth的值只增不減。第8行對i單元的〔四個節(jié)點(diǎn)循環(huán),查找其中每個節(jié)點(diǎn)是不是到了消元時間〔由NodeDeactiveTime給出。如果是的,第10行的邏輯變量deactive的值為"1”,并在第12行將該節(jié)點(diǎn)在NodeInFront中的值改為零。這表示該節(jié)點(diǎn)被清除出波前矩陣。

這段程序結(jié)束全部循環(huán)后,便得到maxFrontWidth的值為11,就是波前矩陣的最大階數(shù)。前面說的波前法總剛度矩陣是24行x11列。其中的11,就是這樣得出的。它的含義,就是在整個計(jì)算過程中,某一時刻,同時存在于波前矩陣的節(jié)點(diǎn)數(shù),其值最大為11。以上程序,實(shí)際上模擬了節(jié)點(diǎn)進(jìn)入和離開波前矩陣的裝配和消元過程。下表給出波前矩陣中的節(jié)點(diǎn)隨著單元循環(huán)的整個變化過程。第1列是單元。第2列是將該單元的剛度矩陣裝入波前矩陣后,波前矩陣中的節(jié)點(diǎn)。第3列是這些節(jié)點(diǎn)的數(shù)目。第4列是此時刻可以消元的節(jié)點(diǎn)。第5列是消元后,將消元節(jié)點(diǎn)清除之后,波前矩陣?yán)锸Ｏ碌墓?jié)點(diǎn)數(shù),即第3列減去第4列?？梢钥闯?兩次達(dá)到最大波前,分別當(dāng)循環(huán)到第7單元和第10單元時。讀到這里,對照前面那個有限元網(wǎng)格,您要是能夠理解這個表中所有數(shù)字的來龍去脈,您就理解了波前法！剩下的,是一些技術(shù)上的細(xì)節(jié),可以通過閱讀下列全部程序來最后徹底"摳"懂。這些技術(shù)細(xì)節(jié)里,比較難懂的,是如何將單元剛度矩陣裝配入波前矩陣。單元剛度矩陣是個4階矩陣,即4x4矩陣,代表了單元4個節(jié)點(diǎn)之間的相互影響。例如,第1行里的4列元素,是單元4個節(jié)點(diǎn)對單元第1節(jié)點(diǎn)的影響；第2行里的4列元素,是那4個節(jié)點(diǎn)對單元第2節(jié)點(diǎn)的影響；余類推。

每個單元節(jié)點(diǎn)的局部編號都是1、2、3、4,但它們的總節(jié)點(diǎn)編號是不同的,而且是唯一的。單元節(jié)點(diǎn)的局部編號與其總節(jié)點(diǎn)編號的對應(yīng)關(guān)系,在本文的程序中,由二維數(shù)組ICO給出。在用普通高斯消去法解方程時,任給一個單元剛度矩陣,我們可以通過ICO表,查到該單元4個節(jié)點(diǎn)的總編號。而節(jié)點(diǎn)總編號與總剛度矩陣K的行與列之間具有一一對應(yīng)的關(guān)系。第1節(jié)點(diǎn)占據(jù)第1行第1列,第2節(jié)點(diǎn)占據(jù)第2行第2列,……,第n節(jié)點(diǎn)占據(jù)第n行第n列。所以,使用普通高斯消去法解方程時,把單元剛度矩陣裝配進(jìn)總剛度矩陣的算法很簡單,如本文第一個程序的第122行所示,只一行程序便可實(shí)現(xiàn)。打個比方,這種情況下的總剛度矩陣K,就像占地廣闊的大觀園,有許多房子,《紅樓夢》里的所有人等,例如金陵十二釵,各有各的住所,還有許多空地〔相當(dāng)于K里的零元素,可供黛玉葬花,姐妹們賞月吟詩。波前法的裝配算法則要復(fù)雜得多。因?yàn)椴ㄇ熬仃囆?不能同時裝下所有節(jié)點(diǎn)的所有元素。這種情況,就如同旅店。世界上所有旅店是住不下裝得下世界上所有居民的,旅客必須有進(jìn)有出才行。節(jié)點(diǎn)進(jìn)出波前矩陣,就如同旅客進(jìn)出旅館。有的旅客住的時間長,有的短。節(jié)點(diǎn)在波前矩陣呆的時間也一樣,有長有短。隨著節(jié)點(diǎn)的進(jìn)進(jìn)出出,波前矩陣的每一行和列都可能先后被不同節(jié)點(diǎn)占用,就像旅館里的每個房間都會被不同旅客先后住過。所以,波前矩陣的行與列,與總節(jié)點(diǎn)編號之間,不再有像使用普通高斯消去法解方程時那樣的一一對應(yīng)的固定關(guān)系。單元矩陣的某一行、某一列的元素,應(yīng)該放到波前矩陣的哪一行、哪一列,是動態(tài)分配的。有兩種可能：如果某節(jié)點(diǎn)已經(jīng)存在于波前矩陣中,那么就把該節(jié)點(diǎn)在單元矩陣中的元素加到波前矩陣的相應(yīng)元素上〔因而需要知道它在哪里；反之,如果某節(jié)點(diǎn)還沒進(jìn)入過波前矩陣,那么就在波前矩陣給它分配一個自由的行與列。在下面的程序中,我們使用一維數(shù)組freeLines來記錄波前矩陣每一行〔同時也是每一列。我們假定行數(shù)等于列數(shù),也就是說,一個節(jié)點(diǎn)占據(jù)了第n行,它也就占據(jù)了第n列的狀態(tài)：0表示自由；大于

0表示已占用〔即占用該行的節(jié)點(diǎn)>。我們還使用一維數(shù)組GlobalID2FrontId來記錄每一節(jié)點(diǎn)在波前矩陣占據(jù)的行數(shù)〔=列數(shù)。也就是說,給出一個節(jié)點(diǎn)的總編號,就能在GlobalID2FrontId中查到它在波前矩陣中的位置。這個位置如果是零,就表示該節(jié)點(diǎn)還沒進(jìn)入過波前矩陣。這些算法的具體實(shí)施可見下列程序中的第148-168行。程序四：使用波前法解線性方程組的二維熱傳導(dǎo)有限元程序消元過程消元也是波前法里比較難于理解的技術(shù)細(xì)節(jié)。我在調(diào)試波前法程序過程中,所花時間最多的,就是消元。而現(xiàn)在,當(dāng)我想解釋消元過程時,發(fā)現(xiàn)遇到了更大的困難。因?yàn)槲也恢缿?yīng)該從哪里說起。消元是個細(xì)活兒,必得有研究針尖上能站立多少天使的牛角尖精神才能弄懂。不是所有人都有這種精神。而為寫此文,已經(jīng)花了如此多工夫,幾乎超過了開發(fā)程序的時間！

所以,還是算了吧。一簡對三繁,我就把消去第一個可消節(jié)點(diǎn)〔第4節(jié)點(diǎn)的過程列表于下。您若看得明白,就算徹底懂了波前法。

為方便說明問題,我將原問題的壁邊界條件換為絕熱,外壁換為溫度給定〔程序四第63至64行須作相應(yīng)修改。這樣,第4節(jié)點(diǎn)就不經(jīng)過邊界處理那一段程序〔第183至184行,而是進(jìn)入消元那一段〔第186至190行。表1,將第1單元的單元矩陣裝配入波前矩陣后的波前矩陣然后計(jì)算第2單元的單元剛度矩陣,結(jié)果如下：注意,第2單元的單元剛度矩陣恰好與第1單元相等,這純屬巧合。巧合的原因是：1.兩單元節(jié)點(diǎn)關(guān)于45°線對稱；2.介質(zhì)各向同性。表3,將第2單元的單元矩陣裝配入波前矩陣后的波前矩陣

表4,第4節(jié)點(diǎn)消元后的波前矩陣紅字顯示因?yàn)橄霈F(xiàn)的非零元素。例如,第7、8兩節(jié)點(diǎn),本來與第1節(jié)點(diǎn)不共單元,對第1節(jié)點(diǎn)沒有影響,但現(xiàn)在有了。在消元中,我們用第1行減去第4行〔乘除兩個系數(shù)后,因而第4行的元素全部出現(xiàn)在第1行中。

消元的結(jié)果,是第4行主元素所在的列,除主元素外,其余元素全被消成零。此時,就可以將第4行存入指定文件中,如程序第203行所示。程序第205行將此時的波前節(jié)點(diǎn)號寫入另一指定文件。以后回代需要。程序第204行寫方程Ax=B的等號右端自由項(xiàng)。

有關(guān)數(shù)據(jù)存入文件后,就可以將第4行清零了〔程序第206至208行。此時波前矩陣如下表所示。表5,第4節(jié)點(diǎn)清除后的波前矩陣波前法與高斯消去法的效率之比較比較結(jié)果見下表。第1列是比較所用的6種網(wǎng)格。Cdiv代表沿圓周的分割,Tdiv代表沿厚度方向的分割。第2列是節(jié)點(diǎn)數(shù)。第3列是波前矩陣的最大寬度。第4列和第5列分別是使用程序四〔波前法和程序一〔高斯消去法解題的時鐘時間〔elapsedtime。由上表可看出,當(dāng)節(jié)點(diǎn)數(shù)很少〔24和80時,兩種方法算題的時間差不多,但當(dāng)節(jié)點(diǎn)數(shù)達(dá)到288時,差別開始顯示出來,高斯消去法需時比波前法將近多一倍。64x16的網(wǎng)格〔1088節(jié)點(diǎn),時間差別達(dá)到15倍。96x24〔2400節(jié)電,差別到了60倍。到了最后一個網(wǎng)格,128x32〔4224節(jié)點(diǎn),運(yùn)轉(zhuǎn)程序一〔高斯消去法需存563MB,這已經(jīng)超出了我的計(jì)算機(jī)存〔512MB,因而所需時間可認(rèn)為無窮大。小結(jié)：波前法的過去、現(xiàn)在和未來在人類計(jì)算技術(shù)發(fā)展史上,有限元方法可以說是個奇跡。從來沒有任何一種計(jì)算方法,