有限元法在流體力學(xué)中的應(yīng)用

1、第五章 有限元法在流體力學(xué)中的應(yīng)用 本章介紹有限元法在求解理想流體在粘性流體運(yùn)動(dòng)中的應(yīng)用。討論了繞圓柱體、翼型和軸對(duì)稱物體的勢(shì)流，分析了求解粘性流動(dòng)的流函數(shù)渦度法流函數(shù)法和速度壓力法，同時(shí)導(dǎo)出粘性不可壓流體的虛功原理。§1 不可壓無(wú)粘流動(dòng)真實(shí)流體是有粘性和可壓縮的，理想不可壓流體模型使數(shù)學(xué)問(wèn)題簡(jiǎn)化，又能較好地反映許多流動(dòng)現(xiàn)象。1. 圓柱繞流本節(jié)詳細(xì)討論有限無(wú)法的解題步驟?？紤]兩平板間的圓柱繞流如圖51所示。為了減小計(jì)算工作量，根據(jù)流動(dòng)的對(duì)稱性可取左上方的l4流動(dòng)區(qū)域作為計(jì)算區(qū)域。 選用流函數(shù)方法，則流函數(shù)應(yīng)滿足以下Laplace方程和邊界條件 (5-1)將計(jì)算區(qū)域劃分成10個(gè)三角形單

2、元。單元序號(hào)、總體結(jié)點(diǎn)號(hào)和局部結(jié)點(diǎn)號(hào)都按規(guī)律編排如圖52所示。從剖分圖上所表示的總體結(jié)點(diǎn)號(hào)與單元結(jié)點(diǎn)號(hào)的關(guān)系，可以建立聯(lián)綴表于下元素序號(hào)12345678910總體結(jié)點(diǎn) 號(hào)n11444226655n2459865710109n322593637810 表5-1各結(jié)點(diǎn)的坐標(biāo)值可在圖52上讀出。如果要輸入計(jì)算機(jī)運(yùn)算必須列表。本質(zhì)邊界結(jié)點(diǎn)號(hào)與該點(diǎn)的流函數(shù)值列于下表邊界結(jié)點(diǎn)號(hào)n12348910流函數(shù)2221000 表5-2選用平面線性三角形元素，插值函數(shù)為(315)式。對(duì)二維Laplace方程進(jìn)行元素分析，得到了單元系數(shù)矩陣計(jì)算公式(319)和輸入向量計(jì)算公式(320)?，F(xiàn)在對(duì)全部元素逐個(gè)計(jì)算系數(shù)矩陣。

3、例如元素1，其結(jié)點(diǎn)坐標(biāo)為=0, =2; =0, =1; =2.5, =2.由(315)式可得; ,; ; 從(319)式可計(jì)算出依次可計(jì)算出全部子矩陣根據(jù)聯(lián)綴表把元素矩陣組合成總體系數(shù)矩陣A=矩陣中零元素沒(méi)有一一寫出，下三角部分與上三角部分對(duì)稱。從(320)式計(jì)算元素輸入向量，由于流函數(shù)滿足齊次的自然邊界條件,所以輸入向量為零，總體輸入向量也為零，這樣就得了總體有限元方程.式中: 用縮減方程的重新編號(hào)修正方法施加邊界條件，本質(zhì)邊界結(jié)點(diǎn)的函數(shù)值是已知的。把它們代入方程，修正右端項(xiàng)，再減去相應(yīng)的方程，整理得 解方程得到 0.845，=1.241，=1.121 這樣求出了全部結(jié)點(diǎn)上的流函數(shù)。為了求出

4、每個(gè)單元形心處的速度，可以由單元的流函數(shù)近似表達(dá)式求導(dǎo)計(jì)算。對(duì)元素e來(lái)說(shuō)，有 例如單元=2, =1, =3，這樣計(jì)算得到的速度為u=1,0。二維繞圓柱流動(dòng)還可以用勢(shì)函數(shù)求解，則定解問(wèn)題可寫成表示勢(shì)函數(shù)，為了使數(shù)值解唯一必須在部分邊界上給定本質(zhì)邊界條件。勢(shì)函數(shù)邊界同樣標(biāo)記在圖5l上。勢(shì)因數(shù)滿足Laplace方程和相應(yīng)的邊界條件，與流函數(shù)不同僅在于有非齊次的自然邊界條件。采用與流函數(shù)方法完全一樣的網(wǎng)格劃分，可知計(jì)算得到的單元系數(shù)矩陣是完全一樣的，總體矩陣也是完全一樣的。元素1和4具有非齊次自然邊界條件應(yīng)該用(320)式計(jì)算輸入向量。元素l, 。元素4, ?？傮w合成得到，這樣就得到方程組巳知，消去相

