流體力學(xué)課第十五章

1、流流 體體 力力 學(xué)學(xué) 退 出中國科學(xué)文化出版社目 錄流體力學(xué)基礎(chǔ)第一篇第二篇流體動力學(xué)基本原理及流體工程退 出第三篇計算流體動力學(xué) 第一篇 流體力學(xué)基礎(chǔ) 緒論 場論與正交曲線坐標(biāo) 流體靜力學(xué) 流體運動學(xué)第一章第二章第三章第四章退 出返 回 第二篇 流體動力學(xué)基本原理及流體工程 流體動力學(xué)微分形式基本方程 流體動力學(xué)積分形式基本方程 伯努利方程及其應(yīng)用 量綱分析和相似原理 流動阻力與管道計算 邊界層理論 流體繞過物體的流動 氣體動力學(xué)基礎(chǔ) 第五章第六章第七章第八章第九章退 出返 回第十章第十一章第十二章第三篇 計算流體動力學(xué) 計算流體動力學(xué)數(shù)學(xué)物理基礎(chǔ) 流體動力學(xué)問題的有限差分解法 流體動力學(xué)

2、問題的有限元解法 第十三章第十四章第十五章退 出返 回第十五章流體動力學(xué)問題的有限元解法 有限元法的基本思想與區(qū)域離散化 有限元法中代數(shù)方程的建立 二維邊值問題有限元法求解舉例 有限分析法介紹 第一節(jié)第二節(jié)退 出返 回第三節(jié)第四節(jié)第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第1頁頁 在上一章中，對求解流體動力學(xué)問題的有限差分方法進行了比較仔細(xì)的討論。有限差分法的優(yōu)點是原理簡單，便于實施，對非線性比較強的對流換熱問題有比較好的適應(yīng)性。其弱點是對復(fù)雜幾何形狀的區(qū)域適應(yīng)性比較差，采用近二十年發(fā)展起來的網(wǎng)格生成技術(shù)后可以克服這一弱點，但增加了計算工作量

3、。而有限元法對復(fù)雜幾何形狀的區(qū)域具有較強的適應(yīng)性。因而這兩種方法已在擴散方程的求解中得到廣泛的應(yīng)用，同時研究工作正在向流場求解方面深入展開。近二十年內(nèi)還發(fā)展起一種稱為有限分析法的數(shù)值方法。本節(jié)中將對有限元法和有限分析法的基本思想作簡要介紹，并以擴散方程為例說明其實施過程的主要步驟。 用有限元法求解物理問題時，總的解題步驟仍如圖13.3所示，它與有限差分法的區(qū)別主要在于區(qū)域離散化的方式不同，建立代數(shù)方程所依據(jù)的原則或方法不同，以及由此而引起的代數(shù)方程求解方法的不同。本節(jié)中主要討論有限元法的基本思想及其區(qū)域離散化（包括插值函數(shù)）等問題，下一節(jié)中再討論代數(shù)方程的生成。 退 出返 回 第一節(jié) 有限元法

