有限元二維熱傳導(dǎo)波前法MATLAB程序

1、有限元二維熱傳導(dǎo)波前法MATLAB程序 HYPERLINK /blog/ l m=0&t=1&c=fks_087069093081084064086085086095085085082075087086084074093 科技論文 2009-12-28 02:40:17 閱讀226 評論15 字號：大中小 江巖聲二維熱傳導(dǎo)有限元使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序波前法的基本概念與算法使用波前法解線性方程組的二維熱傳導(dǎo)有限元程序消元過程波前法與高斯消去法的效率之比較小結(jié)：波前法的過去、現(xiàn)在和未來波前法是求解線性方程組的一種方法，廣泛用于有限元程序。它最初由英國人（？）B.M. Ir

2、ons于1970在“國際工程計(jì)算方法雜志”上發(fā)表。30多年來，波前法有了不少變種。本文所用算法，采于法國人Pascal JOLY所著Mise en Oeuvre de la Mthode des Elments Finis。這本書是我1993年在比利時(shí)一家書店買的，書中有一節(jié)波前法，六頁紙，解釋了基本概念和算法，但沒有程序，也沒有細(xì)節(jié)討論。我曾花了兩個(gè)半天的時(shí)間，在網(wǎng)上尋找波前法程序，或更詳細(xì)的資料，沒有找到（需要花錢才能看的文獻(xiàn)不算）。倒是看到不少中國人，也在尋找。一些人說，波前法程序太難懂了。通過自己編寫程序，我同意這些人的說法，確實(shí)難。我還真很少編如此耗費(fèi)腦力的程序。完工之后，我曾對朋友

3、老王說，有了算法，編程序還這么難，當(dāng)初想出算法的人，真是了不起。現(xiàn)將我對波前法的理解和編程體會(huì)解說如下，供感興趣的網(wǎng)友參考，也為填補(bǔ)網(wǎng)絡(luò)上波前法空白。二維熱傳導(dǎo)有限元波前法和有限元密不可分。因而，在編寫波前法程序之前，必須有個(gè)有限元程序。為了簡化問題，最好是能解算一個(gè)節(jié)點(diǎn)上只有一個(gè)自由度的問題的有限元程序，而且要盡可能地簡單。我手邊現(xiàn)有的有限元程序都不符合這個(gè)要求。就決定先開發(fā)一個(gè)解算二維熱傳導(dǎo)問題的MATLAB有限元程序。二維熱傳導(dǎo)問題的微分方程是其中 T 是溫度，Kx 和 Ky 分別是 x 和 y 方向上的熱傳導(dǎo)系數(shù)，q 是熱源。對于這樣的比較經(jīng)典的二階微分方程，如何導(dǎo)出有限元表達(dá)式？這個(gè)

4、問題，是有限元的首要問題！我相信，許許多多學(xué)過有限元的人，甚至每天都在用有限元的人，并不真的十分明了。我自己曾經(jīng)就是這樣。在我于2005年3月到巴西教書之前，我搞過20年有限元，其中有六年，還是在比利時(shí)一個(gè)專門開發(fā)有限元程序的公司里度過的。但是，如果您那時(shí)問我，任給一個(gè)二階偏微分方程，例如上述熱傳導(dǎo)方程，如何導(dǎo)出有限元表達(dá)？說老實(shí)話，不看書，我還真就不會(huì)！直到2006年，我教了一遍有限元后，才弄明白（惟教書才是最好的學(xué)習(xí)）。其實(shí)簡單極了：只需將那二階偏微分方程，乘上一個(gè)任意標(biāo)量函數(shù)，然后在某個(gè)有限單元上積分。請看下列推導(dǎo)：即其中， e 是單元面積； 是任意標(biāo)量函數(shù)。 注意在以上積分中，溫度要遭