5、應(yīng)的三個(gè)方程得到一個(gè)7×7的代數(shù)方程組，解得單元形心處的速度可以用下列公式計(jì)算式中是單元的結(jié)點(diǎn)勢(shì)函數(shù)向量。對(duì)于元素1來(lái)說(shuō)， ,這樣計(jì)算得到u=1.033，v=-0.05。這結(jié)果與流函數(shù)方法得到的結(jié)果近似相等。如果加密網(wǎng)格，就可以得到更好的結(jié)果。2. 升力問(wèn)題 考慮圖53(a)所示的機(jī)翼繞流。均勻來(lái)流平行于x軸，機(jī)翼邊界為，后緣尖點(diǎn)為，流場(chǎng)外邊界取在離機(jī)翼足夠遠(yuǎn)處。流函數(shù) 滿足以下方程和邊界條件。 (5-3)其中a,b是特定系數(shù)，h是上下邊界之間的距離。機(jī)翼繞流的后駐點(diǎn)應(yīng)位于后緣尖點(diǎn)處，在后緣T點(diǎn)滿足Kutta條件 ； （5-4）由于方程和邊界條件是線性的,可用疊加原理求解，令 （5-

6、5） 其中，和：分別是下列問(wèn)題的解用有限元方法分別解以上三個(gè)問(wèn)題，得到各結(jié)點(diǎn)的、和，代入(55)式得到疊加解。顯然它滿足問(wèn)題(53)的全部方程和邊界條件，特定常數(shù)a,b可利用Kutta條件(54)定出。首先由流函數(shù)、和分別求出各個(gè)結(jié)點(diǎn)上的速度，和，然后在后緣點(diǎn)T處利用Kutta條件，應(yīng)有解之可得到a和b。 圖53(b)上給出了NACA4412具型以攻角置于均勻流場(chǎng)中所引起的流動(dòng)圖案，計(jì)算中采用了三角形單元。與無(wú)升力體繞流一樣，機(jī)具繞流也可以采用速度勢(shì)函數(shù)求解.3軸對(duì)稱問(wèn)題 考慮圓管內(nèi)繞軸對(duì)稱物體的無(wú)旋流動(dòng)，如圖54(a)所示。采用柱坐標(biāo)系(r，,z)，其勢(shì)函數(shù)滿足Laplace方程。 (5-6

7、) 寫出與微分問(wèn)題相應(yīng)的伽遼金積分表達(dá)=分部積分上式的左邊并整理得到弱解積分形式=式中L是元素的邊長(zhǎng)，L繞軸旋轉(zhuǎn)一周形成元素的邊界面。 采用圖54(b)所示軸對(duì)稱的環(huán)形線性元素，它是將平面線性三角形元素繞對(duì)稱軸旋轉(zhuǎn)一周形成的環(huán)形體。采用斜坐標(biāo)系，那么插值函數(shù)可寫成元素結(jié)點(diǎn)上勢(shì)函數(shù)向量為則逼近函數(shù)為總體坐標(biāo)和斜坐標(biāo)系的關(guān)系為式中。，是元素結(jié)點(diǎn)總體坐標(biāo)向量。 將逼近函數(shù)表達(dá)式代入伽遼金公式，推導(dǎo)出元素有限元方程式中影響系數(shù)矩陣和輸入向量分別為K=P= 求出插值函數(shù)向量的偏導(dǎo)數(shù)和，代入上式得影響系數(shù)矩陣K= （5-7）式中 ； i1，2，3時(shí)J=2，3，1;k31，2。 ， 三角形元素面積。假設(shè)元素