4、的基本思想與區(qū)域離散化 第第2頁頁假設(shè)我們要在某一個區(qū)域 內(nèi)求解一個偏微分方程（例如二維穩(wěn)態(tài)的導(dǎo)熱方程）：圖15.1 三角形單元網(wǎng)格0)( fl)(l)(l2222yx123aa xa y一、有限元法的基本思想其中 表示一個微分算子，例如對直角坐標(biāo)中的二維穩(wěn)態(tài)導(dǎo)熱問題， 就代表 。首先把求解區(qū)域劃分成許多子區(qū)域，稱為單元。對二維問題，單元的形狀可以是矩形，三角形或四邊形，單元形狀的這種多樣性使有限元法對求解區(qū)域的幾何形狀有很好的適應(yīng)性。（15.1）在每個單元中取定幾個點作為節(jié)點，例如對三角形單元一般取其三個頂點作為節(jié)點。然后對于單元內(nèi)的被求函數(shù)的局部變化特征作出假設(shè)，也就是選定型線或插值函數(shù)，

5、一般選用多項式作為插值函數(shù)。例如對三角形單元（圖15.1），可設(shè)：第十五章 流體動力學(xué)問題的有限元解法 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第3頁頁式中，系數(shù) ， 及 可以用三個節(jié)點 的值來表示，符號 表示被求函數(shù) 的近似表達式。顯然這一多項式不可能恰好是式（15.1）真正的解，即若把它代入到該式中，則其右側(cè)不會等于零。將不等于零的部分作為余量，記為r，利用加權(quán)余量法，要求余量在某種意義上為最小，即要求：式中，a為區(qū)域的面積，w為權(quán)函數(shù)。上式要求余量r與權(quán)函數(shù)在a區(qū)域上的內(nèi)積為零。由于近似解 是用末知節(jié)點上的函數(shù)值來表示的，因而式（15

6、.2）給出了這些未知節(jié)點上函數(shù)值之間的代數(shù)關(guān)系式，即有限元法的離散方程式。求解這些方程，可得到有限元法的數(shù)值解。在有限元法中導(dǎo)出離散方程的方法較多，有變分法，最小二乘方法，加權(quán)余量法等，其中加權(quán)余量法應(yīng)用范圍較廣，本書中僅介紹這一種方法。 1a2a3a 0ddaaarwawfl（15.2）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第4頁頁二、有限元法區(qū)域的離散化在二維問題中最常用的單元為三角形單元，其三個頂點可作為節(jié)點，如圖15.2(b)所示，是二維問題的線性單元。它對不規(guī)則區(qū)域的適應(yīng)性較好（參見圖15.1）。圖15.2 一維和二維線性單元(

7、b)ixjy(a)xkh12312對于一維問題，單元都是直線段，每個單元上的節(jié)點數(shù)取決于所選定的型線。若選用線性函數(shù)作單元上的函數(shù)逼近，則在該單元上只需兩個節(jié)點，線性方程中的兩個未知量（截距與斜率）可由這兩個節(jié)點上的未知函數(shù)值確定，具有這種特性的單元稱為線性元，如圖15.2(a)所示。如果在直線段的單元上選取多于兩個以上的節(jié)點，就為非線性單元。本章中僅介紹線性單元。第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第5頁頁在二維問題中最常用的單元為三角形單元，其三個頂點可作為節(jié)點，如圖15.2(b)所示，是二維問題的線性單元。它對不規(guī)則區(qū)域的適應(yīng)性較

8、好（參見圖15.1）。圖15.3 稀疏矩陣及帶寬將一個求解區(qū)域劃分成許多相連接又不重疊的子區(qū)域的過程就是區(qū)域離散化，子區(qū)域就是單元。一般先將子區(qū)域分成四邊形與三角形的組合，然后再將四邊形細(xì)分成三角形。要注意不能把一個三角形的頂點取在任一相鄰三角形一條邊的中間位置上。單元的尺寸及疏密程度據(jù)物理問題的性質(zhì)及對計算精度的要求而定。一般而言，物理量變化劇烈的地方單元的尺寸要小一些、排列要密集一些。有限元法的計算精度受到單元內(nèi)最長邊與最短邊長度比的影響，盡量不要把三角形劃分成鈍角三角形，因為那樣會使長短邊之比增加而使計算精度下降。第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本

9、思想與區(qū)域離散化 第第6頁頁區(qū)域離散化過程中需給單元與節(jié)點編號。整個計算區(qū)域中的單元是統(tǒng)一編號的，單元號的字符用e表示，其值從1開始，順序增加。與有限差分法不同，有限元法中的節(jié)點有兩個編號，即單元節(jié)點號（局部的）和總體節(jié)點號。在單元中節(jié)點一般用i，j和k（或1，2，3）按逆時針方向編號（圖15.2(b)），整個計算區(qū)域內(nèi)的節(jié)點則按一定的順序統(tǒng)一編號?？傮w節(jié)點編號的原則是盡可能使同一單元內(nèi)各節(jié)點的編號相近，因為單元節(jié)點號的差值，決定了所形成的代數(shù)方程系數(shù)矩陣的特性。有限元法所生成的代數(shù)方程的系數(shù)矩陣是一個稀疏矩陣，即系數(shù)矩陣中有相當(dāng)多的元素為零。如果總體節(jié)點編號合適，同一單元中各節(jié)點的編號相差較

10、小，可以使非零元素相對集中地分布在系數(shù)矩陣的對角線附近。如圖15.3所示，從對角線到非零元素所在區(qū)邊界之間的距離稱為帶寬。在用直接解法求解代數(shù)方程組時（有限元法所生成的代數(shù)方程組多用直接解法求解），在采用一定的處理方法后，只需把非零元素輸入計算機進行計算。在一定的總節(jié)點數(shù)下，帶寬越窄，需送入計算機的非零元素越少，所占用的計算機內(nèi)存就越小。 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第7頁頁(a)(b)圖15.4 節(jié)點整體編號可以證明，如果計算區(qū)域中每個單元三個節(jié)點編號數(shù)的最大值為r，則帶寬與 成正比。圖15.4所示為同一計算區(qū)域節(jié)點的兩種編號

11、方式，方式(a)可比方式(b)節(jié)省一半以上的計算機內(nèi)存。1r第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第8頁頁三、單元的插值函數(shù)對于一維問題的線性單元（圖15.2(a)），設(shè)插值函數(shù)為由1、2兩節(jié)點上的函數(shù)值 、 可得 的表達式為：式中為書寫方便略去了表示近似值的符號。將該式代入式（15.3）得xaa211221, aahahxxaijijji21,jiijijijjihxxhxxxhhxx（15.3）（15.4）（15.5a）這里h是線性單元的長度。顯然這一插值函數(shù)對計算區(qū)域中的各個單元都是適用的。取其中任一單元為e，則上式可寫成為：jej

12、ieinn)()( ( ),ejeiijxxxxnnhh（15.5b）（15.5c）其中：和 稱為單元e的形狀函數(shù)。 ( )ein ejn第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第9頁頁對于三角形線性元（圖15.2(b)），若討論導(dǎo)熱問題，可假設(shè)單元中的溫度為x，y的線性函數(shù)，則有：其中待定常數(shù) 、 及 可由節(jié)點上的溫度值來表示。為此將三個節(jié)點i，j和k的坐標(biāo)代入上式（區(qū)域離散化后，各個節(jié)點的位置坐標(biāo)及單元的面積均為已知），將所得的三元一次代數(shù)方程組寫成矩陣形式，有：yaxaa3211a2a3akjikkjjiitttaaayxyxyx32

13、1111（15.6）（15.7）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 為書寫方便，令：第第10頁頁利用矩陣求逆的方法，可得a1，a2及a3：11231111111iiijjjkkkijkkjkiikijjijkkjijjiikjjkjikjjkkxytaaxytaxyttx yy xx yy xx yx yyyyyyytxyxxxxxxtxyxy,ijkkjijkikjjkiikjkijikkijjikijkjiax yx ybyy cxxax yx ybyy cxxax yx ybyy cxx第十五章 流體動力學(xué)問題的有限元解法 退 出返

14、 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第11頁頁再將分母中的行列式展開：可證明 之值等于該三角形單元面積的兩倍。把三角形單元的面積記為 ，則有：ijjikikikjkjkkkjkjkikikkjjiicbcbxxyyxxyyyxyyxxyyxxyxyxyx100111ijjicbcbkjikjikjikjitttcccbbbaaaaaa21321（15.8）將a1、a2及a3的表達式代入（15.6），可得單元函數(shù) 的插值計算式：kekjejieitntntn)()()(（15.9a）其中 、 為單元形狀函數(shù)，其計算式為：( )ein( )( )eejknn及)(21)(21)(21)