5、受二階偏導(dǎo)，這很不好。有限元的精髓，在于通過分步積分，將其中一階偏導(dǎo)轉(zhuǎn)移到那個(gè)任意標(biāo)量函數(shù) 上，這樣就可降低一階溫度偏導(dǎo)，改善對它的苛刻待遇。這得用到您在高等數(shù)學(xué)最后一章里學(xué)過的散度定理（Theorem of divergence）：其中， 是面積的邊界； （反） 是梯度算子；F 和 G 是任意兩個(gè)矢量函數(shù)。 對于二維問題，上述散度定理可寫為：將此散度定理應(yīng)用于（2）式中的第一項(xiàng)積分，令得：將此積分結(jié)果帶入（3）式，得到熱傳導(dǎo)偏微分方程的弱化表達(dá)（weak form）：所謂“弱化”，是指對溫度函數(shù)的可導(dǎo)階數(shù)要求降低了。在原熱傳導(dǎo)偏微分方程（1）中，溫度函數(shù)必須是二階可偏導(dǎo)函數(shù)，而在弱化表達(dá)（6

6、）中，溫度函數(shù)只要一階偏導(dǎo)就行了。有限單元法就是以偏微分方程的弱化表達(dá)式為出發(fā)點(diǎn)的?，F(xiàn)在，到了說明那個(gè)任意標(biāo)量函數(shù)，是何方神圣的時(shí)候了。有限單元法有各種各樣的變種，而最常見的，應(yīng)用最廣泛的，是所謂迦遼金法（Galerkin），即將這個(gè)任意標(biāo)量函數(shù)，等同于，有限元的插值函數(shù)。迦遼金法的優(yōu)點(diǎn)是可以最終形成對稱勁度矩陣，從而具有較好的數(shù)值穩(wěn)定性。我們知道，在一個(gè)有限單元中，任意一點(diǎn)的值，例如溫度，是用節(jié)點(diǎn)上的溫度表達(dá)的。對于一個(gè)四節(jié)點(diǎn)雙線性單元來說，設(shè)四個(gè)節(jié)點(diǎn)的溫度分別是T1,T2,T3,T4，則單元內(nèi)任意一點(diǎn)的溫度 T 可表達(dá)為其中 1，2，3，4，為插值函數(shù)，也稱為形函數(shù)，定義如下：其中 和 稱

7、為形坐標(biāo)，取值區(qū)間為 -1,1?；谑剑?），對溫度的偏導(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è)方程合并寫成矩陣的形式： HYPERLINK /photo/2We1ncQaz19_4yqE6j_iYA=/5398127103357919631.jpg t _blank 使用下標(biāo)表達(dá)，式（13）可進(jìn)一步縮寫為：式中對應(yīng)于下標(biāo) i = 14 的每一個(gè)值，下標(biāo) j 取遍14。式（14）就是熱傳

8、導(dǎo)偏微分方程（1）的四節(jié)點(diǎn)雙線性有限元表達(dá)，也是我們接下來編寫有限元程序的出發(fā)點(diǎn)。 使用高斯消去法解線性方程組的二維熱傳導(dǎo)有限元程序這個(gè)程序是專為解一個(gè)特殊的熱傳導(dǎo)問題而設(shè)計(jì)的。所解問題是：已知一個(gè)無限長圓筒，內(nèi)徑100毫米，外徑200毫米，筒內(nèi)壁表面溫度 1C，外壁絕熱，求圓筒截面上的溫度分布。圓筒材料各向均質(zhì)，熱傳導(dǎo)系數(shù)為 1 （單位還得查一下，但也無所謂）。問題的解答很簡單：均布，截面各點(diǎn)溫度都是 1C。為什么？因?yàn)橥獠拷^熱，沒有熱量損失。溫度只能是均布。而內(nèi)壁溫度為 1C，所以到處都是1C。因?yàn)閱栴}的幾何圖形簡單，有限元網(wǎng)格便容易自動(dòng)生成。在以下程序中，第12行至第51行，生成四節(jié)點(diǎn)單

9、元的有限元網(wǎng)格?？刂谱兞坑袃蓚€(gè)，Cdiv 和 Tdiv。 前者定義沿圓周分成多少單元；后者定義沿徑向即筒壁厚度方向分成多少單元。如果 Cdiv = 8，Tdiv = 2， 則所生成有限元網(wǎng)格如下（由第52行子程序DrawFEM畫出）：第64行使用MATLAB命令 syms 定義了兩個(gè)符號變量 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ù)值。本程序中，每個(gè)單元有四個(gè)

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

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

