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

1、有限元二維熱傳導(dǎo)波前法 MATLAB程序 二維熱傳導(dǎo)有限元 使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序 波前法的基本概念與算法 使用波前法解線性方程組的二維熱傳導(dǎo)有限元程序 消元過程 波前法與高斯消去法的效率之比較 小結(jié)：波前法的過去、現(xiàn)在和未來 波前法是求解線性方程組的一種方法，廣泛用于有限元程序。它最初由英國(guó)人(? ) B.M. Irons于1970在國(guó)際工程計(jì)算方法雜志”上發(fā)表。30多年來，波前 法有了不少變種。本文所用算法，采于法國(guó)人Pascal JOLY所著Mise en Oeuvre de la M thode des El ments Finis 。這本書是我1993年在比

2、利時(shí)一家書店買 的，書中有一節(jié)波前法，六頁紙，解釋了基本概念和算法，但沒有程序，也沒 有細(xì)節(jié)討論。我曾花了兩個(gè)半天的時(shí)間，在網(wǎng)上尋找波前法程序，或更詳細(xì)的資 料，沒有找到(需要花錢才能看的文獻(xiàn)不算)。倒是看到不少中國(guó)人，也在尋找。一些人說，波前法程序太難懂了。通過自己編寫程序，我同意這些人的說法，確實(shí)難。我還真很少編如此耗費(fèi)腦力 程程序。完工之后，我曾對(duì)朋友老王說，有了算法，編程序還這么難，當(dāng)初想出 算法的人，真是了不起?，F(xiàn)將我對(duì)波前法的理解和編程體會(huì)解說如下，供感興趣的網(wǎng)友參考，也為填補(bǔ)網(wǎng) 絡(luò)上波前法空白。二維熱傳導(dǎo)有限元波前法和有限元密不可分。因而，在編寫波前法程序之前，必須有個(gè)有限元程序

3、。 為了簡(jiǎn)化問題，最好是能解算一個(gè)節(jié)點(diǎn)上只有一個(gè)自由度的問題的有限元程序， 而且要盡可能地簡(jiǎn)單。我手邊現(xiàn)有的有限元程序都不符合這個(gè)要求。就決定先開發(fā)一個(gè)解算二維熱傳導(dǎo)問題的 MATLAB有限元程序。二維熱傳導(dǎo)問題的微分方程是其中T是溫度，Kx和Ky分別是x和y方向上的熱傳導(dǎo)系數(shù)，q是熱源 對(duì)于這樣的比較經(jīng)典的二階微分方程，如何導(dǎo)出有限元表達(dá)式？這個(gè)問題，是有限元的首要問題！我相信，許許多多學(xué)過有限元的人，甚至每天都在用有限元的人，并不真的十分 明了。我自己曾經(jīng)就是這樣。在我于 2005年3月到巴西教書之前，我搞過20年有限 元，其中有六年，還是在比利時(shí)一個(gè)專門開發(fā)有限元程序的公司里度過的 。但