15、()()(ycxbanycxbanycxbankkkekjjjejiiiei（15.9b）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第12頁頁在獲得了單元形狀函數(shù)后就可以構(gòu)造總體的形狀函數(shù)，即構(gòu)造整個求解區(qū)域上被求函數(shù)的一種近似形式。我們以一維問題為例來說明。如圖15.5(a)所示，該計算區(qū)域有三個單元，圖中在括號內(nèi)的數(shù)字表示單元編號，水平線以上的數(shù)字表示單元節(jié)點，水平線以下的數(shù)字表示總體節(jié)點號。注意到同一個節(jié)點的相鄰的兩個單元中的局部編號是不同的。例如總體編號為2的節(jié)點在單元中為編號2而在單元中則為1。在每個單元中兩個插值函數(shù)是線性函數(shù)，如

16、圖15.5(b)所示。把各個單元內(nèi)的插值函數(shù)疊加起來，有： 式中上角標(biāo)表示單元編號，下角標(biāo)為該單元中節(jié)點的編號?？紤]到：于是得：)3(2)3(2)3(1)3(1)2(2)2(2)2(1)2(1)1(2)1(2)1(1)1(1nnnnnn4)3(23)3(1)2(22)2(1)1(21)1(1,443322114)3(23)3(1)2(22)2(1)1(21)1(1nnnnnnnnnn（15.10）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第一節(jié) 有限元法的基本思想與區(qū)域離散化 第第13頁頁nnnn(e)(d)(f )(g)1234111112134(a)(b)(c)1圖15.5 整體

17、形狀函數(shù)1121234122式中n1，n2，n3及n4就是相對于節(jié)點1，2，3及4的整體形狀函數(shù)，而上式也就是在整個求解區(qū)域中函數(shù) 的近似表達式。注意每個整體形狀函數(shù)在相應(yīng)的節(jié)點上都取得“1”的值，而在其余節(jié)點上為零，這使上述的整體函數(shù)近似表達式能夠滿足在不同節(jié)點上取得該節(jié)點函數(shù)值的要求。在圖15.5(c)中畫出了各個單元中的局部形狀函數(shù)，它們的疊加所形成四個整體形狀函數(shù)如圖15.5(d)15.5(g)所示。 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第二節(jié) 有限元法中代數(shù)方程的建立第第1頁頁本節(jié)主要介紹有限元法中離散方程建立的原理與過程，包括galerkin（伽遼金）加權(quán)余量法的原

18、理，單元矩陣的生成，總體矩陣合成及邊界條件的處理等內(nèi)容。一、加權(quán)余量法加權(quán)余量法是獲得微分方程近似解的一種有效方法。假設(shè)微分方程（15.1）中的未知函數(shù) 可以近似地表示成為：miiicn1（15.11）式中， 是一些所選定的線性獨立的函數(shù)（即其中任意一個函數(shù)都不能由其它函數(shù)經(jīng)過線性運算而得出），ci為未知的變量，m為未知變量的個數(shù)。加權(quán)余量法要求將（15.11）代入微分方程后所得到的余量在整個計算區(qū)域上與所選定的權(quán)函數(shù) 的內(nèi)積為零，即滿足：inmiwi, 1 miawrawfliaia, 2, 1, 0dd（15.12）也就是在某種平均意義上要求余量為零。 第十五章 流體動力學(xué)問題的有限元解法

19、 退 出返 回 第第2頁頁第二節(jié) 有限元法中代數(shù)方程的建立加權(quán)函數(shù)的不同選擇導(dǎo)致多種不同的加權(quán)余量法，其中應(yīng)用較廣的galerkin余量法選定權(quán)函數(shù) 。因而galerkin余量法要求：即要求余量（誤差）的加權(quán)平均值在整個計算區(qū)域上應(yīng)等于零。這樣m個積分式就產(chǎn)生m個代數(shù)方程，從而可以解出m個未知量 。由于式（15.11）中的形狀函數(shù)都是坐標(biāo)的線性函數(shù)，而導(dǎo)熱問題控制方程的最高階導(dǎo)數(shù)為二階，線性函數(shù)的二階導(dǎo)數(shù)為零，因此不能直接將式（15.11）代入上式進行積分計算。為了克服這一困難可對式（15.13）作分部積分，對于一維問題，在區(qū)域 上有：,1, 2,iiwnim mianflia, 2, 1,

20、0dicbxabababauvuvvudd（15.13）（15.14a）式中，u相當(dāng)于式（15.13）中的權(quán)函數(shù) （即wi），而dv則相當(dāng)于微分算子l。分部積分的結(jié)果就可以把包括在算子符號內(nèi)的函數(shù)的導(dǎo)數(shù)降低一階。對于二維及三維的問題，分部積分相當(dāng)于應(yīng)用gauss降維定理。設(shè)有一空間區(qū)域，其體積為，表面積為a，則gauss降維定理為：in第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第3頁頁第二節(jié) 有限元法中代數(shù)方程的建立式中，u相當(dāng)于權(quán)函數(shù)； 相當(dāng)于微分算子，為矢量運算 ；n為邊界外法線上的單位矢量。ddduauduavnvvvdkjizyx（15.14b）12a1a2a3443a圖1

21、5.6 分段線性近似解分部積分引起了兩點變化：（1）對函數(shù)的近似表達式的要求降低了，只要求一階導(dǎo)數(shù)存在即可；（2）引入了邊界條件，式（15.14a）中積分的上下限就是引入了邊界條件。由于式（15.14a）右端第二項的積分 在計算時是按單元分段進行的，因而只需在每個單元上近似解的一階導(dǎo)數(shù)存在且在單元邊界上函數(shù)值連續(xù) 即 可 。 例 如 對 一 維 問 題 的 近 似 解 式（15.10），可用分段線性函數(shù)近似表示（圖15.6）bavdu第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第4頁頁第二節(jié) 有限元法中代數(shù)方程的建立在得出了分部積分的表達式后，可將所假定的近似解（如式（15.11）所

22、示）代入進行積分計算，對每一個整體形狀函數(shù)都要在整個計算區(qū)域內(nèi)作積分，每完成這樣一個積分就得出一個代數(shù)方程式。但是注意到每個整體形狀函數(shù)實際上只在一個單元或相似的幾個單元之內(nèi)值不為零，在其它單元上其值均為零（如圖15.5(d)(g)所示），且在整體形狀函數(shù)不等于零的單元內(nèi)，整體形狀函數(shù)之值就等于該單元的形狀函數(shù)之值。因而整體形狀函數(shù)在整個求解區(qū)域內(nèi)的積分形成代數(shù)方程的過程為：對每個單元按單元形狀函數(shù)作積分，然后把共享一個節(jié)點的各單元的積分結(jié)果按一定方式相加。由于每個單元的單元形狀函數(shù)都是一樣的，因而只要對一個代表性單元作積分即可。對每個單元積分所形成的該單元內(nèi)的代數(shù)方程的矩陣稱為單元矩陣，其實

23、施過程稱為單元分析。把共享一個節(jié)點的各單元的結(jié)果疊加以形成總體代數(shù)方程的過程叫作總體合成。單元分析與總體合成是有限元法中建立離散方程的兩個重要步驟。 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第5頁頁第二節(jié) 有限元法中代數(shù)方程的建立二、單元矩陣的生成下面我們以一維擴散問題為例來說明單元分析與總體合成的過程。設(shè)式（15.1）中的微分算子為一二階導(dǎo)數(shù)，即 ，任一單元e的形狀函數(shù)為 ，則對該單元內(nèi)的積分有：即對每個單元，余量的加權(quán)平均值也等于零。式中x1，x2為單元e的兩個端點。對該式做分部積分，可得： xfxfl22dd)2, 1( inei 2, 1, 0ddd21)(22ixxnx

24、fxxxei 2, 1, 0dddddddd212121ixxnxfxxnxxnxxeixxeixxei（15.15a）（15.15b）利用 的表達式（15.5b），有： eeieeeeiieixnxnxnxnxnxnxdddddddddddddd2121221121（15.15c）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第6頁頁第二節(jié) 有限元法中代數(shù)方程的建立式中，符號 表示行矢量， 表示列矢量。將式（15.15c）代入式（15.15b），得：對上式右端的第一項計算如下： xxnxfxnxxnxnxxeixxeiexxeieidddddddd212121)(dd)(dd)()(

25、dd)(dd, 2)(dd)(dd)()(dd)(dd, 1211)(222)(2)(111)(122)(1)(2121xxxxxnxxxnxnixxxxxnxxxnxnieexxeieexxei（15.15d）這里已把圖15.5(b)中單元形狀函數(shù)的兩個端點值代入。若采用下列符號： xxnxffxxnxffxxnxnkxxnxnkxxnxnkxxnxnkexxeexxexxeeexxeeexxeeexxeeedddddddddddddddddddddd22112222122121121111212121212121（15.15e）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第7頁頁

26、第二節(jié) 有限元法中代數(shù)方程的建立則式（15.15d）可寫成：該式即單元e的代數(shù)方程，矩陣 及 稱為單元矩陣。在獲得了單元代數(shù)方程后就應(yīng)將共享一個節(jié)點的各個單元的矩陣元素按一定規(guī)則疊加，以形成整個計算區(qū)域的代數(shù)方程。在疊加過程中，每個單元代數(shù)方程（式（15.15f）右端第一項的疊加結(jié)果，除了整個計算區(qū)兩端單元外，其余都將相互抵消為零。而這涉及到該問題的邊界條件。 )(21)(21)(11)(22211211ddddeeeeffxxxxkkkk ek ef（15.15f）三、總體矩陣的合成與邊界條件的處理關(guān)于邊界條件引入到總體代數(shù)方程中去的方法在介紹了總體合成規(guī)則以后再討論。為了便于敘述，以具有三

27、個單元的一維擴散問題為例說明?？傮w合成的步驟如下：第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第8頁頁第二節(jié) 有限元法中代數(shù)方程的建立（1）列出各單元節(jié)點的局部編號與總體編號間的對照表。對如圖15.5a所示的三個單元的區(qū)域，這一對應(yīng)關(guān)系如表15.1所示； 單元節(jié)點號局部總體112122122331234（2）將每個單元方程中矩陣元素的下標(biāo)由局部節(jié)點號換成為總體節(jié)點號；（3）將每個單元的矩陣元素按其經(jīng)過轉(zhuǎn)換的下標(biāo)號插入在到相應(yīng)的總體矩陣的相應(yīng)位置上。本例中有四個節(jié)點，因而矩陣k就為44方陣，而有關(guān)邊界條件及源項的兩列矢量均有四個分量；表15.1 局部節(jié)點編號與總體節(jié)點編號的對照第十五章

28、 流體動力學(xué)問題的有限元解法 退 出返 回 第第9頁頁第二節(jié) 有限元法中代數(shù)方程的建立（4）將位于總體矩陣同一位置上的各個單元矩陣的元素相加，形成總體矩陣的元素。這一疊加過程表明總體矩陣中的元素代表了共享一個節(jié)點的各單元所作出的貢獻之和。 3433232212114321344343334333233232223222122121112111dd000dd000000fffffflxxkkkkkkkkkkkk對所舉的例子，總體合成所得結(jié)果如下列矩陣方程所示：式中l(wèi)為計算區(qū)域的長度。為書寫簡便起見，把近似函數(shù) 直接寫為 。等號右端的兩個列矢量都是已知值，是代數(shù)方程組中的常數(shù)項，可合并記為列矢量r

29、。于是此總體代數(shù)方程組便可表示為下列矩陣形式：第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第10頁頁第二節(jié) 有限元法中代數(shù)方程的建立 式中，n為計算區(qū)域中的總節(jié)點數(shù)；矩陣k可稱為熱導(dǎo)矩陣，相當(dāng)于有限差分法中系數(shù) 等所形成的矩陣，反映了熱傳導(dǎo)的機理對溫度場的影響；列矢量r稱為熱載荷矢量，代表了熱邊界條件、內(nèi)熱源項等的影響。以下討論如何把邊界條件應(yīng)用于矩陣方程（式（15.16）中。對第一類邊界條件，可以采用有限差分法中對某點賦給定值的方法。例如設(shè)區(qū)域的左端點為給定值 ，則可以用一個很大的數(shù)，如1020，乘k中的第一個對角元，然后將r中的第一個元素以 來代替，結(jié)果形成以下矩陣方程： 11

30、nn nnkreagk201110（15.16） 2020111111213142212223242331323334344142434441010gkkkkkkkkkrkkkkrkkkkr第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第11頁頁第二節(jié) 有限元法中代數(shù)方程的建立顯然求解可得 。對于給定熱密度的邊界條件，只要把給定的 之值代入r即可。對第三類邊界條件，若設(shè)區(qū)域的右端點為第三類邊界條件，則有：式中對流換熱系數(shù) 、周圍流體溫度 均為已知?？紤]到導(dǎo)出總體矩陣方程時所依據(jù)的控制方程式 相當(dāng)于導(dǎo)熱系數(shù)為1的情形，為協(xié)調(diào)起見，這里也取 （ 時分析的方法完全一樣），因而有 。顯然只要把