8、的“l(fā)一2”邊落在自然邊界上且q為常數(shù)，則可得轉(zhuǎn)入向量計(jì)算公式 （5-8）式中 ：是“l(fā)一2”邊的邊長(zhǎng)。計(jì)算了各元素的K和P，然后總體合成，代入本質(zhì)邊界條件就可以解總體方程。為了計(jì)算其它物理量，下面給出了相應(yīng)的公式。元素形心處的速度： （5-9）附加質(zhì)量m：m= (5-10)式中，i=1,2,.N. N是物體表面上所劃的單元數(shù)。p是輸入向量P在元素結(jié)點(diǎn)上的值。在文獻(xiàn)4中，以圓球作為例子計(jì)算了三種狀態(tài)。球在無(wú)阻空間中運(yùn)動(dòng)，計(jì)算的附加質(zhì)量系數(shù)0.4671，理論值是0.5。計(jì)算值小于理論值是符合第二章2節(jié)例4中證明的附加質(zhì)量極大值原理的.在圓管中球作勻速運(yùn)動(dòng)，計(jì)算的結(jié)果與TJChung在參考書(shū)(9)

9、中給出的結(jié)果比較，雖然我們采用了較少的結(jié)點(diǎn)，但達(dá)到了相同的精度。Chung用流函數(shù)方法，采用軸對(duì)稱四邊形單元計(jì)算。由于流函數(shù)滿足Stokes方程，是非自伴的，這樣行成的影響系數(shù)矩陣是非對(duì)稱的不僅計(jì)算麻煩而且不能利用半帶寬存儲(chǔ)。四邊形單元的短陣元素計(jì)算須用Guass數(shù)值積分,計(jì)算量大且有誤差。而采用勢(shì)函數(shù)方法和三角形元素恰好克服了以上兩個(gè)缺點(diǎn)。第三種狀態(tài)計(jì)算了圓球在半盲管(一端封死)中的運(yùn)動(dòng)。附加質(zhì)量系數(shù)0.897，這等于無(wú)限空間中附加質(zhì)量系數(shù)的1.6倍。這使我們想到，在計(jì)算水下管中發(fā)射彈道問(wèn)題時(shí)，應(yīng)考慮物體在管道中的附加質(zhì)量系數(shù)。軸對(duì)稱不可壓無(wú)粘流動(dòng)也存在看流函數(shù)提法，流函數(shù)應(yīng)滿足以Stoke

10、s方程，而不是Laplace方程。 （5-11） 應(yīng)特別注意的是,Stokes算子是非自伴的。為了寫出相應(yīng)的迦遼金弱解積分表達(dá)式,先可將方程改寫其伽遼金積分表達(dá)為將上式第一項(xiàng)分部積分，并代入自然邊界條件，得到=假設(shè)近似解可表示成那么 ， ，將以上各式代入弱積分表達(dá)式，得到單元有限元方程式中影響系數(shù)矩陣為 K式中第二項(xiàng)是非對(duì)稱項(xiàng)，它使得K成為非對(duì)稱矩陣。輸入向量為 P= 如果采用前面已用到的三角形洄旋環(huán)狀體元素，則是K和P可以導(dǎo)出，得到相應(yīng)的計(jì)算公式。看來(lái)流函數(shù)方法不及勢(shì)函數(shù)方法簡(jiǎn)便。§2 不可壓粘性流動(dòng) 不可壓粘性流體運(yùn)動(dòng)由速度散度為零的連續(xù)方程及NavierStokes方程描述。二

11、維問(wèn)題，引進(jìn)流函數(shù)可導(dǎo)出流函數(shù)渦量方程和四階的流函數(shù)方程。粘性流動(dòng)中存在粘性應(yīng)力，固壁邊界上必須滿足無(wú)滑條件，使得流體的運(yùn)動(dòng)一股是有旋的. 1流函數(shù)問(wèn)量法以流函數(shù)和渦量。表示的粘性不可壓流動(dòng)方程和自然邊界條件是 （5-12） （5-13）邊界上還應(yīng)給出本質(zhì)邊界條件，即和的函數(shù)值。根據(jù)具體邊界不難繪出流函數(shù)的邊界值。而固壁上的渦量，則要通過(guò)區(qū)域內(nèi)的流函數(shù)值，由渦量邊界條件來(lái)確定。流函數(shù)、渦量和速度的關(guān)系為分別寫出流函數(shù)方程和渦量方程的伽遼金積分表達(dá) 對(duì)以上兩式中Laplace算子進(jìn)行分部積分并代入自然邊界條件：得弱解積分形式=將求解區(qū)域剖分，單元中流函數(shù)、渦量。選用相同的插值展開(kāi)形式，則有 ，其