12、組的系數(shù)矩陣，即剛度矩陣，是稀疏矩陣，也就是說，矩陣?yán)锖性S許多多的零元素。有限元網(wǎng)格節(jié)點(diǎn)數(shù)目越巨大，非零元素與零元素的比值越小，剛度矩陣越稀疏。用普通高斯消去法求解這樣的線性方程組，完全不考慮矩陣的稀疏性，對大量的零元素也進(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.8 Mb。 這還沒問題。但若要計(jì)算一萬個(gè)節(jié)點(diǎn)的問題，就需要 10 x 10 x 7.8 = 78

13、0 Mb 來存剛度矩陣。內(nèi)存為 1 Gb的計(jì)算機(jī)已經(jīng)不能計(jì)算這樣的問題了，因?yàn)槲④浺暣暗绕渌到y(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)為是無限長！在這種情況下，因?yàn)槠胀ǜ咚瓜シㄐ枰獦O其頻繁地使用虛擬內(nèi)存，它的解算時(shí)間也就無限地延長了。但如果您在這樣的計(jì)算機(jī)上運(yùn)行ANSYS，或任何需要花錢買的有限元程序，就會(huì)發(fā)現(xiàn)，解算同樣的問題，只需幾分鐘。您甚至可以毫無困難地解算十萬

14、個(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)然可以任意選取。然而我們沒有“任意”地、胡亂地選取，而是居心叵測地，讓它恰恰等于有限元的插值函數(shù)！而這些插值函數(shù)，恰恰只在給定單元內(nèi)部非零，在單元外邊一律為零！換句話說，插值函數(shù)只是些局部函數(shù)。我們讓等于插值函數(shù)，它也就

15、具有了這種局部性。正是的這種局部性，使得一個(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份。圖中括弧內(nèi)是單元號碼，其余數(shù)字為節(jié)點(diǎn)號碼?？梢钥闯?，第1節(jié)點(diǎn)只與第1和第8單元相連，其影響也就只限于這兩個(gè)單元。這里所說的影響，就是在剛度矩陣中，第1和第8單元的所有節(jié)點(diǎn)，即第1、2、5、

16、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個(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

17、、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)存里，波前法只需要放一個(gè) 11 x 11 的所謂波前矩陣就行。那么，什么是波前矩陣呢？就是在某一時(shí)刻，彼此發(fā)生關(guān)系的節(jié)點(diǎn)的剛度系數(shù)組成

18、的矩陣。這個(gè)矩陣是方的，其中含有最多非零元素的那一行就叫波前。什么叫某一時(shí)刻？就是某一單元！如前面MATLAB程序所示，計(jì)算有限元?jiǎng)偠染仃囉袃芍匮h(huán)，最外面那層循環(huán)，是對單元循環(huán)，即從編號為第一的單元，到編號最大的單元。在波前法中，當(dāng)循環(huán)到某一單元時(shí)，在計(jì)算該單元?jiǎng)偠染仃囈院?，還要看看哪一個(gè)節(jié)點(diǎn)可以消去了，也就是消元。被消元的節(jié)點(diǎn)，對以后其它單元?jiǎng)偠染仃嚲筒辉儆杏绊懥?，該?jié)點(diǎn)的剛度系數(shù)就可以存入硬盤指定文件中，而這些系數(shù)就可以從波前矩陣中清除掉，以便空出位置來，存儲(chǔ)其它節(jié)點(diǎn)信息。因此，一個(gè)節(jié)點(diǎn)可以被消元的時(shí)間，是可以用單元循環(huán)的進(jìn)度來度量的。那么，一個(gè)節(jié)點(diǎn)，什么時(shí)候可以消元了？就是與該節(jié)點(diǎn)相連

