有限體積法介紹

有限體積法1有限體積法基本原理上一章講到的有限差分法將數(shù)值網(wǎng)格的節(jié)點(diǎn)上定義為計(jì)算節(jié)點(diǎn)，并在網(wǎng)格節(jié)點(diǎn)上對(duì)微分形式的流體基本方程進(jìn)行離散，用網(wǎng)格節(jié)點(diǎn)上的物理量的代數(shù)方程作為原PDE的近似。在本章所要學(xué)習(xí)的有限體積法則采用了不同的離散形式。首先，有限體積法離散的是積分形式的流體力學(xué)基本方程：(1)Jp°v.nds=\「n.V°ds+J(1)計(jì)算域用數(shù)值網(wǎng)格劃分成若干小控制體。和有限差分法不同的是，有限體積法的網(wǎng)格定義了控制體的邊界，而不是計(jì)算節(jié)點(diǎn)。有限體積法的計(jì)算節(jié)點(diǎn)定義在小控制體內(nèi)部。一般有限體積法的計(jì)算節(jié)點(diǎn)有兩種定義方法，一種是將網(wǎng)格節(jié)點(diǎn)定義在控制體的中心，另一種方法中，相鄰兩個(gè)控制體的計(jì)算節(jié)點(diǎn)到公共邊界的距離相等。第一種方法的優(yōu)點(diǎn)在于用計(jì)算節(jié)點(diǎn)的值作為控制體上物理量的平均值具有二階的精度；第二種方法的好處是在控制體邊界上的中心差分格式具有較高的精度。ooooGoooOQoQoooooooQoooOoQQo0Ooooo積分形式的守恒方程在小控制體和計(jì)算域上都是成立的。為了獲得每一個(gè)控制體上的代數(shù)方程，面積分和體積分需要用求面積公式來(lái)近似。2面積分的近似采用結(jié)構(gòu)化網(wǎng)格，在二維情況下，每一個(gè)控制體有4個(gè)面，二維情況，每一個(gè)控制體有6個(gè)表面。計(jì)算節(jié)點(diǎn)用大寫(xiě)字母表示，控制體邊界和節(jié)點(diǎn)用小寫(xiě)字母表示。為了保證守恒性，控制體不能重疊，每一個(gè)面都是相鄰兩個(gè)控制體的唯一公共邊界?？刂企w邊界上的積分等于控制體個(gè)表面的積分的和：(2)JfdS=ZJfds(2)S kSk出上式中，f可以表示p°un或「總o0NNooGNWnwN。nneNE00wWW。w swp,ec seE?EE。sWSE”o◎oSS°0c顯然，為了獲得邊界上的積分，必須知道f在邊界上的詳細(xì)分布情況，這是不可能實(shí)現(xiàn)的，由于只是計(jì)算節(jié)點(diǎn)上的函數(shù)值，因此必須采用近似的方法來(lái)計(jì)算積分。整個(gè)近似過(guò)程分成兩步第一步：用邊界上幾個(gè)點(diǎn)的近似積分公式第二步：邊界點(diǎn)上的函數(shù)值用計(jì)算節(jié)點(diǎn)函數(shù)值的插值函數(shù)近似面積分可采用以下不同精度的積分公式：二階精度積分：F=』fds=fS氏fS (3)eS eeeee上式中fe為邊界中點(diǎn)出的函數(shù)值。近似為方格中心點(diǎn)的值乘以方格的面積。三階精度積分：(4)四階精度積分：F=Jfds六'J4L+fS (5)eS 6 ee應(yīng)該注意的是，采用不同精度的積分公式，在相應(yīng)的邊界點(diǎn)的插值時(shí)也應(yīng)采用相應(yīng)精度的插值函數(shù)。積分公式的精度越高，近似公式就越復(fù)雜。3體積分的近似和面積分相似，體積分也有不同精度的近似公式階精度積分公式(6)Q=Jqds=qS氏qAQ(6)Se采用雙二次樣條函數(shù)+ax+ay+ax2+ay2+axy+ax2y+axy2+ax2y2123 4 5 6 7 8可以得到四階精度的積分公式：AQ。Q=JAQ。Q=Jqdsxu6s36SeqP+4q+4q+4q+4q+4q +4q+4q+4qs wn s se swne nw(8)4函數(shù)的插值在上節(jié)講到的積分的近似公式中用到了非計(jì)算節(jié)點(diǎn)上的函數(shù)值，被積函數(shù)f中包含了多個(gè)物理量及其偏微分，如對(duì)流項(xiàng)fc=pGv-n，擴(kuò)散項(xiàng)fd=rn-V。，在源項(xiàng)中也有類(lèi)似情況，這里假定流場(chǎng)和流體的物性參數(shù)是已知的，物理量。及其偏導(dǎo)數(shù)在控制面上的值需要通過(guò)計(jì)算節(jié)點(diǎn)上物理量的插值得到。下面已e面為例進(jìn)行討論。迎風(fēng)插值（UDS）巾e用上游計(jì)算節(jié)點(diǎn)的函數(shù)值近似相當(dāng)于對(duì)一階偏導(dǎo)數(shù)采用迎風(fēng)格式，因此用UDS來(lái)表示這種近似方法，在UDS中:(9)if(v.n)>0(9)if(v-n)<0eUDS是唯一無(wú)條件滿足有界性要求的近似格式，在數(shù)值過(guò)程中不會(huì)產(chǎn)生數(shù)值振蕩。UDS存在數(shù)值粘性。根據(jù)Taylor公式，該格式具有一階精度，并具有數(shù)值粘性：rnum=（pU）Ax/2 （10）在多維問(wèn)題中，如果流動(dòng)方向和網(wǎng)格是斜交的，截?cái)嗾`差會(huì)在垂直于流動(dòng)方向以及流線方向產(chǎn)生擴(kuò)散，這是一種非常嚴(yán)重的誤差，函數(shù)的峰值或函數(shù)值的快速變化會(huì)被抹平，為了得到高精度結(jié)果需要采用非常精細(xì)的網(wǎng)格。線性插值（CDS）(11)(12)線性插值具有二階精度，線性插值相當(dāng)于FDM中的CDS格式，因此用CDS表示。CDS格式會(huì)產(chǎn)生數(shù)值振蕩。對(duì)于擴(kuò)散項(xiàng)(13)