12、中 。是結(jié)點(diǎn)未知數(shù)向量。和分別是單元中i結(jié)點(diǎn)在t時(shí)刻流函數(shù)和渦量值。這樣就有，將流函數(shù)和渦量的插值表達(dá)式及其變分代入弱解積分形式，得=0整理得到單元的有限元方程 流函數(shù)的單元有限元方程是代數(shù)方程，而渦量的有限元方程是常微分方程組，必須在時(shí)間步進(jìn)中求解。式中各矩陣分別為質(zhì)量矩陣 損耗矩陣 對(duì)流矩陣 輸入向量為 是t的函數(shù)，對(duì)流矩陣A是與時(shí)間t有關(guān)的非對(duì)稱矩陣，因此每一時(shí)步必須重新計(jì)算。將所有單元進(jìn)行總體合成，得到總體有限元方程 （5-14） （5-15）式中和是總體結(jié)點(diǎn)未知數(shù)向量。 總體有限元方程是代數(shù)方程和非線性常微分方程的聯(lián)立方程組。將渦量力程用顯式差分貉式離散，可用交替迭代的方法求解耦合方

13、程組。 用流函數(shù)渦量法，可以計(jì)算低雷諾數(shù)時(shí)繞障礙，如圓柱、圓球的流動(dòng)，還可以模擬渦街的形成過(guò)程。 2流函數(shù)方法 二維不可壓粘性流動(dòng)可用四階的流函數(shù)方程描述 （5-16）寫出相應(yīng)的伽遼金積分表達(dá)利用分部積分=-可得弱解積分表達(dá)式 = 假定在全部邊界上給定本質(zhì)邊界條件，對(duì)于四階算子，即給定，那么在S上有及和等于零，則伽遼金弱積分形式可簡(jiǎn)化為設(shè)單元的近似函數(shù)為 式中將近似函數(shù)代入弱解表達(dá)式，可以推導(dǎo)出單元有限元方程。在推導(dǎo)中首先注意到全導(dǎo)數(shù)的運(yùn)算然后把近似解代入積分表達(dá)式中整理得到 （5-17）式中矩陣包含著未知向量 ，這是對(duì)流項(xiàng)的非線性影響產(chǎn)生的。如果寫成下標(biāo)表示的方程組，則可以看出方程的第二項(xiàng)包

14、含的二次項(xiàng)，所以得到的是非線性的常微分方程組.由于流函數(shù)方程是四階偏微分方程，因此函數(shù)及其一階導(dǎo)數(shù)、，是本質(zhì)變量，應(yīng)取為結(jié)點(diǎn)未知數(shù)并保持其在元素內(nèi)和元素之間的連續(xù)性。為此，應(yīng)選用函數(shù)及一階導(dǎo)數(shù)都連續(xù)的Hermite插值函數(shù)要構(gòu)造相容的Hermite插值函數(shù)有時(shí)是困難的。另外單元近似函數(shù)要用高階插值，元案自由度大，這使得總體有限元方程的系數(shù)矩陣階次很高要占用大量的計(jì)算機(jī)內(nèi)存.這種方法只具有一個(gè)未知函數(shù)，所以計(jì)算程序簡(jiǎn)單，也比較省時(shí)。用流函數(shù)方法計(jì)算圓柱繞流，得到的表面壓力分布曲線與有限差分解比較，符合得很好。§3. 流體力學(xué)虛功率原理及速度壓力法不可壓粘性流動(dòng)基本方程如下 連續(xù)方程 (

15、5-18)動(dòng)量方程 k=1,2 (5-19)本構(gòu)方程 (5-20)方程中重復(fù)指標(biāo)暗指疊加運(yùn)算。本構(gòu)方程給出了應(yīng)力和應(yīng)變率之間的關(guān)系。其中應(yīng)力張量是二階張量，表示作用在流體微團(tuán)小六面體的j面上k方向的應(yīng)力分量。由于流體微團(tuán)動(dòng)量矩平衡的原因它是對(duì)稱張量，即有是粘性應(yīng)力張量。 稱為Kronecker 是質(zhì)量力分量。流動(dòng)區(qū)域邊界可分成二類，一類是指定速度的，另一類是指定應(yīng)力。在指定速度的邊界上 k=1,2如圖55所示如果邊界的單位外法向矢量為n，其方向余弦j=1,2. 那么其法向速度應(yīng)滿足 (5-21)如固壁邊界，流體質(zhì)點(diǎn)粘附在固壁上，流體速度等于固壁運(yùn)動(dòng)速度，滿足無(wú)滑條件。若固壁靜止，則其值為零。在