31、r第四個元素中的 用 來代替，并把 中的 用 來代替即可把對流邊界條件考慮在內(nèi)。g1xddflxttx4ddft0)(dd22xfx11 4ddttlxf)(ddlxft k44k44k至此，可得采用galerkin加權(quán)余量法導(dǎo)出有限元法離散方程的步驟：1寫出所求解問題的控制方程及其邊界條件的表達式；2將galerkin加權(quán)余量法的表達式作一次分部積分，得出降維后的積分表達式；第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第6頁頁第二節(jié) 有限元法中代數(shù)方程的建立3區(qū)域離散化，剖分單元，確定節(jié)點及其編號；4選定單元插值函數(shù)，導(dǎo)出單元形狀函數(shù)；5單元分析，將所選定的近似函數(shù)代入降維后的積分

32、表達式，導(dǎo)出單元矩陣；6總體合成，將單元矩陣的元素按一定的方式疊加形成總體有限元方程；7處理邊界條件。至于所形成的代數(shù)方程的求解與有限差分法中的代數(shù)方程求解一樣，也有直接解法與迭代法兩大類。在有限元法中，只要節(jié)點的總體編號合適，可以把總體代數(shù)方程組系數(shù)矩陣中的非零元素集中在對角線附近區(qū)域，同時采用壓縮存儲的技術(shù)，可以使總體熱導(dǎo)矩陣中需要存儲的元素個數(shù)遠(yuǎn)遠(yuǎn)小于該矩陣的元素總數(shù)，因而多采用直接解法。有關(guān)壓縮存儲技術(shù)的內(nèi)容可參見相關(guān)文獻。 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第1頁頁第三節(jié) 二維邊值問題有限元法求解舉例 例例15.1設(shè)有兩塊形狀為等邊三角形的平行平板，相距為h，其間

33、充滿了一種粘性不可壓縮的流體，粘度為 。該兩平板以速度w向前運動（圖15.7(a)），同時粘性流體則從板中被擠出（軸承潤滑是這一問題的工程背景）。由流體力學(xué)可知，此時兩平板中壓力場的分布由poisson方程描述：(a)xzhywxy=132(b)14(x,y)32y=2-3x3ooxy圖15.7 例15.2的圖示為了使讀者進一步熟悉有限元法的求解步驟，本節(jié)中以一個二維邊值問題為例給出有限元法求解的全過程。為便于手工計算，單元數(shù)取得很少，主要目的是使讀者進一步熟悉單元矩陣的生成及總體矩陣的合成過程。第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第2頁頁第三節(jié) 二維邊值問題有限元法求解舉例

34、 同時，在平板的邊界上流體壓力為零。三角形的邊長為 。試用有限元法確定三角形區(qū)域中任一點的流體壓力值。解解：為能進行手工計算在流場（平面流場，沿高度方向壓力是均勻的）中只取四個節(jié)點，即三角形區(qū)域的三個頂點和一個內(nèi)點4，其坐標(biāo)為 。并且為使計算結(jié)果能應(yīng)用于內(nèi)部各點，這里x，y暫不取定具體數(shù)值。節(jié)點的總體編號及單元編號如圖15.7(b)所示。單元矩陣的元素為： const1232222ghwypxp3/2yx, kjijkjiiyxynynxnxnkediiiieji,dd,及 kjiiyxgnrediei,ddkjii,式中 表示下標(biāo)i作ijk的輪換。計算形函數(shù)的導(dǎo)數(shù)項，并代入上式，得： kji

35、jkjiiccbbyxccbbkjijiedjijieejie,41dd412,及第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 對單元3有：第第3頁頁第三節(jié) 二維邊值問題有限元法求解舉例 式中 為單元面積。將三個單元矩陣中的節(jié)點號均轉(zhuǎn)換為總體節(jié)點號，如表15.2所示。 3ddddedidieigyxngyxgnree e于是，對單元1有：對單元2有： (1)(1)(1)1112141212224241424333kkkpgkkkpkkkp(2)(2)(2)2223242323334342434443kkkpgkkkpkkkp(3)(3)(3)3331343131114144341443k

36、kkpgkkkppkkk第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第4頁頁第三節(jié) 二維邊值問題有限元法求解舉例 將上述三個單元方程的元素裝入到總體方程矩陣中：在上式中為表明元素來自哪一個單元，在其右上角用小括號標(biāo)明了單元號。注意式中每個元素均已為整體坐標(biāo)。(1)(3)(1)(3)(1)(3)1111112131414(1)(1)(2)(2)(1)(2)2212222232424(3)(2)(2)(3)(2)(3)3313233333434(1)(3)(1)(2)(2)(3)(1)(2)(3)4414142424343444444pkkkkkkpkkkkkkpkkkkkkpkkkkk

37、kkkk(1)(3)(1)(2)(2)(3)(1)(2)(3)3g 下面來確定邊界條件。由于上式中1、2、3節(jié)點處的流體壓力 ， 和 均為已知（等于零），因而采用消行修正法后只得出一個關(guān)于 的方程：或1p2p3p4p(1)(2)(3)(1)(2)(3)44444443gpkkk (1)(2)(3)4(1)(2)(3)4444443gpkkk 第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第5頁頁第三節(jié) 二維邊值問題有限元法求解舉例 (1)44k 244k 344k(1)222244441221(1)(1)2(1)(1)11()()() 44121()433kbcyyxx(2)22444

38、4(2)(2)11()43kbc(3)224444(3)(3)11()43kbc其中， 、 和 的值計算如下：代入p4的計算式得：而 的關(guān)系可由三角形面積的計算公式求得：(1)(2)(3)(1)(2)(3)444(2)(3)(1)(3)(1)(2)()()(,)p xyg ),(44)3()2()1(yx與及、)3(321)32(3213213244)3(44)2(44)1(yxxyyy第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第6頁頁第三節(jié) 二維邊值問題有限元法求解舉例 代入p4計算式并整理有：可知在給定了 后即可求得p4。有限元法是50年代從固體力學(xué)中發(fā)展起來的，60年代逐漸推

39、廣到傳熱學(xué)與流體力學(xué)的計算中。雖然有限元法在區(qū)域的離散（常稱剖分）代數(shù)方程的形成等方面工作量比較大，但由于它對不規(guī)則區(qū)域的適應(yīng)性強而得到重視，目前已是計算流體力學(xué)中應(yīng)用較廣的一種數(shù)值方法，并且隨著計算機技術(shù)的高速發(fā)展，將在計算流體動力學(xué)中發(fā)揮更大的作用。 )( 3)3(2 3)3)(32(),(24244444444444yxxyyxyxygyxp44,yx第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第1頁頁第四節(jié) 有限分析法介紹 有限分析法是80年代初發(fā)展起來的一種數(shù)值計算方法。在計算流場方面它可以克服在高re數(shù)下的有限差分法容易振蕩或發(fā)散的缺點。本節(jié)中以二維穩(wěn)態(tài)導(dǎo)熱問題為例來介紹

40、它的基本思想。用有限分析法求解邊值問題的主要步驟如下：（1）把求解區(qū)域用一系列與區(qū)域邊界平行的網(wǎng)格線進行離散，兩根網(wǎng)格線的交點為計算節(jié)點，每一個節(jié)點與其相鄰的四個網(wǎng)格組成一個單元，即每一個單元由一個內(nèi)點和八個鄰點組成（圖15.8），圖中有陰影線的部分即為p單元。應(yīng)指出，這樣定義的單元在每個相鄰單元間有一部分區(qū)域是相互重疊的。（2）在單元內(nèi)將控制方程的非線性項局部線性化，并對單元邊界上未知函數(shù)的變化型線作出選擇，把所選定型線表達式中的常數(shù)或系數(shù)項用單元邊界節(jié)點上的值來表示。（3）找出上述條件下單元內(nèi)未知函數(shù)的分析解（相當(dāng)于邊界上未知函數(shù)為給定的第一類邊界條件）。 第十五章 流體動力學(xué)問題的有限元

41、解法 退 出返 回 第第2頁頁第四節(jié) 有限分析法介紹 （4）選用所假定的分析解找出其內(nèi)節(jié)點及八個鄰點的未知函數(shù)之間的關(guān)系式。即建立該內(nèi)點的離散方程。（5）逐一地對求解區(qū)域內(nèi)每一節(jié)點建立起離散方程，對非第一類邊界條件的邊界節(jié)點，利用與該節(jié)點相鄰單元中的分析解建立一個補充方程。（6）求解所建立的代數(shù)方程組，在獲得了求解區(qū)域每一個點上的函數(shù)值后再進行其它量（如邊界上法向?qū)?shù)）的計算。一個單元上九個節(jié)點的命名方式如圖15.9所示。下面研究在這一單元上如何建立相應(yīng)于laplace方程 =0的離散方程。 2首先確定單元邊界上用以逼近被求變量 的近似函數(shù)的形狀。一般有三種選擇。以右邊界為例來寫出：（1）分段