4、是, 如果您那時(shí)問我，任給一個(gè)二階偏微分方程，例如上述熱傳導(dǎo)方程，如何導(dǎo)出有限元表達(dá)？說老實(shí)話，不看書，我還真就不會(huì)！直到2006年，我教了一遍有限元后，才弄明白（惟教書才是最好的學(xué)習(xí)）。其實(shí)簡(jiǎn)單極了：只需將那二階偏微分方程，乘上一個(gè)任意標(biāo)量函數(shù)，然后在某個(gè) 有限單元上積分。請(qǐng)看下列推導(dǎo)：即其中，Qe是單元面積；是任意標(biāo)量函數(shù)。注意在以上積分中，溫度要遭受二階偏導(dǎo)，這很不好。有限元的精髓，在于通過分步積分，將其中一階偏導(dǎo)轉(zhuǎn)移到那個(gè)任意標(biāo)量函數(shù)上，這樣就可降低一階溫度偏導(dǎo)，改善對(duì)它的苛刻待遇。這得用到您在高等數(shù)學(xué)最后一章里學(xué)過的散度定理（Theorem of divergence ）:j(5 f

5、) curf(VF)GdQ = -(VG)Tdft +J其中，r是面積。的邊界；（反）A是梯度算子；F和G是任意兩個(gè)矢量函 數(shù)。對(duì)于二維問題，上述散度定理可寫為：r dF clF. dG dGL Kh士次徹 J JI dx ffy ） J 1將此散度定理應(yīng)用于（2）式中的第一項(xiàng)積分，令Gx =GV =0日.將此積分結(jié)果帶入（3）式，得到熱傳導(dǎo)偏微分方程的弱化表達(dá)（weak form ）:所謂弱化；是指對(duì)溫度函數(shù)的可導(dǎo)階數(shù)要求降低了 。在原熱傳導(dǎo)偏微分方程（1） 中，溫度函數(shù)必須是二階可偏導(dǎo)函數(shù)，而在弱化表達(dá)（6）中，溫度函數(shù)只要一階偏導(dǎo)就行了。有限單元法就是以偏微分方程的弱化表達(dá)式為出發(fā)點(diǎn)的。

6、現(xiàn)在，到了說明那個(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)雙線性單元來說，設(shè)四個(gè)節(jié)點(diǎn)的溫度分別是Tl,T2,T3,T4,則單元內(nèi)任意一點(diǎn)的溫度 T可表達(dá)為其中W1, W2, W3, W4,為插值函數(shù)，也稱為形函數(shù)，定義如下:V2 =4（l + Xl-T|）4明二。十自火斗士4% -g + n） 4其中己

7、和“稱為形坐標(biāo)，取值區(qū)間為-1,1基于式（7）,對(duì)溫度的偏導(dǎo)數(shù)可寫為：T?n將此二偏導(dǎo)數(shù)代入弱化表達(dá)式（6）,該式就轉(zhuǎn)化為以節(jié)點(diǎn)溫度為變量的代數(shù)方 程：+ 0- d+ 忡，qdQ = 0到此，為得到原偏微分方程的有限元表達(dá)，我們只需將任意標(biāo)量函數(shù)，依次取為四個(gè)插值函數(shù)，W1-4,就得到四個(gè)代數(shù)方程：注意到式（9）-（12）中下標(biāo)的規(guī)律，我們可將這四個(gè)方程合并寫成矩陣的形式:使用下標(biāo)表達(dá)，式（13）可進(jìn)一步縮寫為:式中對(duì)應(yīng)于下標(biāo)i = 14的每一個(gè)值，下標(biāo)j取遍這。式（14）就是熱傳導(dǎo)偏微分方程（1）的四節(jié)點(diǎn)雙線性有限元表達(dá)，也是我們接 下來編寫有限元程序的出發(fā)點(diǎn)。使用高斯消去法解線性方程組的二

8、維熱傳導(dǎo)有限元程序這個(gè)程序是專為解一個(gè)特殊的熱傳導(dǎo)問題而設(shè)計(jì)的。所解問題是：已知一個(gè)無限 長(zhǎng)圓筒，內(nèi)徑100毫米，外徑200毫米，筒內(nèi)壁表面溫度1,外壁絕熱，求 圓筒截面上的溫度分布。圓筒材料各向均質(zhì)，熱傳導(dǎo)系數(shù)為 1 （單位還得查一 下，但也無所謂）。問題的解答很簡(jiǎn)單：均布，截面各點(diǎn)溫度都是10為什么？因?yàn)橥獠拷^熱，沒有熱量損失。溫度只能是均布。而內(nèi)壁溫度為1C,所以到處都是10-200Q11 f潮 .1009 Q因?yàn)閱栴}的幾何圖形簡(jiǎn)單，有限元網(wǎng)格便容易自動(dòng)生成。在以下程序中，第 12 行至第51行，生成四節(jié)點(diǎn)單元的有限元網(wǎng)格?？刂谱兞坑袃蓚€(gè)，Cdiv和Tdivo 前者定義沿圓周分成多少單元