19、的所有單元都循環(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)便可以被消元了。算法很簡單，就是查找每一個(gè)節(jié)點(diǎn)所連接的最大單元編號。程序如下：程序二：計(jì)算節(jié)點(diǎn)消元時(shí)間以上程序中，第三行的ICO變量是個(gè)兩維數(shù)組，18行，4列。它的每一行代表一個(gè)單元，該行的4列給出該單元的四個(gè)

20、節(jié)點(diǎn)。這段程序執(zhí)行的結(jié)果，是在一維數(shù)組變量NodeDeactiveTime中定義每個(gè)節(jié)點(diǎn)可以消元的時(shí)間（即單元號）。此時(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（單元）；余類推。請注意，消元最早的節(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ì)算波前矩陣的最大階

21、數(shù)了。程序如下：程序三：計(jì)算波前矩陣的最大階數(shù)在以上程序中，第1行開了一維輔助數(shù)組，NodeInFront，用于記錄每一個(gè)節(jié)點(diǎn)是不是在波前矩陣中，“1”表示在，“0”表示不在。第2行將我們要計(jì)算的波前矩陣的最大階數(shù)maxFrontWidth的值初始為零。第3行開始對單元循環(huán)。對編號為 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)在波前矩陣中了，即它們在NodeInFront中的值已經(jīng)是“1”。但這沒有關(guān)系，現(xiàn)在只是重新再植一次“1”，再一次表示該節(jié)點(diǎn)

22、在波前矩陣中。第6行計(jì)算 maxFrontWidth。就是將NodeInFront中的“1”相加，再與當(dāng)前maxFrontWidth比較，誰的值更大，maxFrontWidth就取誰的值。也就是說，maxFrontWidth的值只增不減。第8行對 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ù)。前面說

23、的波前法總剛度矩陣是24 行 x 11列。其中的11，就是這樣得出的。它的含義，就是在整個(gè)計(jì)算過程中，某一時(shí)刻，同時(shí)存在于波前矩陣的節(jié)點(diǎn)數(shù)，其值最大為11。 以上程序，實(shí)際上模擬了節(jié)點(diǎn)進(jìn)入和離開波前矩陣的裝配和消元過程。下表給出波前矩陣中的節(jié)點(diǎn)隨著單元循環(huán)的整個(gè)變化過程。第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í)。讀到這里，對照前面那個(gè)有限元網(wǎng)格，您要是能夠理解這個(gè)

24、表中所有數(shù)字的來龍去脈，您就理解了波前法！剩下的，是一些技術(shù)上的細(xì)節(jié)，可以通過閱讀下列全部程序來最后徹底“摳”懂。 這些技術(shù)細(xì)節(jié)里，比較難懂的，是如何將單元?jiǎng)偠染仃囇b配入波前矩陣。單元?jiǎng)偠染仃囀莻€(gè)4階矩陣，即4x4矩陣，代表了單元4個(gè)節(jié)點(diǎn)之間的相互影響。例如，第1行里的4列元素，是單元4個(gè)節(jié)點(diǎn)對單元第1節(jié)點(diǎn)的影響；第2行里的4列元素，是那4個(gè)節(jié)點(diǎn)對單元第2節(jié)點(diǎn)的影響；余類推。每個(gè)單元節(jié)點(diǎn)的局部編號都是1、2、3、4，但它們的總節(jié)點(diǎn)編號是不同的，而且是唯一的。單元節(jié)點(diǎn)的局部編號與其總節(jié)點(diǎn)編號的對應(yīng)關(guān)系，在本文的程序中，由二維數(shù)組ICO給出。在用普通高斯消去法解方程時(shí)，任給一個(gè)單元?jiǎng)偠染仃?，我們?/p>