三階迎風(fēng)格式（QUICK）和UDS類(lèi)似，QUICK格式也和流動(dòng)方向有關(guān)g1aE-g20W+(I-g1+gg1aE-g20W+(I-g1+g2)°P

g。一g0+(1-g+g)。v3P4EE 3 4Eif(v.n)>0if(v-n)<0e其中：g—G-九e,W2P.gg1 1+九「九W’ 2l+九l-九)g= e,W e,Pg3 1+九-九e,E e,P(1-九％-九) e,P e,W-1+入e，「入eW九2九g= e:Ee:P41+匕e'p(14)(15a)(15b)4.4高階格式（4階精度CDS）采用三次曲線可擬合出四階精度的中心插值公式，在均勻網(wǎng)格中，四階公式為:270P+270「30「30”.P E48""W E(16)「。0) 270-270+0-0 = E P W EE、Sx) 24Axe(17)5邊界的處理對(duì)于對(duì)流項(xiàng)，在入口處一般給出了流量或函數(shù)值，在邊界和對(duì)稱(chēng)面上流量為零，在出口處假設(shè)和出口的法向坐標(biāo)無(wú)關(guān)，因此可采用迎風(fēng)格式。對(duì)于擴(kuò)散項(xiàng)則可能需要采用偏心格式。6有限體積法應(yīng)用舉例例：考慮一標(biāo)量在已知流場(chǎng)中的輸運(yùn)過(guò)程（如圖4.4所示），輸運(yùn)方程為:Jp0v.ndS=Jrv0-ndSS S邊界條件：(18)0=0；北部入口邊界0=1-y；西部壁面邊界對(duì)稱(chēng)條件；南部邊界梯度為0；東部出口條件對(duì)流項(xiàng):對(duì)流項(xiàng):em=Jpv-ndS=(pu)△y為質(zhì)量通量。Imax(m,0)巾+min(rn,0)巾 forUDS(20)(21)Fc=4 eP eE(20)(21)e Im(1—九)巾+m九。 forCDSIeePeEE若采用UDS格式，代數(shù)方程組中各項(xiàng)系數(shù)為:Ac=min(m,0)；Ac=min(m,0)E e W wAc=min(m,0)；Ac=min(m,0)N n S sAc=—(Ac+Ac+Ac+Ac)P EWNS若采用CDS格式，代數(shù)方程組中各項(xiàng)系數(shù)為Ac=m九；Ac=m'EeeWwwTOC\o"1-5"\h\zAc=m九；Ac=m九 (22)N nn S ss\o"CurrentDocument"Ac=—(Ac+Ac+Ac+Ac)P E W NS根據(jù)連續(xù)性方程：m+m+m+m=0 (23)ewns相鄰CV之間的關(guān)系：m =—m；九=1—X (24)其余相鄰CV有類(lèi)似關(guān)系擴(kuò)散項(xiàng)采用CDS格式