9、；后者定義沿徑向即筒壁厚度方向分成多少單元。如果Cdiv = 8, Tdiv = 2,則所生成有限元網(wǎng)格如下（由第52行子程序DrawFEM 畫出）：2Q。*二 中-伽- 穴、-/ 、一11QQ-/?、/ I/ / /sn 4 / j J / z/ / 0 tw -5Q 人t _第64行使用MATLAB命令syms定義了兩個(gè)符號(hào)變量 ksi（即前面公式中的 E）,eta（ 0 n 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ì)算這些

10、導(dǎo)數(shù)在高斯積分點(diǎn)的數(shù)值。本程序中，每個(gè)單元有四個(gè)高斯積分點(diǎn)，也就是說，在ksi和eta兩個(gè)方向上，各有二個(gè)高斯積分點(diǎn)。第101至124行，根據(jù)式（14）計(jì)算單元?jiǎng)哦染仃?，Kelem ,并將其裝配入總 勁度矩陣Ko本問題沒有熱源，所以在單元循環(huán)水平上，不需計(jì)算式（14）中 的熱源積分項(xiàng)。第127至139行，施加邊界條件。圓筒外壁是絕熱條件，即法向熱流等于零。 在有限元中，這是自然滿足的，所以式（14）中的邊界熱流積分項(xiàng)為零，不必 考慮。唯一需要考慮的，是圓筒內(nèi)壁溫度等于 1。這種溫度給定的邊界條件， 在數(shù)學(xué)上叫第一類邊界條件。在有限元技術(shù)中，有各種各樣的方法施加第一類邊 界條件。主要是考慮提高內(nèi)

11、存效率。鑒于本程序的目的是進(jìn)一步開發(fā)波前法，不 需要仔細(xì)考慮如何更有效地施加這種溫度給定的邊界條件，因而所用的方法是最 簡(jiǎn)單的一種：即將內(nèi)壁邊界節(jié)點(diǎn)的各行方程，全部換為T = Tin （Tin = 1 C）。相應(yīng)地，將這些行的主元素所占據(jù)的列，乘以Tin后，移到等號(hào)右邊。這樣處理后，就得到待解的線性方程組：K x = F。第141行，使用本人自編的高斯消去法函數(shù)，egauss,解出為未知量x,也就 是個(gè)節(jié)點(diǎn)的溫度，者B等于1。這一行，也可直接用MATLAB命令，x = K F,取 代。我使用egauss的目的，是為下一步與波前法解方程比較效率。程序一：使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限

12、元程序1 -12% Purpoie:3ISolution of the heat conduction in a thick cyliner with FEM.41 Calling special funtions：6%dravFEILa for dravini finite tleaent aeih6%cgausa.b for the aolution of the systen of linear equations.7% Author： Yanshenf Jiang8% May 20, 2007Q%一11710%11clear all;12Cdiv = inputCniaber of d