25、以通過ICO表，查到該單元4個(gè)節(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列。所以，使用普通高斯消去法解方程時(shí)，把單元?jiǎng)偠染仃囇b配進(jìn)總剛度矩陣的算法很簡單，如本文第一個(gè)程序的第122行所示，只一行程序便可實(shí)現(xiàn)。打個(gè)比方，這種情況下的總剛度矩陣 K，就像占地廣闊的大觀園，有許多房子，紅樓夢里的所有人等，例如金陵十二釵，各有各的住所，還有許多空地（相當(dāng)于K里的零元素），可供黛玉葬花，姐妹們賞月吟詩。 波前法的裝配算法則要復(fù)雜得多。因?yàn)椴ㄇ熬仃囆。荒芡瑫r(shí)裝下所有節(jié)點(diǎn)的所有元素。這種情況，就如同

26、旅店。世界上所有旅店是住不下裝得下世界上所有居民的，旅客必須有進(jìn)有出才行。節(jié)點(diǎn)進(jìn)出波前矩陣，就如同旅客進(jìn)出旅館。有的旅客住的時(shí)間長，有的短。節(jié)點(diǎn)在波前矩陣呆的時(shí)間也一樣，有長有短。隨著節(jié)點(diǎn)的進(jìn)進(jìn)出出，波前矩陣的每一行和列都可能先后被不同節(jié)點(diǎn)占用，就像旅館里的每個(gè)房間都會(huì)被不同旅客先后住過。所以，波前矩陣的行與列，與總節(jié)點(diǎn)編號之間，不再有像使用普通高斯消去法解方程時(shí)那樣的一一對應(yīng)的固定關(guān)系。單元矩陣的某一行、某一列的元素，應(yīng)該放到波前矩陣的哪一行、哪一列，是動(dòng)態(tài)分配的。有兩種可能：如果某節(jié)點(diǎn)已經(jīng)存在于波前矩陣中，那么就把該節(jié)點(diǎn)在單元矩陣中的元素加到波前矩陣的相應(yīng)元素上（因而需要知道它在哪里）；反

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

28、68行。程序四：使用波前法解線性方程組的二維熱傳導(dǎo)有限元程序消元過程消元也是波前法里比較難于理解的技術(shù)細(xì)節(jié)。我在調(diào)試波前法程序過程中，所花時(shí)間最多的，就是消元。而現(xiàn)在，當(dāng)我想解釋消元過程時(shí)，發(fā)現(xiàn)遇到了更大的困難。因?yàn)槲也恢缿?yīng)該從哪里說起。消元是個(gè)細(xì)活兒，必得有研究針尖上能站立多少天使的牛角尖精神才能弄懂。不是所有人都有這種精神。而為寫此文，已經(jīng)花了如此多工夫，幾乎超過了開發(fā)程序的時(shí)間！所以，還是算了吧。一簡對三繁，我就把消去第一個(gè)可消節(jié)點(diǎn)（第4節(jié)點(diǎn)）的過程列表于下。您若看得明白，就算徹底懂了波前法。為方便說明問題，我將原問題的內(nèi)壁邊界條件換為絕熱，外壁換為溫度給定（程序四第63至64行須作相

29、應(yīng)修改）。這樣，第4節(jié)點(diǎn)就不經(jīng)過邊界處理那一段程序（第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線對稱；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行（乘除兩個(gè)系數(shù)后），因而第4行