(25)Fc=』rv^-ndSe(25)esSe代數(shù)方程組中擴(kuò)散項(xiàng)系數(shù)為:- rAyAdEAd二— AdEWx-xAd=NrAxAdWrAx(26)Ad=—(Ad+Ad+AdAd=NrAxAdWrAx(26)Ad=—(Ad+Ad+Ad+Ad)

P EWNS對(duì)于任意控制體A^+A巾+A巾+A巾+A巾=Q(27)A,=Ac+Ad,1為任意指標(biāo)P,E,W,S,No (28)邊界條件的處理：對(duì)于西部和北部邊界，由于給定了函數(shù)值，對(duì)流項(xiàng)可直接代入函數(shù)值而無(wú)需插值，擴(kuò)散項(xiàng)則采用一側(cè)差分(29)這里，W點(diǎn)和P的w邊中點(diǎn)重合。南邊和西邊的梯度為零,以南邊為例,由于梯度為零，6p=6s,代數(shù)方程變?yōu)?A6+(A+A)巾+A6+A6=Q(30)6SIMPLE方法考慮定常不可壓流動(dòng)問(wèn)題，控制方程為:連續(xù)性方程：Jpv-ndS=0SV動(dòng)量方程：Jpvv-ndS=JNn-VvdS—JpndS+JpbdQ(31)(32)SV SV SV CV不可壓縮問(wèn)題求解的困難在于壓力場(chǎng)的求解。主要原因在于壓力p沒(méi)有獨(dú)立的方程組。先考慮一維問(wèn)題：對(duì)于動(dòng)量方程：(puu)—(puu)Su(su1N—I—P+P\SxJ eww(33)若采用CDS格式(puu)+(p若采用CDS格式(puu)+(puu) (puu)+(puu)p——2u-uu-u=u—E P—u—P WAx Ax簡(jiǎn)化后得：(puu)-(puu)E——2wAx-2u ppE+pw22W(34)根據(jù)連續(xù)性方程，u「丫ui-i=根據(jù)連續(xù)性方程，u「丫ui-i=c系方程，容易造成壓力交錯(cuò)現(xiàn)象。為了解決這一問(wèn)題，可采用交錯(cuò)網(wǎng)格技術(shù)，即速度場(chǎng)和壓力場(chǎng)采用不同的網(wǎng)格。主控制體為壓力控制體（黑色實(shí)線網(wǎng)格），u的控制體（紅色虛線網(wǎng)格）的計(jì)算節(jié)點(diǎn)在主控制體的e邊，控制體的e，w邊界通過(guò)主控制體的計(jì)算節(jié)點(diǎn)，v控制體（藍(lán)色雙點(diǎn)劃線網(wǎng)格）的計(jì)算節(jié)點(diǎn)在主控制體的n邊，該控制體的n，s面經(jīng)過(guò)主控制體的計(jì)算節(jié)點(diǎn)。在u的控制體中，采用有限體積法離散可得u的代數(shù)方程：au=Eanvn=Zau=Eanvn=Zau+Q+(p-p)AnbnbpEeav+Q+(p-p)A(35)nbnbpNe壓力場(chǎng)的求解采用壓力校正方法。即采用預(yù)估的壓力場(chǎng)求速度，再用連續(xù)性方程校正壓力場(chǎng)。當(dāng)連續(xù)性方程得到滿足時(shí)，壓力場(chǎng)就是真實(shí)的壓力場(chǎng)。具體步驟如下1．2．1．2．3．預(yù)測(cè)壓力場(chǎng)p將預(yù)測(cè)壓力場(chǎng)代入動(dòng)量方程，分別求解速度場(chǎng)u,Vau =Za u +Q+ （p -p ）Aee nbnb pEeav =2a v +Q+ （p -p ）A用連續(xù)性方程校正壓力(36)設(shè)方程的精確解為u，v，p(37)u=u+u'；v=v+v'；p=p+p'(37)其中u',v',p'為校正量。則校正量滿足方程:au'ee__Javau'+Q+(p'

nbnb Pav，+Q+(p'nbnb P—p))AEe—p))ANe(38)略去相鄰節(jié)點(diǎn)速度校正量的影響，可得:A()()U=e\p—p7=d\p—p)

eaPEePEe(39a)M=土(p，—p,)=d(p，—p,)

naPNnPNn代入連續(xù)性方程：(39b)(U+u)A—(U+u)A+(v+vr)A—(v+vr)A=0 (40)eeewwwnnnsss整理得：dA (p' —p')—dA (p' —p' )+d A (p'—p')—dA (p' —p')e,P e,P PE e,W e,W WP n,P e,