13、ivisions in the circimferencial direction - ,);13Tdiv = inputCnimber of divisions in the thickneM direction =)；14%15tie16a 1 100: Wnntr rtdiuii ,17b - 200; MXiter radiust an18NumcroDeElcmcntos = Cdiv * Tdiv;19NumeroDeNos 二(Cdiv)*(Tdiv+l);20% Coordonnedades21X = zeros(NumeroDeNos, 1);22Y = zeros(Nume

14、roDeNos, 1);23dalpha = 2*pi/Cdiv;24dT= (b-a)/Tdiv;25delta n = Tdiv+1;26m = delta_n (Cdiv - 1);27for i = l：TdivH28r = a +29nl = i;30nfin = m + nl：31X(nl:delta_n：nfin) = r*cos(0：dalpha：2*pi-dalpha);32Y(nl:delta_n:nfin) = r*sin(0:dalpha:2*piwdalpha);33end34% Connect ividade de elementos35TCO = zeros(Nu

15、meroDeElementos, 4),36for i = l:Tdiv37nel = (i-l)*Cdiv 1;38nl = i;39n2 = nl + 1;40n3 = n2 + delta_n;41n4 = n3 - 1：42ICO(nel, 1:4) = nl n2 n3 n4：-13ne = nel;44for j = 2:Cdiv-l4Sne = ne + 1;46IC0(ne1:4) = IC0(ne-l, 1:4) + delta-n;47end48nefin = ne + 1：49IC0(nefin,l：2) = ICO(ne, 4) IC0(ne, 3)1;50ICO(ne

16、fin, 3:4) = Ln2 nlj,IFend52,Draw finite element mesh53drawFEM(NumeroDeElemento8i X Y,ICO);54%bb b657 Muterial datakconx - k、Coefficient of hot c.6b 66 % n4 oo n3w/%oo70%nln27172% psi = 0. 2b*( (I kai) *(l*etA)1-ksi) nunod-j nvod._pt j.); for i = l：nun causspoints_per ele97dpisileta_ma(i ：)iGuBPoint(

17、i I) GuPoint(， 2);99end100%101for ttlto = 1 :匕口舊總口?？贓lt：;呢ntnc102，numeros global 3 dom nos do 9lesiento103Noa = ICO(eleaentoi :);104xe = X(ob);105ye - Y(Nob)；106Ke lets2eros (nuM-nodeS-per.elem, num_nodes_per_elea).107for igauss - 1 ：nun_saixspoint8_j)er_eleii108JacohiMn109J - Ldpsidksi_niuD(iBaussi

18、 :1*xe ,Noe) = K(NoarXos) + Kelex;；123clear J invj detJ xe ye Kelen;124end125 Applv essential boundarv conditions usine Reddys alcorithm, p. 173126% boudary conditions： temperature at the inner surface of cylinder = 1G127Tin = i-128, normal thermal flux dl/dn 0 at the outer aurafee.129% nodes on the

19、 inner surface *here teger&tur& Tin is applied.130ebedof 11:deltH n：m+1j,131for i = 1:length(ebedof)132F = F - tKebcdoffi) * Tin;133end134F(ebcdof)工 Tin:136Ktebcdof, :) = 0;136K(： ebedof) = 0;137for i - l:lenith(eWof)133Kebedo(f (i) t ebedof () = 1;13gend140% Solution for nodal teoiperature141X =日耳a

20、usE(K, F);：4toe波前法的基本概念與算法有限元形成的線性方程組的系數(shù)矩陣，即剛度矩陣，是稀疏矩陣，也就是說，矩 陣?yán)锖性S許多多的零元素。有限元網(wǎng)格節(jié)點(diǎn)數(shù)目越巨大，非零元素與零元素的 比值越小，剛度矩陣越稀疏。用普通高斯消去法求解這樣的線性方程組，完全不 考慮矩陣的稀疏性，對(duì)大量的零元素也進(jìn)行加減乘除，結(jié)果浪費(fèi)了大量時(shí)間。不 僅如此，應(yīng)用高斯消去法，因?yàn)樾枰獙⒄麄€(gè)剛度矩陣存在計(jì)算機(jī)內(nèi)存里，所需計(jì) 算機(jī)內(nèi)存量與有限元網(wǎng)格節(jié)點(diǎn)未知量總數(shù)成平方的關(guān)系。以熱傳導(dǎo)問題為例。一 千個(gè)節(jié)點(diǎn)，光存剛度矩陣，就需內(nèi)存 1000 x 1000 x 8 /(1024 x 1024 )= 7.8Mb。這還

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

22、需要花錢買的有限元程序，就 會(huì)發(fā)現(xiàn)，解算同樣的問題，只需幾分鐘。您甚至可以毫無困難地解算十萬個(gè)節(jié)點(diǎn) 的熱傳導(dǎo)問題。秘密就在于，這些商業(yè)有限元軟件，在求解線性方程組時(shí)，都盡可能地利用剛度 矩陣的稀疏性。波前法就是這樣一種充分考慮了剛度矩陣的稀疏性求解方程的方 法。剛度矩陣為什么會(huì)稀疏？因?yàn)樵谟邢拊?，一個(gè)節(jié)點(diǎn)的影響，只限于它所連接的那些單元。為什么？就是 因?yàn)樵谇懊妫覀兺茖?dǎo)有限元時(shí)，在式(2)中，將熱傳導(dǎo)偏微分方程乘以的那 個(gè)神奇函數(shù)。我們說過，是任意標(biāo)量函數(shù)。既然是任意的，當(dāng)然可以任意選 取。然而我們沒有任意”地、胡亂地選取，而是居心叵測(cè)地，讓它恰恰等于有限 元的插值函數(shù)！而這些插值函數(shù)，恰

