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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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