有限體積法的基本理論與應(yīng)用

有限體積法的基本理論與應(yīng)用

0fvm的發(fā)展有限體積法是20世紀(jì)80年代以來(lái)開(kāi)發(fā)的一種新的微分方程色散方法，具有獨(dú)特的優(yōu)勢(shì)。現(xiàn)在它是計(jì)算位移問(wèn)題、數(shù)值計(jì)算和數(shù)值模擬的重要方法之一。長(zhǎng)期以來(lái),微分方程數(shù)值求解的離散方法中,使用得最多的是有限差分法(FDM)和有限元方法(FEM).有限差分法著眼于求解區(qū)域剖分的節(jié)點(diǎn)上的函數(shù)值,方法簡(jiǎn)便、靈活,離散的格式豐富多樣,在收斂性、穩(wěn)定性等理論研究方面也比較完善.但由于計(jì)算中對(duì)求值節(jié)點(diǎn)的分布要求比較規(guī)則,不能適應(yīng)復(fù)雜的幾何求解域.另外,它不能求許多實(shí)際物理問(wèn)題中出現(xiàn)的弱解.有限元方法基于微分方程的弱解形式和廣義變分原理,利用能適應(yīng)復(fù)雜的幾何形狀的求解區(qū)域網(wǎng)格剖分,在剖分單元上用形函數(shù)插值逼近來(lái)求解,但有時(shí)需要求解大型的線性方程組,計(jì)算上沒(méi)有差分方法那樣靈活方便,并且處理大變形間斷問(wèn)題較困難.FVM在一定程度上吸收了FDM與FEM的長(zhǎng)處,克服了它們的缺點(diǎn).FVM從控制體的積分形式出發(fā),對(duì)求解區(qū)域的剖分同F(xiàn)EM一樣具有單元特征,能適應(yīng)復(fù)雜的求解區(qū)域,離散方法具有差分方法的靈活性,間斷解的適應(yīng)性.FVM的原始雛形從60年代開(kāi)始出現(xiàn),如Harlow等人提出的PIC、MAC和FLIC方法,70年代的MacCormark和Patankar等人的FVM思想.80年代以來(lái),網(wǎng)格生成技術(shù)、特別是無(wú)結(jié)構(gòu)網(wǎng)格生成技術(shù)的發(fā)展,給FVM注入了強(qiáng)大的生命力,使FVM進(jìn)入了一個(gè)快速發(fā)展時(shí)期.同時(shí),FVM結(jié)合其他一些數(shù)值方法,如有限元方法,Upwind(迎風(fēng))格式,Runge-Kutta方法,還有新發(fā)展起來(lái)的MUSCL、Roe、TVD、ENO等方法,產(chǎn)生了一系列新穎有效的方法.1維堅(jiān)守型問(wèn)題的控制方程守恒型問(wèn)題的FVM是典型的FVM.我們以二維守恒型問(wèn)題為例,介紹典型的FVM的離散過(guò)程、基本格式,給出FVM的一個(gè)基本輪廓.二維守恒型問(wèn)題的控制方程的形式為ut+(f(u))x+(g(u))y=0(1)這里u、f、g是標(biāo)量或向量,u是待求變量.二維Euler方程和守恒型淺水波方程都是這種形式.1.1無(wú)結(jié)構(gòu)網(wǎng)格的應(yīng)用有限體積法首先要對(duì)求解區(qū)域進(jìn)行網(wǎng)格剖分,網(wǎng)格分為有結(jié)構(gòu)網(wǎng)格和無(wú)結(jié)構(gòu)網(wǎng)格兩種.有結(jié)構(gòu)網(wǎng)格的生成技術(shù)比較成熟,網(wǎng)格的生成和安排是有序的和按一定結(jié)構(gòu)的,適合于求解區(qū)域有幾何規(guī)律和求解模型比較均衡的問(wèn)題.無(wú)結(jié)構(gòu)網(wǎng)格是當(dāng)前數(shù)值計(jì)算中一個(gè)受人關(guān)注的重要發(fā)展方向,其網(wǎng)格的生成和安排可以是無(wú)序的,可根據(jù)需要加密或減疏,易于構(gòu)造自適應(yīng)網(wǎng)格,這是有結(jié)構(gòu)網(wǎng)格難以做到的.無(wú)結(jié)構(gòu)網(wǎng)格能適應(yīng)復(fù)雜的求解區(qū)域,這一點(diǎn)也比有結(jié)構(gòu)網(wǎng)格優(yōu)越.原則上講,無(wú)結(jié)構(gòu)網(wǎng)格可以具有任意結(jié)構(gòu),但為了實(shí)現(xiàn)和計(jì)算上的方便,多采用同一類(lèi)型的網(wǎng)格,對(duì)二維問(wèn)題多為三角形,三維問(wèn)題多為四面體.我們這里將以二維問(wèn)題的三角剖分來(lái)說(shuō)明問(wèn)題.1.2維三角剖分網(wǎng)絡(luò)結(jié)構(gòu)作網(wǎng)格剖分之后,就要根據(jù)問(wèn)題的特點(diǎn)和需要來(lái)確定控制元.然后,在控制元上積分原方程并進(jìn)行離散和數(shù)值計(jì)算.控制元的類(lèi)型有兩種,一種是將單一的網(wǎng)格單元(這里我們選取的是三角元)作為控制元,另一種是將共角點(diǎn)的網(wǎng)格各取一部分合在一起作為控制元.二維三角剖分的共角點(diǎn)的網(wǎng)格組合的控制元又分為:(1)統(tǒng)一型(圖1(a));(2)垂心型(圖1(b));(3)完全中心型(圖1(c));(4)部分中心型(圖1(d)).在二維情形,控制元可以統(tǒng)一描述為k(k≥3)個(gè)頂點(diǎn)的多邊形.1.3u3000正當(dāng)化.記控制元為V,且設(shè)V是k個(gè)頂點(diǎn)的多邊形,在V上對(duì)方程(1)積分得∫V?u?tdxdy+∫V(?f?x+?g?y)dxdy=0.(2)一種簡(jiǎn)單的考慮是將?u/?t在V上看作常數(shù),由Green公式可得Ac?u?t+∮?Vfdy-gdx=0?(3a)或Ac?u?t+∮?VF?ndl=0?(3b)其中Ac為控制元V的有向面積,F=(f,g),n=(nx,ny).我們首先介紹根據(jù)式(3a)得出的一階正則格式或單調(diào)格式.在式(3a)中,沿V的邊線用梯形公式近似計(jì)算積分,有∮?Vfdy-gdx=k∑j=1[12(fi+fj-1)(yj-yj-1)-12(gj+gj-1)(xj-xj-1)].記δxj=12(xj+1-xj-1),δyj=12(yj+1-yj-1),并約定x0=xk,xk+1=x1,y0=yk,yk+1=y1,則∑δxj=0,∑δyj=0.這時(shí)∮?Vfdy-gdx=k∑j=1(fjδyj-gjδxj).(4)將式(4)代入式(3a),并引入f、g在控制元共角點(diǎn)或中心點(diǎn)o的值fo、go.因?yàn)椤苂oδxj=0,∑foδyj=0,可得?u?t=-1Ack∑j=1[(fj-fo)δyj-(gj-go)δxj]?(5)若定義aj={-(fj-fo)δyj-(gj-go)δxjuj-uo?uj≠u(mài)o,-(?f?uδyj-?g?uδxj)|u=uo,uj=uo,(6)則(5)式成為?u?t=1Ack∑j=1aj(uj-uo).(7)對(duì)時(shí)間導(dǎo)數(shù)項(xiàng)采用向前差分離散可得un+1o=uno+ΔtAck∑j=1aj(unj-uno)=(1-ΔtAck∑j=1aj)uno+ΔtAck∑j=1ajunj.(8)記c0=1-ΔtAc∑kj=1aj,cj=ΔtAcaj,即得un+1o=c0uno+k∑j=1cjunj.(9)若cj>0(j=0,…,k),則式(7)為正則格式.正則格式必為單調(diào)格式.若式(9)不是正則格式,可采用補(bǔ)償形式,un+1o=c0uno+k∑j=1cjunj+k∑j=0αjunj=(c0+α0)uno+k∑j=1(cj+αj)unj?(10)使式(10)成為正則格式.2fvm數(shù)值解中間斷解的應(yīng)用非線性的雙曲型守恒方程(組),其解可能存在間斷,這給方程的求解帶來(lái)非常大的困難.20世紀(jì)中期以來(lái),間斷解問(wèn)題的研究一直是偏微分方程數(shù)值解中最引人注目的領(lǐng)域之一.近年來(lái),FVM在離散方案中引入當(dāng)前極為活躍的一些新方法,產(chǎn)生了非常豐富的求解間斷解問(wèn)題的新算法.2.1數(shù)值流fj、gj的計(jì)算假定方程(1)是非線性雙曲型守恒方程(組),則方程的求解必須考慮到間斷解的處理.現(xiàn)假定Vi是k條邊圍成的控制元,lij是Vi的第j條邊,處理間斷解的一種途徑是從(3a)式出發(fā),將(3a)式中的線積分化為∫?Vifdy-gdx=k∑j=1[fj?(yj-yj-1)-gj?(xj-xj-1)]?(11)其中數(shù)值流fj、gj是f、g在邊lij上的近似值.在這種離散方式下,問(wèn)題的關(guān)鍵就在于如何準(zhǔn)確地計(jì)算數(shù)值流fj、gj.下面以fj為例介紹(gj同理).在間斷解情況下,由Godunov思想,lij兩側(cè)的函數(shù)值可以認(rèn)為是不同的,即有跳躍,因而fj的計(jì)算必須考慮這種間斷,即fj=f*(uL,uR)?其中f*是數(shù)值流函數(shù),uL,uR分別為u在邊lij逆時(shí)針?lè)较蜃蟆⒂覀?cè)的近似值.數(shù)值流fj的計(jì)算有兩個(gè)重要的方面,一是數(shù)值流函數(shù)f*的形式,二是uL,uR的計(jì)算(通常稱(chēng)為重構(gòu)(reconstruction)).2.1.1值流函數(shù)的結(jié)構(gòu)中小型f全卷或f#全卷/fur/ur型f*(uL,uR)=12(f(uR)+f(uL))?或f*(uL,uR)=f(uR+uL2)?(12)這是一種最簡(jiǎn)單的格式,效果相對(duì)較差.r-um-3流函數(shù)f*(uL,uR)=12(f(uR)+f(uL)-α(uR-uL))?(13)其中α=maxu|?f/?u|,最大值在u的相關(guān)范圍中取.另外一些單調(diào)流函數(shù),如Godunov流函數(shù),Engquist-Osher流函數(shù)等,詳見(jiàn)文.roe的線性化替代矩陣f*(uL,uR)=12[f(uR)+f(uL)-|?a|(uR-uL)]?(14)其中?a是a=?f/?u在uL,uR的Roe平均上的值.當(dāng)方程(1)是方程組時(shí),則?a是Jacobian矩陣a=?f/?u的Roe平均線性化替代矩陣.這時(shí),如果設(shè)?a的特征值和特征向量為:?λk??ek?k=1?2???m?將未知函數(shù)跳躍量表示成特征向量的線性組合:uR-uL=m∑k=1?αk?ek?(15)(15)式代入(14)式,就得到Roe的特征向量展開(kāi)方法的公式為f*(uL,uR)=12[f(uR)+f(uL)-m∑k=1?αk|?λk|?ek]?(16)這一形式在實(shí)際計(jì)算中廣泛使用,效果很好.sl、sr的計(jì)算f*(uL,uR)=1SR-SL[SRf(uL)-SLf(uR)+SRSL(uR-uL)].(17)這里SL、SR是左、右傳波的波速,SL、SR的一種近似計(jì)算方法見(jiàn)文.這一形式的計(jì)算比③中的Roe形式稍簡(jiǎn)單一點(diǎn),在規(guī)則網(wǎng)格上能取得較好的效果,但在無(wú)結(jié)構(gòu)網(wǎng)格上的效果明顯不及Roe形式.2.1.2中心型fvmuL,uR的重構(gòu)是FVM中重要的環(huán)節(jié),決定方法的空間精度、分辨率.為介紹方法,將配置在控制元Vi的中心的守恒量uo也記為ui,且不妨設(shè)與Vi有公共邊lij的控制元為Vi+1,ui+1為Vi+1的中心的守恒量.uL,uR的重構(gòu)有下面這些不同的方法:①分片常數(shù)逼近:uL=ui,uR=ui+1.(18)這是一階精度方法,穩(wěn)定性較好,但精度不高.②分片線性逼近的MUSCL格式:由于線性逼近要引進(jìn)參考量ui-1、ui+2,故這種方法以控制元為單一網(wǎng)格時(shí)計(jì)算較為簡(jiǎn)便.如果控制元為單一三角元,且守恒量配置在控制元的中心點(diǎn),記ui-1、ui+2分別是Vi、Vi+1的與邊lij相對(duì)的角點(diǎn)上的守恒量,如圖2所示.角點(diǎn)上的守恒量的值可以利用相鄰中心點(diǎn)的值作距離倒數(shù)加權(quán)平均來(lái)求出.這時(shí){uR=ui+1-12φ(ri+1)?δui+12?uL=ui+12φ(ri)?δui-12.(19)這里δui+12=ui+1-ui?δui-12=ui-ui-1?ri=δui+12/δui-12?ri+1=δui+32/δui+12?φ是限制器(Limiter)函數(shù).這是二階精度方法,計(jì)算量中等,計(jì)算結(jié)果較好.③Upwind(迎風(fēng))型格式:u(x,y)=uo+Δuo?r?(20)這里的r為點(diǎn)o到點(diǎn)(x,y)的向量,Δuo為u在o點(diǎn)的梯度,其近似值可以通過(guò)公式[ΗS2]Δuo=1A∮?Aundl(21)來(lái)計(jì)算,其中積分路徑?A是點(diǎn)o周?chē)阎鋟值的點(diǎn)的連線所組成,A是?A所圍區(qū)域的有向面積.分別在Vi和Vi+1上使用公式(20)、(21),就可以得出uL、uR.為確保重構(gòu)的單調(diào)性,實(shí)際計(jì)算中需要在公式(20)中加入限制器φ,即[ΗS1*3/4]u(x,y)=uo+φ(r)Δuo?r?(22)這也是二階精度方法,而且是標(biāo)準(zhǔn)的中心型FVM,計(jì)算結(jié)果比MUSCL格式還好,但計(jì)算量較大.④ENO或WENO重構(gòu):ENO、WENO重構(gòu)采用避開(kāi)間斷的思想,在均勻矩形網(wǎng)格上可以構(gòu)造出高階逼近的uL、uR,同時(shí)式(3a)中線積分的離散采用相應(yīng)的高階格式,能得到很好的計(jì)算結(jié)果,但計(jì)算代價(jià)非常大.在非均勻網(wǎng)格上的計(jì)算則更加困難.采用不同形式的流函數(shù)和不同的uL、uR重構(gòu)方法,就可以得到各種各樣的FVM,例如,流函數(shù)采用Roe的Riemann解算子逼近形式,uL、uR的重構(gòu)采用MUSCL方法,就得到Roe型MUSCL方法的FVM.2.2復(fù)合形式的fvm求解非線性雙曲型守恒方程(組)(1)的FVM的另一途徑是從式(3b)式出發(fā),將式(3b)式中的線積分化為∮?ViF?ndl≈k∑j=1(F?n)ijlij?(23)這里也用lij記邊lij的長(zhǎng)度,數(shù)值流(F·n)ij是外法向流向量函數(shù)F·n在邊lij上的近似值.數(shù)值流的計(jì)算是這一類(lèi)方法的設(shè)計(jì)和構(gòu)造的關(guān)鍵所在,構(gòu)造方法與上面介紹的fj、gj的構(gòu)造方法相同,只是要采用復(fù)合的形式.例如,數(shù)值流函數(shù)的Roe的Riemann解算子逼近形式為(F?n)ij=12[(F(uR)+F(uL))?nij-|A?|(uR-uL)]?(24)其中A?是矩陣A=?(F·n)/?u的Roe的線性化替代矩陣在邊lij上的值.復(fù)合形式的各種構(gòu)造方法在水動(dòng)力學(xué)的二維淺水波方程的求解中使用得較多,取得了較好的效果.我們也就二維潰壩問(wèn)題對(duì)FVM的復(fù)合形式的幾種方法進(jìn)行了數(shù)值實(shí)驗(yàn).下面給出實(shí)驗(yàn)的結(jié)果.二維潰壩問(wèn)題:考慮一個(gè)面積為200m×200m的區(qū)域,正中有一條大壩將水庫(kù)和下游一分為二,上游水深10m,下游水深5m.大壩距離水庫(kù)一邊95m處,突然有一段75m長(zhǎng)的壩體倒塌,要求模擬這部分壩體倒塌后水流的情況.問(wèn)題的平面圖見(jiàn)圖3.控制方程是二維淺水波方程:Ut+fx+gy=0?U=[h?uh?vh]Τ?f=[uh?u2h+gh2/2?uvh]Τ?g=[vh?uvh?v2h+gh2/2]Τ.實(shí)驗(yàn)結(jié)果見(jiàn)圖4、圖5.幾種方法的精度、計(jì)算速度和穩(wěn)定性的分析比較,我們將另文給出.2.3種各樣限制器函數(shù)的比較為自適應(yīng)地調(diào)節(jié)方法的數(shù)值耗散與數(shù)值色散效應(yīng),保持格式的單調(diào)性,達(dá)到數(shù)值結(jié)果的高分辨率和高保真的目的,通常需要在上述二階類(lèi)方法,如MUSCL、Upwind等方法中使用限制器(Limiter).有各種各樣的限制器函數(shù)見(jiàn)諸于文獻(xiàn)中,其中最有名的限制器有:Minmod:φ=max[0,min(r,1)],Superbee:φ=max[0,min(2r,1),min(r,2)],vanLeer:φ=(r+|r|)/(1+r).理論分析和數(shù)值實(shí)驗(yàn)的結(jié)論是,在上面三種限制器中,Minmod數(shù)值耗散性最強(qiáng),vanLeer的次之,Superbee數(shù)值耗散性最弱.而Superbee數(shù)值色散性最強(qiáng),Minmod數(shù)值色散性最弱.因此,通常Superbee給出一個(gè)較陡的剖面圖,Minmod產(chǎn)生一個(gè)相對(duì)平坦的剖面圖,vanLeer的居中.在實(shí)際計(jì)算時(shí),我們可根據(jù)不同的需要來(lái)選用它們.3維淡水波方程的數(shù)值方法及實(shí)際應(yīng)用80年代以來(lái),FVM在實(shí)際應(yīng)用中取得了很大的成功,在解決實(shí)際問(wèn)題的同時(shí),也發(fā)展了FVM的數(shù)值方法.這方面的進(jìn)展大致可以分為以下幾個(gè)方面.①Euler方程的求解:JamesonA等人在有結(jié)構(gòu)網(wǎng)格上,結(jié)合Runge-Kutta方法,用FVM對(duì)二維Euler方程進(jìn)行了數(shù)值求解.WangJCT等人對(duì)Euler方程的求解構(gòu)造了TVD型FVM,并進(jìn)行了爆炸流場(chǎng)的數(shù)值模擬.PanD等人對(duì)Euler方程構(gòu)造了無(wú)結(jié)構(gòu)網(wǎng)格上的二階迎風(fēng)型FVM.LevequeRJ給出了利用波傳播(Wavepropagation)的高精度FVM.這方面其它一些工作見(jiàn)文.②淺水波問(wèn)題的求解:由AlcrudeF等人給出求解二維淺水波方程的Roe型MUSCL方法的FVM,并進(jìn)行了水動(dòng)力學(xué)中三個(gè)經(jīng)典的間斷解問(wèn)題的數(shù)值模擬,呈現(xiàn)的結(jié)果令人鼓舞.AnastasiouK等人給出了求解二維淺水波方程的無(wú)結(jié)構(gòu)三角網(wǎng)格上的Roe型迎風(fēng)方法的FVM,并成功地進(jìn)行了包括二維潰壩問(wèn)題在內(nèi)的五種水力學(xué)數(shù)值實(shí)驗(yàn).HuK等人給出了二維淺水波方程求解的HLL型MUSCL方法的FVM,成功地進(jìn)行了四種開(kāi)管流(Open-channelflows)的數(shù)值實(shí)驗(yàn).這方面另外一些工作見(jiàn)文.③對(duì)流—擴(kuò)散問(wèn)題的求解:PatankarSV建立了傳熱和流體流動(dòng)研究中的FVM.FeistauerM等人研究了非線性對(duì)流—擴(kuò)散問(wèn)題的結(jié)合有限體積和有限元的數(shù)值方法.CorderoE等人研究了對(duì)流—擴(kuò)散方程求解的一種新的FVM.④湍流問(wèn)題的研究:HalloL等人給出了三維湍流問(wèn)題的一種混合有限體積和有限元的隱式數(shù)值方法,成功地進(jìn)行了低Mach數(shù)噴射湍流(turbulentlow-Mach-numberjet),超音速混合層流(supersonicmixinglayer)和三維管流(three-dimensionalpipe)等數(shù)值模擬.文中也用FVM研究了湍流模型.⑤一般雙曲問(wèn)題的求解:SonarT研究了一般雙曲守恒律問(wèn)題的ENO型FVM,系統(tǒng)地研究了適合于無(wú)結(jié)構(gòu)網(wǎng)格上FVM的ENO重構(gòu)的基函數(shù)(recoveryfunctions).BallandR等人給出了變系數(shù)雙曲問(wèn)題求解的單元—頂點(diǎn)型FVM.⑥一般Navier-Stokes(N-S)方程的求解:PanD等人給出了求解N-S方程的無(wú)結(jié)構(gòu)三角網(wǎng)格上的迎風(fēng)型FVM.DespotisGK等人給出了求解不可壓N-S方程的無(wú)結(jié)構(gòu)三角網(wǎng)格上的分?jǐn)?shù)步FVM.曾揚(yáng)兵等人給出了N-S方程在無(wú)結(jié)構(gòu)網(wǎng)格上的一種FVM求解方法.另外,在橢圓型問(wèn)題,兩相流問(wèn)題以及其它一些問(wèn)題的FVM數(shù)值方法和實(shí)際應(yīng)用的進(jìn)展見(jiàn)文[40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58].4fvm的收斂性FVM理論分析方面的工作有一定難度,特別是在無(wú)結(jié)構(gòu)網(wǎng)格FVM方面尤其困難.盡管如此,近年來(lái),這方面工作