23、恰只在給定單元內(nèi)部非零，在單元外邊一律 為零!換句話說，插信函數(shù)只是些局部函數(shù)。我們讓等于插值函數(shù)，它也就具有了這種局部性。正是 的這種局部性，使得一個(gè)節(jié)點(diǎn)的影響，只限于它所 連接的單元。有限元方法，之所以能夠在計(jì)算力學(xué)領(lǐng)域里令人眼花繚亂的各種各 樣的計(jì)算方法中，獨(dú)領(lǐng)風(fēng)騷，傲視群雄，鶴立雞群，至今幾達(dá)50年之久，其全部奧妙，皆出于此！為進(jìn)一步考察這些影響到底有多少，我們來看下面的例子。使用前面的MATLAB 有限元程序，給 Cdiv的值輸入8, Tdiv輸入2 ,就會(huì)生成以下網(wǎng)格。它將圓 周分成8份，厚度分成2份。21圖中括弧內(nèi)是單元號(hào)碼，其余數(shù)字為節(jié)點(diǎn)號(hào)碼?？梢钥闯?，第1節(jié)點(diǎn)只與第1和第8單

24、元相連，其影響也就只限于這兩個(gè)單元。這里所說的影響，就是在剛度 矩陣中，第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ù)等于24o在第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é)

25、點(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è)元素為零。如此，可想而知，總剛度矩陣的每一行 的大部分元素都是零?，F(xiàn)在我們要考慮怎樣利用這些零元素了。在以上網(wǎng)格中，共有16個(gè)單元，24個(gè)節(jié)點(diǎn)。使用高斯消去法，會(huì)生成24 x 24的總剛度矩陣，即有24行，24歹限而 使用波前法，總剛度矩陣的行數(shù)雖然不變，也是24,但列數(shù)僅為11 （至于為什么是11 ,下面要詳細(xì)討論）。最重要的，在本例情況下，是波前法根本就不 需要將24 x 11的總剛度矩陣存在內(nèi)存中，而是存在硬盤上的。內(nèi)存

26、里，波前 法只需要放一個(gè)11 x 11的所謂波前矩陣就行。那么，什么是波前矩陣呢?就是在某一時(shí)刻，彼此發(fā)生關(guān)系的節(jié)點(diǎn)的剛度系數(shù)組成的矩陣。這個(gè)矩陣是方 的，其中含有最多非零元素的那一行就叫波前。什么叫某一時(shí)刻？就是某一單元！如前面 MATLAB程序所示，計(jì)算有限元?jiǎng)偠?矩陣有兩重循環(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ù)就可以從波前矩陣中清除掉，以便空

27、出位置來，存儲(chǔ)其它節(jié)點(diǎn)信息。因此，一個(gè)節(jié)點(diǎn)可以被消元的時(shí)間，是可以用單元循環(huán)的進(jìn)度來度量的。那么，一個(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)所連接

28、的最大單元編號(hào)。程序如下：程序二：計(jì)算節(jié)點(diǎn)消元時(shí)間1mJi|i*fTilIH=函口加hiwiMIeksJJ；2for i = uDKrolkEkmentiK3、m=KOiih：；45fnip j = 1: XhNwIics6CnnwecteclMl? -1if i X iMk, 1ifTiinci m t Leri cd h chJ l* l8NodeDeactiTtTimConiHTtedNndc - i;9mlKriidnr nd以上程序中，第三行的ICO變量是個(gè)兩維數(shù)組，18行，4歹I。它的每一行代表 一個(gè)單元，該行的4列給出該單元的四個(gè)節(jié)點(diǎn)。這段程序執(zhí)行的結(jié)果，是在一維 數(shù)組變量Nod