42、線性型線，以ec為原點分為上、下兩支，即：)0()0()(22yhycahyybayeeeee（15.17a）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第3頁頁第四節(jié) 有限分析法介紹 （2）二項式分布：（3）多項式指數(shù)函數(shù)：其中，系數(shù) 可采用邊界上三點的未知值 、 表示，而系數(shù)a則在求解過程中確定，上述函數(shù)都是單調(diào)變化的。如果在求解前能預(yù)期到相鄰兩邊界點間被求函數(shù)不是作單調(diào)變化時，應(yīng)在該兩節(jié)點間設(shè)置更多節(jié)點，以使每兩鄰節(jié)點間被求函數(shù)呈單調(diào)變化。 )()(222hyhycybayeeee2( )ayeeeeyab yc eeeecba、neseec和（15.17b）（15.17c）對

43、于分段線性及二項式分布這兩種型線采用多項式擬合法來導(dǎo)出系數(shù) 、 與邊界節(jié)點函數(shù)值間的關(guān)系。把坐標(biāo)原點設(shè)在ec上，則對分段線性的型線有：由此得：eaebec2222)(,)(,)0(hcahhbahaeeeseeeeneeeec22,hchbaseeceecneeece（15.18a）第十五章 流體動力學(xué)問題的有限元解法 退 出返 回 第第4頁頁第四節(jié) 有限分析法介紹 對二項式分布有：由此得：可計算當(dāng) 、 、 和 均已知時，這一單元上 的分析解。這是第一類邊界條件下矩形域上的laplace方程，利用分離變量法可求得其分析解為：22222222)()(,)0(hchbahhchbahaeeeeseeeeeneeeecseecneeseneeecehchba221,21,222 ye yw yn ys02222yx（15.18b）1112112122211sinchsh1sinchsh,hxhydhydhyhxchxcyxnnnnnnnnnnnn（15.19）其中： 22d1sin2sh122111hhnwennyhyyfyfhhhc（15.20