16、指定應(yīng)力條件的邊界，上，有 是作用于邊界上，法向?yàn)閚的單位微元面積上的應(yīng)力矢量在k方向的分量見(jiàn)圖55。 設(shè)液體微元j面上的應(yīng)力矢量為 式中表示直角坐標(biāo)系中的兩個(gè)單位矢量。那么等于二個(gè)應(yīng)力矢量在n方向的投影之和，即 那么應(yīng)力邊界條件可寫成 （5-22）如自由液面上，是法向?yàn)閚的自由面上結(jié)定的作用于液體的風(fēng)應(yīng)力，多大氣靜止常壓時(shí)，其值為零.流體入口或出口邊界，可給出速度或應(yīng)力值。不同液體分界面上，一船給出速度，壓力連續(xù)性條件。把本構(gòu)方程(520)式代入動(dòng)量方程(519)式，并利用連續(xù)方程，就可以導(dǎo)出NavierStokes方程。下面推導(dǎo)虛功率原理。 以壓力的變分為權(quán)函數(shù)寫出連續(xù)方程的伽遼金積分公式

17、分部積分得到弱解積分公式 （5-23）以速度的變分為權(quán)函數(shù)，寫出動(dòng)量方程的伽遼金積分公式 （5-24）式中是重復(fù)指標(biāo),指疊加。伽遼金積分公式是兩個(gè)動(dòng)量方程分別乘以各自的權(quán)函數(shù)，作內(nèi)積后的和。式中應(yīng)力張量的導(dǎo)數(shù)包含著速度的二階導(dǎo)數(shù)項(xiàng)。對(duì)524)式中二所導(dǎo)數(shù)項(xiàng)進(jìn)行分部積分，并代入應(yīng)力邊界條件(522)，則=將上式代入(524)式=將(520)式代入上式，并考慮到全導(dǎo)數(shù)運(yùn)算，上式可以寫成= k=1,2 (5-25)這就是不可壓粘性流體的虛功率原理。方程右邊是由表面力和質(zhì)量力所作的外部虛功率，左邊是內(nèi)應(yīng)力的虛機(jī)械功率和慣性力的虛功率。虛功原理是壓分速度法有限元分析的基礎(chǔ)。設(shè)單元中壓力和速度采用相同的插

18、值展開(kāi)式 , , 結(jié)點(diǎn)未知數(shù)向量是時(shí)間t的函數(shù)。將近似解代入方程(623)注意到重復(fù)指標(biāo)暗指累加，得到 （5-26）式中 稱為連續(xù)矩陣，,為邊界流量向量。 將近似解代入(525)式為了簡(jiǎn)便寫成以下形式。 （5-27） （5-28）式中質(zhì)量系數(shù)矩陣 對(duì)流矩陣 耗散矩陣 壓力矩陣 外力向量 由于壓力和速度采用了相同的插值形式，所以連續(xù)矩陣與壓力矩陣相等。對(duì)流矩陣中包含著未知的速度向量，方程(528)，(527)和(528)構(gòu)成了非線性的常微分方程組。速度和壓力采取同一所次插值，速度可獲得較精確的解,但壓力會(huì)產(chǎn)生較大誤差。如果壓力采用線性插值而速度采用二次插值，則二者都可獲得好的結(jié)果。將單元有限元方程總體合成，代入速度邊界條件得到總體有限元