29、eDeactiveTime中定義每個(gè)節(jié)點(diǎn)可以消元的時(shí)間 （即單元號(hào)）。此時(shí)，NodeDeactiveTime 的值是：8161621010311114121251313614147151581616。第1個(gè)數(shù)“8弋表第1節(jié)點(diǎn)的消元時(shí)間是8 （單元）；第2個(gè)數(shù)“1騏表第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ù)i41rllll riHnl jibfciX

30、iiiiix 加3小i：2ni id imil Mill J；L而的M/l;DK，*l 小曲山1：6uudMHijilih * m,ii|Mjih VhIlhd山川 im卜闡岫;1尿卜,仙。麻Bj= l:hXndesT(g ttckddr = gjl;10Lkhiih - Xi J .klhiihjin lijiinJi11iliL JdiM12 HeMmdrd、* th13diifeadlbtnd16i LdtJrhd imik171Hi八Lui itin.nl ： MitUidihjiuxI Widlh ；在以上程序中，第1行開了一維輔助數(shù)組，NodeInFront ,用于記錄每一個(gè)節(jié)點(diǎn)

31、是不是在波前矩陣中，俵示在，”欣示不在。第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。”但這沒有關(guān)系，現(xiàn)在只是重新再植一次“1,”再一次表示該節(jié)點(diǎn)在波前矩陣中。第 6行計(jì)算 maxFrontWidth 。就是將NodeInFront中的“1” 相加，再與當(dāng)前 maxFrontW

32、idth 比較，誰的值更大，maxFrontWidth 就取誰 的值。也就是說，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ù)。前面說的波前法總剛度矩陣是 24行x 11歹I。其中的11,就是這 樣得出的。它的含義，就是在整個(gè)計(jì)算過程中，某一

33、時(shí)刻，同時(shí)存在于波前矩陣的節(jié)點(diǎn)數(shù)，其值最大為11。以上程序，實(shí)際上模擬了節(jié)點(diǎn)進(jìn)入和離開波前矩陣的裝配和消元過程。下表給出 波前矩陣中的節(jié)點(diǎn)隨著單元循環(huán)的整個(gè)變化過程。第 1列是單元號(hào)碼。第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歹I。可以看出，兩次達(dá)到最大波前，分別當(dāng)循環(huán)到第7單元和第10單元時(shí)。單元波前期庫中的節(jié)點(diǎn)波前 七點(diǎn) 融消元節(jié)點(diǎn)消元后波域節(jié)點(diǎn)數(shù)11544無七交化2i254X64w3i2511X710776A2511K14in13