30、的元素全部出現(xiàn)在第1行中。消元的結(jié)果，是第4行主元素所在的列，除主元素外，其余元素全被消成零。此時(shí)，就可以將第4行存入指定文件中，如程序第203行所示。程序第205行將此時(shí)的波前節(jié)點(diǎn)號寫入另一指定文件。以后回代需要。程序第204行寫方程Ax = B的等號右端自由項(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列分別是使用程序四（

31、波前法）和程序一（高斯消去法）解題的時(shí)鐘時(shí)間（elapsed time）。由上表可看出，當(dāng)節(jié)點(diǎn)數(shù)很少（24和80）時(shí)，兩種方法算題的時(shí)間差不多，但當(dāng)節(jié)點(diǎn)數(shù)達(dá)到288時(shí)，差別開始顯示出來，高斯消去法需時(shí)比波前法將近多一倍。64 x 16的網(wǎng)格（1088節(jié)點(diǎn)），時(shí)間差別達(dá)到15倍。96 x 24 （2400節(jié)電），差別到了60倍。到了最后一個(gè)網(wǎng)格，128 x 32（4224節(jié)點(diǎn)），運(yùn)轉(zhuǎn)程序一（高斯消去法）需內(nèi)存563MB，這已經(jīng)超出了我的計(jì)算機(jī)內(nèi)存（512MB），因而所需時(shí)間可認(rèn)為無窮大。 小結(jié)：波前法的過去、現(xiàn)在和未來在人類計(jì)算技術(shù)發(fā)展史上，有限元方法可以說是個(gè)奇跡。從來沒有任何一種計(jì)算方法，能

32、像有限元這樣，深刻而廣泛地影響著人們的日常生活，它使得上百的人成為大師，上千的人成為專家，上萬的人靠它吃飯，上億的人每天乘坐使用有限元程序設(shè)計(jì)的飛機(jī)、汽車、輪船、火車。在有限元半個(gè)世紀(jì)的發(fā)展史當(dāng)中，無數(shù)人寫了無數(shù)關(guān)于有限元的文獻(xiàn)，可用汗牛充棟來形容，時(shí)至今日，已經(jīng)無人能夠全部讀完（也沒必要）。但1970年Irons發(fā)表的那篇關(guān)于波前法的文章，大概要算有限元文獻(xiàn)里被引用次數(shù)最多的一篇了。波前法的提出，突破了制約有限元應(yīng)用的瓶頸，使得當(dāng)時(shí)內(nèi)存量很小的計(jì)算機(jī)，可以解算規(guī)模很大的題目?？梢哉f，沒有波前法，就沒有有限元的大發(fā)展，也就沒有有限元的今天。但是，波前法的使用，也有上限。就是波前矩陣不能超過計(jì)算

33、機(jī)內(nèi)存。以我的計(jì)算機(jī)為例，516MB內(nèi)存。假設(shè)有一半可用來計(jì)算，就是約256MB。這相當(dāng)于256x1024x1024/8個(gè)雙精度實(shí)數(shù)。將此數(shù)開平方，得到約5700。也就是說，在我的計(jì)算機(jī)上，能使用波前法解算的題目，其波前矩陣的最大階數(shù)，約為5000。也就是說，如果是解算三維問題，其變量總數(shù)，大約是五、六萬。所以，今天，經(jīng)典波前法，例如本文程序四，您只能在大學(xué)里搞科研的那些教授自己開發(fā)的有限元恐龍化石里找得到了。今天，在像ANSYS那種商業(yè)有限元程序里，經(jīng)典波前法早已被波前法的變種所取代，就是多波前法。簡單地說，多波前法，就是將網(wǎng)格分成許多子區(qū)域，在每個(gè)子區(qū)域上使用波前法。這樣可以將波前矩陣的階

34、數(shù)降低，以突破內(nèi)存容許的解題上限并提高效率。每個(gè)子區(qū)域間數(shù)據(jù)的依賴關(guān)系與傳遞，是用所謂消去樹來管理的。毋庸置疑，這樣的多波前法，算法更加復(fù)雜，且編程量巨大，可能一萬行都嫌少！當(dāng)然，商業(yè)有限元程序里，不光是多波前法，還有其它方法；不光是直接解法，還有迭代法。這些程序，光解方程的部分，據(jù)介紹，有的已經(jīng)達(dá)到三萬行。但，萬變不離其宗。欲理解有限元程序，欲開發(fā)有限元程序，必須理解波前法。參考文獻(xiàn)：Pascal JOLY，Mise en Oeuvre de la Mthode des Elments Finis。J.N. REDDY, An Introduction to The Finite Eleme

35、nt Method, McGraw-Hill Book Company，1984周洪偉，吳舒，陳璞，有限元分析快速直接求解技術(shù)進(jìn)展，“力學(xué)進(jìn)展”雜志，第37卷，第2期，175-188頁，2007年5月25日。附錄資料：matlab繪圖指令大全繪圖指令1 二維曲線圖1.1 繪制折線圖plot指令圖例Y=1,3,6,5,9,0,2;plot(Y);X=0: pi/10: pi*2;Y=sin(X);plot(X,Y);X=0: pi/10: pi*2;Y1=sin(X);Y2=cos(X);Plot(X,Y1,X,Y2);調(diào)整坐標(biāo)范圍：axisaxis(0,300,0,2)1.2 繪制自定義函數(shù)D

36、rawCircle.mfunction DrawCircle(Point,Radius) Hold on t=0: pi/10: 2*pi; x=Point(1)+ Radius*cos(t); y=Point(2)+ Radius*sin(t); plot(x,y);DrawCircle(10 10,1)DrawCircle(20 10,2)DrawCircle(10 20,3)1.3 繪制符號函數(shù)顯函數(shù)ezplot(sin(x),0,2*pi)隱函數(shù)ezplot(x2+y2-10,-5,5,-6,6)參數(shù)方程ezplot(cos(t)3,sin(t)3,0,2*pi)1.4 繪制自定義函數(shù)

37、function y=myf1(x) y=sqrt(100-x2);fplot(myf1,-15 15)fplot(sin(x) cos(x) myf1(x),-15 15)1.5 圖形修飾 設(shè)置顏色 y m c r g b w k 設(shè)置線型 - : -. - 設(shè)置標(biāo)記 . o x + * 指令圖例Y=1,3,6,5,9,0,2;plot(Y, r-+);X=0: pi/10: pi*2;Y=sin(X);plot(X,Y, b-.);X=0: pi/10: pi*2;Y1=sin(X); Y2=cos(X);plot(X,Y1,r+-, X,Y2,b-*); 在指定坐標(biāo)處，書寫文字：text

38、(3.5, 0.6, 曲線比較);x=1.6*pi, 1.6*pi; y=-0.3, 0.8;s=曲線cos; 曲線sin; text(x,y,s);1.6 更多類型的二維圖指令圖例bar直方圖X=0:pi/10:2*pi;Y=sin(X);bar(X,Y);polar極坐標(biāo)圖T=0: pi/10: 4*pi;R=T;polar(T, R);誤差棒棒圖X=0:pi/10:2*pi;Y=sin(X);e=0.2*rand(size(X);errorbar(X,Y,e);火柴桿圖X=0:pi/10:2*pi; Y=sin(X);stem(X,Y);stairs樓梯圖X=0:pi/10:2*pi;

39、Y=sin(X);stairs(X,Y);多邊形填色圖X=1,2,3,4,5; Y=3,5,2,1,6;fill(X,Y,r);hold on; % 保持圖形plot(X,Y,o)1.7 數(shù)值函數(shù)的二維圖 可用于繪圖，更可用于采樣取點(diǎn)。 fplot(0.5*cos(x),-pi,pi) % 繪圖X,Y = fplot(0.5*cos(x),-pi,pi); % 返回點(diǎn)坐標(biāo)fplot(cos(x),-pi,pi,r-+); % 觀察點(diǎn)的位置控制采樣點(diǎn)的密度fplot(cos(x),-pi,pi,r-+,0.05);fplot(cos(x),-pi,pi,r-+,0.1); 可繪制系統(tǒng)函數(shù)，也可繪

40、制自定義函數(shù)的圖形。2 三維曲線圖2.1 三維曲線plot3指令圖例X=0: 0.1: 8*pi;Y=sin(X);Z=cos(X); plot3(X,Y,Z,r);X=0: 0.1: 8*pi;Y=sin(X);Z1=cos(X);Z2=2*cos(X); plot3(X,Y,Z1,r, X,Y,Z2,b);2.2 三維面填色fill3指令圖例X1=2,2,1;Y1=0,2,1; Z1=0,0,1;fill3(X1,Y1,Z1,r);hold on;X2=1,0,0;Y2=1,2,0;Z2=1,0,0;fill3(X2,Y2,Z2,r);X3=0,2,1;Y3=2,2,1;Z3=0,0,1;

41、fill3(X3,Y3,Z3,b);text(1,1,1,1,1,1);3 曲面圖形3.1 網(wǎng)格點(diǎn)坐標(biāo)的表示x=1:2:7y=2:2:6X,Y = meshgrid(x,y)X = 1 3 5 7 1 3 5 7 1 3 5 7Y = 2 2 2 2 4 4 4 4 6 6 6 63.2 三維網(wǎng)格mesh、meshc、meshz 用途：數(shù)據(jù)場的觀察分析命令圖例隨機(jī)數(shù)據(jù)的網(wǎng)格Z=rand(5,5);mesh(Z);% 設(shè)置顏色colormap(1,0,0);自定義函數(shù)的網(wǎng)格x=-4: 1: 4;y=-5: 1: 5;X,Y=meshgrid(x,y);Z=X.2+Y.2;mesh(X,Y,Z);

42、 消影開關(guān)：hidden on / hidden off 利用peaks(50)作為模擬數(shù)據(jù)矩陣;命令圖例 帶等高線的網(wǎng)格Z=peaks(50);meshc(Z);Z=peaks(50);meshc(Z);colormap(1,0,0);帶基準(zhǔn)面的網(wǎng)格Z=peaks(50);meshz(Z);剪孔Z=peaks(50);Z(30:45,15:30)=NaN*ones(16,16);meshc(Z);3.3 著色表面圖surf、surfc命令圖例表面著色的網(wǎng)格Z=peaks(50);surf(Z);自定義函數(shù)的著色網(wǎng)格x=-2: 0.1: 2;y=-2: 0.1: 2;X,Y=meshgrid(

43、x,y);Z=sqrt(X.2+Y.2);surfc(X,Y,Z);3.4 二元函數(shù)的偽彩色圖pcolor 用途：污染濃度場的觀察分析。命令圖例Z=peaks(50); pcolor(Z);colorbar(hor);colorbar(vec);3.5 等高線contour不僅可用于繪圖，更可以用以求截面數(shù)據(jù)。命令圖例以矩陣下標(biāo)為x、y分量的等高線Z=peaks(50);C=contour(Z);colormap(1,0,0);C：保存了全部等高線上的點(diǎn)坐標(biāo)。均分n條等高線，并標(biāo)注之Z=peaks(50);n=5;C=contour(Z,n);colormap(1,0,0);clabel(C)

44、;在指定高度繪制等高線Z=peaks(50);V=-10: 2: 10;C=contour(Z, V);colormap(1,0,0);clabel(C);完整圖形數(shù)據(jù)的等高線x=-2: 0.1: 2;y=-2: 0.1: 2;X,Y=meshgrid(x,y);Z=sqrt(X.2+Y.2);n=5;C=contour(X,Y,Z,n);三維等高線x=-2:0.1:2;y=-2:0.1:2;X,Y=meshgrid(x,y);Z=sqrt(X.2+Y.2);n=10;C=contour3(X,Y,Z,n); 3.6 矢量場圖quiver用于挖掘數(shù)據(jù)變化趨勢。命令圖例構(gòu)造起伏跌宕的曲面x=-2

45、:0.2:2;y=-1:0.2:1;X,Y=meshgrid(x,y);Z=X.*exp(-X.2-Y.2);mesh(X,Y,Z);xlabel(X軸);ylabel(Y軸);colormap(1,0,0);計(jì)算曲面的梯度px,py= gradient(Z,0.2,0.2);繪制矢量場圖quiver(x,y,px,py);3.7 視角控制view視點(diǎn)控制方式及效果：view(1 1 1)view(2 1 1)view(3 1 1)view(1 1 1)view(1 2 1)view(1 3 1)view(1 1 1)view(1 1 2)view(1 1 3)方位角、仰角控制方式及效果：缺省為(-37.5,30)。view(-37.5,30)view(-17.5,30)view(-5.5,30)view(-37.5,30)view(-37.5,45)view(-37.5,60)3.8 多視區(qū)控制subplotsubplot(2,1,1); mesh(X,Y,Z);subplot(2,1,2); quiver(x,y,px,py);3.9 制作、播放動(dòng)畫x,y,z=peaks(30); surf(x,y,z)%