34、Bin75i25U6iH1417印161910lbq?i5)13M17加23121110*1一5nM14J7帛22to12. 1K9X25TiUIT-sr2V何1無無如化103nXE41720n(i 911H3卞tX14J7歲M1”8 *N12312H15U179IL 1271332IH1?14172*114. IS(|14于 IH217MJU7門.IX5153U21202462IL 2J42與43 lt 24* J3U讀到這里，對(duì)照前面那個(gè)有限元網(wǎng)格，您要是能夠理解這個(gè)表中所有數(shù)字的來龍 去脈，您就理解了波前法！剩下的，是一些技術(shù)上的細(xì)節(jié)，可以通過閱讀下列全 部程序

35、來最后徹底 摳”懂。這些技術(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)偠?矩陣，我們可以通過ICO表，查到該單元4個(gè)節(jié)點(diǎn)的總編號(hào)。而節(jié)點(diǎn)總編號(hào)與 總剛度矩陣K的行與列之間

36、具有一一對(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里的零元素），可供黛玉葬花，姐妹們賞月吟詩。波前法的裝配算法則要復(fù)雜得多。因?yàn)椴ㄇ熬仃囆。荒芡瑫r(shí)裝下所有節(jié)點(diǎn)的所 有元素。這種情況，就如同旅店。世界上所有旅店是住不下裝得下世界上所有居 民的，旅客必須有進(jìn)有出才行。節(jié)

37、點(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ì)被不同旅客先后住過。所以，波前矩陣的行與列，與總節(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)還沒進(jìn)入過波前矩陣，那么就在波前矩陣給它分配一個(gè)自

38、由的行與列。 在下面的程序中，我們使用一維數(shù)組 freeLines來記錄波前矩陣每一行（同時(shí)也 是每一列。我們假定行數(shù)等于列數(shù)，也就是說，一個(gè)節(jié)點(diǎn)占據(jù)了第n行，它也就占據(jù)了第n歹I）的狀態(tài)：0表示自由；大于0表示已占用（即占用該行的節(jié) 點(diǎn)）。我們還使用一維數(shù)組GlobalID2FrontId 來記錄每一節(jié)點(diǎn)在波前矩陣占據(jù) 的行數(shù)（二列數(shù)）。也就是說，給出一個(gè)節(jié)點(diǎn)的總編號(hào),就能在GlobalID2FrontId 中查到它在波前矩陣中的位置。這個(gè)位置如果是零，就表示該節(jié)點(diǎn)還沒進(jìn)入過波前矩陣。這些算法的具體實(shí)施可見下列程序中的第148-168行。程序四：使用波前法解線性方程組的二維熱傳導(dǎo)有限元程1t

39、:2，PtnpiJMT：3、 Solution of the hr*t conduction in a thick cyliiier with FEM and4“a frun tai nitlhvd soher for (Ik4 1 Liiiear equiilio恒5， ilhnghnitiuiis；6%dnni rEM-ni 1or drjwing LiihIl1 tknicnl nn sii7f;15Cdii 二 in|nit ntiniber, iTdiviiniis in Iht cirriinlien-rinal diivrlion = ：J6end17Tdiv = inputsn

40、um her of divisions in the radial direction = 111Sh hik Tdh 119dkpJWn siuall number ff dikinns in thi I hkIkJ direct20Tdi* = inputrnumber ; ft JnrxT ludiust mni21l = 200： rOuler radiuMUinZ5uiiirL n|elJcnkvnt(s =(dh : I ftiv；26NmiHToDcNcft = iCdi、)*Tdiv+l)；27X = zeroMMumroDeNo 1);28Y zrroYgmeroDeNosd

41、);dalphj - 2*pi/Cdiv：30 |nl n2 a3n4;48nc = IK- 1:49forj = 2:Cdiv-I|SOne = nc + 1;51ICO( ned:4) = ICO, n。1.1:4) + deha_n;52rnd53nd in = ne + 1;54ICOncfin,l:2) = ICOncdc = l:dclkcn:m-1-l;63BCcodednner suiTaw nodes) = essential;64R Cy al nd i n ncr.Mi rfact-n xics) = Tin:6bkconx =1:66kcony = kconx;67sy

42、ms bi eta;68psi = 0.25|*( l-eta1 l4-ksi)* I+etal (l-ksi)a1 1 retail：69dpsidksi = dim pMksi);170dpsidcki =(IifTshvllcs_|HT.ckiii = 4:7879num gausspoints |)cr elem = 4;dpsidksLnum = ztroM ntiin_jiausspiiits_|M r_clcm. inim_n(Hlc-s_|KT_clcmi：80dpsideta num = zvroMnum gausspoints per elenu num nodes per

43、 ekin)：81for i = l:num_jiaiiss|Hints_prr_elcin82dpsidksijHiiniL：)=doublc(subs(dp&idksi4kM(ahGaussPoinU(i.l)tGausftPointMi2n):83dpsideta numli.:)二.doubldsubsi dpsideta, ksiArtaJJGaussPointM i.ll.GaussPointsl84end85maElvmPerde = 4;6687NodeDeactiveliine = zt ros(Xumcrul)vA,l h for i = 1 :NunieroDe Elements88Noss ICO(i:):I89bodes = length! Nos);0for j = l:NbNodcs91Cmnrtrl4Mk = X XodvDcactix vU inici Conncctvd