陳澤平,王晨濤,李 坤,馬天寶

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100081)

0 引言

爆炸與沖擊是在高溫高壓和相變等極端條件下,氣液固多介質(zhì)間強(qiáng)耦合作用的瞬態(tài)動(dòng)力學(xué)問題,數(shù)值求解該問題模型對(duì)航空航天等國防工業(yè)及武器裝備的研制開發(fā)均具有重要的基礎(chǔ)應(yīng)用價(jià)值[1-2]。計(jì)算爆炸力學(xué)很重要的一個(gè)問題在于精確捕捉爆轟波的波陣面及波陣面前后物理量的變化,但對(duì)于爆炸沖擊波的雙曲型守恒律方程,其解隨著時(shí)間的演化往往會(huì)出現(xiàn)奇異性即存在強(qiáng)間斷,例如激波、接觸間斷,或者稀疏波之類的弱間斷。

為了降低這種奇異性,提高間斷分辨率,近年來發(fā)展了許多的理論和計(jì)算方法,像譜方法、奇異攝動(dòng)理論、小波分析法[3]以及某些高分辨率的數(shù)值格式等。在早期,數(shù)值計(jì)算捕捉激波通常采取一階精度的差分格式,后來,Van Leer[4]將一階Godunov格式進(jìn)行推廣,率先提出了二階精度MUSCL(monotone upstream-centred schemes for conservation laws),隨后近40年,高精度、高分辨率數(shù)值格式蓬勃發(fā)展,出現(xiàn)諸如TVD、PPM、ENO、RKDG、CE/SE、WENO、WENO-Z及各類雜交格式等。高精度數(shù)值格式一般色散較強(qiáng),容易產(chǎn)生非物理性振蕩,通常需要額外的限制振蕩的方法,比如人工粘性法,可人工粘性的添加有時(shí)會(huì)使數(shù)值解的耗散比較嚴(yán)重,間斷的分辨率不夠清晰,其也非從根本上解決振蕩。

此外,提升Eulerian法數(shù)值求解精度的另一個(gè)角度則是網(wǎng)格加密。若采用均勻網(wǎng)格計(jì)算求解,便需對(duì)整個(gè)計(jì)算域進(jìn)行加密,使得計(jì)算資源極大地浪費(fèi),基于這個(gè)矛盾,網(wǎng)格自適應(yīng)技術(shù)應(yīng)運(yùn)而生。網(wǎng)格自適應(yīng)技術(shù)大致分三類,h-型(額外增加節(jié)點(diǎn))、p-型(增加逼近多項(xiàng)式的階數(shù))和r-型(移動(dòng)網(wǎng)格節(jié)點(diǎn)來達(dá)到網(wǎng)格重新分布),此外還有正在發(fā)展的結(jié)合性方法——h-p方法及h-r方法等,這里著重介紹r-型,也即移動(dòng)網(wǎng)格方法。移動(dòng)網(wǎng)格法發(fā)展至今,其所憑借的網(wǎng)格移動(dòng)策略及網(wǎng)格泛函逐漸多樣化,比如基于變量擴(kuò)展的等勢(shì)方法、調(diào)和映射、坐標(biāo)變換的雅可比矩陣思想、基于所謂等分布和對(duì)齊條件的泛函等。最近幾年,Luo等[5]又將DG(discontinuous galerkin)方法同移動(dòng)網(wǎng)格偏微分方程法相組合,提出了一種針對(duì)雙曲守恒律方程的擬拉格朗日移動(dòng)網(wǎng)格間斷Galerkin法,從舊網(wǎng)格到新網(wǎng)格的物理變量并不需要插值;Lopez等[6]還提出一種基于MMPDE法的并行變分網(wǎng)格改進(jìn)方案。本文中的偽弧長算法可歸結(jié)為r-型方法。

由于偽弧長算法涉及網(wǎng)格的自適應(yīng)移動(dòng),進(jìn)而導(dǎo)致原始均勻正交的物理空間發(fā)生扭曲變形,這給格式的重構(gòu)與插值帶來了困難。為了避免直接在變形的非結(jié)構(gòu)網(wǎng)格中重構(gòu)數(shù)值格式(需要大量的模板,計(jì)算效率較低),根據(jù)弧長映射關(guān)系,借助坐標(biāo)變換將物理空間映射至均勻正交的弧長計(jì)算空間,然后在弧長計(jì)算空間中,基于維數(shù)分裂思想可以用較少的模板完成高精格式的重構(gòu),從而保證了偽弧長算法的計(jì)算效率。偽弧長算法能夠?qū)⑽锢砜臻g中的強(qiáng)間斷問題轉(zhuǎn)變成均分弧長空間中正常的弱奇異流體問題,而計(jì)算空間的轉(zhuǎn)換更保證了原本數(shù)值格式的運(yùn)用,在此過程中還能巧妙地避開虛假振蕩。

算法程序的可靠性一直是數(shù)值模擬領(lǐng)域亟待解決的重難點(diǎn),諸如數(shù)理模型的簡(jiǎn)化、邊界條件的設(shè)置近似以及迭代格式的選取都有可能對(duì)計(jì)算結(jié)果造成偏差[7]。人為解方法(method of manufactured solution,MMS)早期經(jīng)典工作可現(xiàn)于Oberkampf、Trucano及Roy等。Blais等[8]提出了一種針對(duì)VANS方程(the volume-averaged Navier-stokes equations)的人為解方法框架以驗(yàn)證其算法程序,并能同任何CFD技術(shù)相結(jié)合。Choudhary等[9]介紹了基于旋度而滿足無散度約束的人為解方法,以驗(yàn)證兩相不可壓縮流控制方程。劉學(xué)哲等[10]利用其構(gòu)造的一類二維人為解模型針對(duì)性地驗(yàn)證輻射流體力學(xué)程序,該辦法中動(dòng)量、能量方程包含源項(xiàng)而質(zhì)量方程無源項(xiàng)。

本文中研究了二維情況下基于MUSCL的偽弧長算法模型建立過程,針對(duì)網(wǎng)格自適應(yīng)移動(dòng)造成的物理空間扭曲變形而不易構(gòu)造高精度格式的難題,基于坐標(biāo)變換的思想,將變形的物理空間映射至一套正交的弧長計(jì)算空間,從而在弧長計(jì)算空間中實(shí)現(xiàn)控制方程的求解與新網(wǎng)格守恒變量的插值過程,提高了計(jì)算效率。編寫了該算法相應(yīng)的一維、二維程序,在人為解方法基礎(chǔ)上驗(yàn)證其精度,借助某些經(jīng)典算例,將PALM數(shù)值結(jié)果對(duì)比分析于有限體積法。

1 偽弧長算法

本節(jié)偽弧長算法包含了2個(gè)部分。第一部分,借助有限體積法給出物理空間中控制方程的時(shí)空離散格式。如前描述,在考慮網(wǎng)格尺度影響的非均勻網(wǎng)格下,格式重構(gòu)不易,尤其在二維及其以上的情況,因此,第二部分在網(wǎng)格自適應(yīng)移動(dòng)之后,采取坐標(biāo)變換的策略將物理量映射至弧長計(jì)算空間中,并在該空間進(jìn)行格式重構(gòu),求解完之后再逆變換映射回原物理空間。

1.1 離散格式

考慮如下雙曲守恒系統(tǒng):

(1)

式中:w為質(zhì)量、動(dòng)量、能量組成的守恒變量;F(w),G(w),S(w)為w的函數(shù)。

式(1)在網(wǎng)格單元Ki, j上進(jìn)行積分,得到:

(2)

式中:dσ為面積微元。

(3)

式中:|Ki, j|為網(wǎng)格Ki, j的面積;?Ki, j為網(wǎng)格的邊界;ds為網(wǎng)格邊界長度微元;ni, j為邊界?Ki, j的單位外法向量;δ(w)=(F(w),G(w))。

數(shù)值通量借助局部Lax-Friedrichs格式,定義為:

(4)

h(u,v,n)=-h(u,v,-n),h(u,u,n)=δ(u)·n

(5)

式(3)寫成半離散格式,有:

(6)

圖1 網(wǎng)格單元計(jì)算示意圖Fig.1 Schematic diagram of grid cell calculation

(7)

式中:ψ為非線性限制器函數(shù)。在計(jì)算中ψ取[4]:

(8)

式中:ε為小量,可取ε=10-9。

對(duì)于時(shí)間的離散,忽略掉網(wǎng)格索引i,j,采取如下的三階TVD-RK格式:

(9)

1.2 網(wǎng)格自適應(yīng)移動(dòng)與計(jì)算空間的轉(zhuǎn)換

以x=(x,y)和ζ=(ξ,η)分別代表物理空間坐標(biāo)和弧長空間坐標(biāo),(xi-1, j-1,xi-1, j,xi, j,xi, j-1)為網(wǎng)格單元Ki, j的4個(gè)頂點(diǎn),從計(jì)算空間Ωc到物理空間Ωp的一一對(duì)應(yīng)坐標(biāo)映射為(x,y)=(x(ξ,η),y(ξ,η))。多維空間要同時(shí)照顧到網(wǎng)格的移動(dòng)速度、方向和尺寸問題,且期望網(wǎng)格朝物理量梯度大之處進(jìn)行移動(dòng),借助變分原理可知,網(wǎng)格的自適應(yīng)移動(dòng)將受泛函E(x,y)最小值影響[12]。

(10)

式中:G1,G2為給定的正定對(duì)稱矩陣;▽=(?ξ,?η)T,且有?ξ=?/?ξ,?η=?/?η,則式(10)對(duì)應(yīng)的Euler-Lagrange方程[13]為:

?ξ(G1?ξx)+?η(G1?ηx)=0

?ξ(G2?ξy)+?η(G2?ηy)=0

(11)

它對(duì)應(yīng)于泛函臨界點(diǎn)。

(12)

具體光滑次數(shù)視情況而定,這里不再贅述。令G=wI,如此,進(jìn)一步簡(jiǎn)化式(11),有:

▽·(φ▽x)=0

▽·(φ▽y)=0

(13)

對(duì)式(13)采用Gauss-Seidel迭代法進(jìn)行計(jì)算,具體步驟如下:

(14)

(15)

網(wǎng)格移動(dòng)完之后,需要得到物理量在新網(wǎng)格上的值,這一步稱作物理量重映,本文中采用基于通量的重映方法。如圖2二維新舊網(wǎng)格轉(zhuǎn)化圖所示。

圖2 二維新舊網(wǎng)格轉(zhuǎn)換示意圖Fig.2 Schematic diagram of old grid transforming to new grid in 2D

設(shè)Dk為經(jīng)過一次迭代物理空間中網(wǎng)格變化導(dǎo)致邊?kKi, j所掃過的面積,則有:

(16)

基于前人的工作[14],考慮到新網(wǎng)格x[k+1]與舊網(wǎng)格x[k]之間通量守恒,得到如下守恒插值格式:

(17)

Fk(m,n)=max{Dk,0}·n+min{Dk,0}·m

(18)

引入弧長坐標(biāo)空間變換式(1)原始雙曲守恒方程,得到:

(19)

式中:τ=t,U=ξt+uξx+vξy和V=ηt+uηx+vηy為弧長空間的逆變速度,ψ=w/J,J為空間轉(zhuǎn)換的Jacobian行列式,它的表達(dá)式為:

xξyη-xηyξ

(20)

式中:xξ,xη,yη,yξ,xτ,yτ為度量系數(shù),其中,xτ和yτ為網(wǎng)格的移動(dòng)速度,表達(dá)式見式(15)。Jacobian矩陣J的計(jì)算可借助幾何守恒條件來進(jìn)行,若能給出J的初值,J的推進(jìn)便可采取式(21)來計(jì)算:

(21)

xξ的計(jì)算過程如式(22),其余3個(gè)度量系數(shù)的計(jì)算相類似:

(22)

這樣通過式(19)求解完物理量之后,再進(jìn)行逆變換ψ=w/J,便可得到原始空間物理量。

綜上所述,可以給出基于二階MUSCL的偽弧長算法具體步驟:

第2步:求解弧長監(jiān)控函數(shù)φ,并根據(jù)式(12)進(jìn)行光滑濾波打磨。

第5步:控制方程的時(shí)間步迭代更新,通過式(19)、式(21)及式(22)得到弧長計(jì)算坐標(biāo)系中的控制方程,通過式(9)更新下一時(shí)刻物理量,求解完之后再根據(jù)逆變換將物理量映射回原物理空間。

第6步:若求解達(dá)到終點(diǎn)時(shí)間,程序終止,否則轉(zhuǎn)到第3步。

2 精度驗(yàn)證

出于方程跟問題的復(fù)雜性,精確解很難獲取,因此,這里采用人為解方法對(duì)一維、二維程序進(jìn)行驗(yàn)證。

2.1 一維情況

第1章的理論基于二維情況,對(duì)一維而言需降維處理計(jì)算,這里不再贅述。針對(duì)一維情況構(gòu)造一組人為解(算例中物理量為無量綱,下同),計(jì)算域取[0,2π],計(jì)算終止時(shí)間tend=2.0,初值條件為:

ρ(x,0)=1.0+0.2sin(x)

u(x,0)=1.0,p(x,0)=1.0

(23)

計(jì)算采用周期性邊界條件,則其精確解如下:

ρ(x,t)=1.0+0.2sin(x-t)

u(x,t)=1.0,p(x,t)=1.0

(24)

偽弧長控制函數(shù)取(a1,a2)=(10,20.25),對(duì)于一維、二維情況范數(shù)誤差公式及收斂階計(jì)算式如下[15]:

(25)

驗(yàn)證結(jié)果如表1所示。可以看出偽弧長算法并沒有因網(wǎng)格變形而損失太多的精度,隨著計(jì)算網(wǎng)格數(shù)量的增多,其L1誤差在逐漸減小,計(jì)算結(jié)果逐漸接近于精確解,并且,基于MUSCL的偽弧長算法收斂階也為二階。將偽弧長算法的誤差和收斂階同有限體積法的進(jìn)行比較,能看到相同網(wǎng)格數(shù)目下,PALM的L1誤差一般比固定網(wǎng)格的要小。

表1 有限體積法與偽弧長算法在不同網(wǎng)格數(shù)時(shí)的誤差和精度(1D)Table 1 Error and accuracy of FVM and PALM(1D)

2.2 二維情況

基于兩步化學(xué)反應(yīng)流模型,針對(duì)式(1)的雙曲守恒系統(tǒng),令:

(26)

式中:ωα,ωβ為化學(xué)反應(yīng)變化速率,且有:

(27)

式中:α,β分別代表未活化的物質(zhì)占所有物質(zhì)的比例和放熱反應(yīng)進(jìn)行的程度;Q為熱釋放率;R為氣體常數(shù);kα,kβ為反應(yīng)率常數(shù);Eα,Eβ為激活能。

狀態(tài)方程:

(28)

式中:γ為比熱比。計(jì)算域取[0,2π]×[0,2π],計(jì)算終止時(shí)間tend=2π,給定初值條件:

ρ(x,0)=1.0+0.5sin(x+y)

u(x,0)=v(x,0)=1.0,p(x,0)=1.0

α(x,0)=0.5+0.5sin(x+y)

β(x,0)=0.5+0.5sin(x+y)

(29)

計(jì)算采用周期性邊界條件,其精確解如下:

ρ(x,t)=1.0+0.5sin(x+y-t)

u(x,t)=v(x,t)=1.0,p(x,t)=1.0

α(x,t)=0.5+0.5sin(x+y-t)

β(x,0)=0.5+0.5sin(x+y-t)

(30)

偽弧長控制函數(shù)取(a1,a2)=(6,16)。

驗(yàn)證結(jié)果如表2所示。二維情況下偽弧長算法比之固定網(wǎng)格方法同樣能適當(dāng)提高計(jì)算精度并減小L1誤差。須知,PALM并未增加網(wǎng)格節(jié)點(diǎn),它是根據(jù)物理量的分布狀況致使網(wǎng)格自適應(yīng)移動(dòng)變形,從而在物理量梯度大的地方集聚了較多網(wǎng)格,減小了總體均值誤差。隨著網(wǎng)格單元數(shù)量增多,PALM的L1誤差在逐漸減小,且其收斂速度較快,收斂階更接近于2。

表2 有限體積法與偽弧長算法在不同網(wǎng)格數(shù)時(shí)的誤差和精度(2D)Table 2 Error and accuracy of FVM and PALM(2D)

3 數(shù)值算例

下面基于偽弧長算法開展算例驗(yàn)證。

3.1 Sod激波管問題

初值條件:

(31)

計(jì)算域?yàn)閇-5,5],計(jì)算終止時(shí)間tmax=2.0,邊界條件為自由輸出邊界條件,γ取1.4,偽弧長控制函數(shù)為(a1,a2)=(1.7,200),網(wǎng)格數(shù)N=150。計(jì)算結(jié)果如圖3所示,分別以Palm、Fixed來代表偽弧長算法和固定網(wǎng)格下有限體積法的數(shù)值標(biāo)識(shí)。

該算例參照解(Reference)在5 000個(gè)固定網(wǎng)格下獲得。對(duì)比有限體積法,偽弧長算法能以較少的網(wǎng)格數(shù)獲得更好的界面分辨率,在相同分辨率的要求下,PALM顯得更加高效;同時(shí)也發(fā)現(xiàn),要想提高計(jì)算精度,相當(dāng)數(shù)目的網(wǎng)格量是必須的。圖3(b)給出了移動(dòng)網(wǎng)格軌跡線圖,表征偽弧長算法很好地實(shí)現(xiàn)了對(duì)奇異間斷處的自適應(yīng)捕捉,在弧長參數(shù)作用下,網(wǎng)格點(diǎn)可聚集在激波、接觸間斷處甚至是稀疏波的頭部與尾部。

圖3 Sod激波管圖Fig.3 Results of Sod shock tube

3.2 一維爆炸波

該雙爆轟波碰撞問題初值條件如下:

(32)

計(jì)算域取[0,1],計(jì)算終止時(shí)間為tmax=0.038,置CFL系數(shù)為0.8,邊界條件為反射邊界,γ取1.4。計(jì)算網(wǎng)格數(shù)N=250,偽弧長控制函數(shù)以兩套不同的參數(shù)進(jìn)行比較,分別為:(2,20)和(2,200),在圖4中標(biāo)識(shí)為Palm-2-1和Palm-2-2。由于稀疏波的作用,計(jì)算過程中可能會(huì)出現(xiàn)負(fù)密度、負(fù)壓力,從而導(dǎo)致程序終止,因此對(duì)該算法附加了保正性條件,詳細(xì)理論見文獻(xiàn)[16]。

此算例參照解在10 000個(gè)固定網(wǎng)格下獲得。圖4將3種數(shù)值結(jié)果同參照解進(jìn)行比較,說明在網(wǎng)格總數(shù)較少情況下偽弧長算法的數(shù)值結(jié)果顯著優(yōu)于有限體積法,而固定網(wǎng)格算法數(shù)值耗散較大,對(duì)極值點(diǎn)的捕捉能力很弱。在相同分辨率的要求下,PALM比固定網(wǎng)格法更為高效。圖5的網(wǎng)格軌跡線刻畫出沖擊波在t=0.027時(shí)刻左右相撞,之后又沿2個(gè)方向繼續(xù)傳播,偽弧長算法對(duì)解變化劇烈區(qū)域之模擬效果更為逼真。另外能夠驗(yàn)證,不同偽弧長監(jiān)控函數(shù)模型對(duì)激波的捕捉能力存在差異,針對(duì)本算例第二套系數(shù)刻畫的分辨率的確要優(yōu)于第一套,其網(wǎng)格移動(dòng)得更劇烈。

圖4 爆炸波問題數(shù)值結(jié)果比較Fig.4 Numerical results of blast-wave problem

圖5 不同參數(shù)下網(wǎng)格軌跡圖Fig.5 Grid trajectories under different parameters

3.3 二維Riemann問題

在計(jì)算域[0,1]×[0,1]內(nèi),二維Riemann問題初始分布如圖6所示,其邊界條件均為流入流出條件。

圖6 二維Riemann問題初始區(qū)域分布圖Fig.6 Initial region distribution of two-dimensional Riemann problem

對(duì)于這樣一個(gè)Riemann問題——各有2個(gè)正負(fù)向滑移線的接觸間斷,初值條件為:

(33)

計(jì)算終止時(shí)間為tmax=0.3,偽弧長控制函數(shù)取(a1,a2)=(2.7,1 000)。數(shù)值結(jié)果如圖7所示。

圖7 二維Riemann問題在t=0.3時(shí)刻的密度云圖和網(wǎng)格自適應(yīng)分布圖Fig.7 Density cloud map and mesh adaptive distribution map of two-dimensional Riemann problem at t=0.3

借助波傳播的形態(tài)可以發(fā)現(xiàn),采用200×200的固定網(wǎng)格算法,其模擬的精度較差,密度分布梯度線比較粗糙,而偽弧長算法能更銳利地追蹤梯度界面。雖然網(wǎng)格自適應(yīng)移動(dòng)會(huì)導(dǎo)致額外的局部時(shí)間迭代,但移動(dòng)網(wǎng)格的計(jì)算代價(jià)比之h-型方法直接加密網(wǎng)格的代價(jià)要小,對(duì)于求解較大尺度問題,在保證精度的同時(shí),還應(yīng)減少CPU運(yùn)行時(shí)間,合適的弧長參數(shù)下,偽弧長算法的模擬效率要優(yōu)于固定網(wǎng)格方法。

3.4 雙馬赫反射問題

該算例是一個(gè)波系結(jié)構(gòu)十分復(fù)雜的流場(chǎng)問題。為簡(jiǎn)化模型,將計(jì)算域設(shè)置在規(guī)則區(qū)域當(dāng)中:[0,4]×[0,1]。其中底部(x>1/6為固壁邊界,x<1/6為來流)一Mach數(shù)為10的斜強(qiáng)激波擱置于x=1/6,y=0處,與x軸成60°角,其他壁面為反射邊界條件,模型介紹這里不再贅述。初值條件:

(34)

計(jì)算終止時(shí)間tmax=0.2,偽弧長控制函數(shù)取(a1,a2)=(1,15),該問題在320×80網(wǎng)格數(shù)下計(jì)算,如圖8所示。

圖8 雙馬赫反射問題在t=0.2時(shí)刻的密度云圖和網(wǎng)格自適應(yīng)分布圖(包含局部細(xì)節(jié))Fig.8 Density cloud diagram and grid adaptive distribution diagram (including local details) of dual Mach reflection problem at t=0.2

能夠驗(yàn)證二者在流場(chǎng)的大尺度現(xiàn)象(例如跳變點(diǎn)、附壁射流、馬赫桿及大致滑移面等)上模擬得相差無幾,而圖8(b)能更真實(shí)地反映流場(chǎng)內(nèi)波傳播的強(qiáng)度與小尺度結(jié)構(gòu),諸如沿滑移線的小漩渦匯聚現(xiàn)象、第二個(gè)三波點(diǎn)等。從局部放大圖可以看到,相對(duì)固定網(wǎng)格的計(jì)算,偽弧長算法對(duì)渦核附近以及x=2.5激波相匯區(qū)域的計(jì)算更為細(xì)化。PALM更好地模擬了高馬赫數(shù)的傳播與反射問題。

4 結(jié)論

1) 相對(duì)于固定網(wǎng)格方法,偽弧長算法在提高爆轟波強(qiáng)間斷分辨率上彰顯了很好的優(yōu)越性,本文中基于坐標(biāo)變換策略發(fā)展的弧長算法理論,讓爆轟波的求解刻畫有效地避開雙曲守恒系統(tǒng)的奇異性問題,使得MUSCL格式在弧長計(jì)算空間中得到很好地運(yùn)用,提高了數(shù)值求解的分辨率。

2) 對(duì)比分析不同偽弧長控制函數(shù)對(duì)計(jì)算精度的影響,驗(yàn)證了監(jiān)控函數(shù)模型的選擇將直接影響數(shù)值模擬的效果,弧長參數(shù)決定網(wǎng)格移動(dòng)的合理性,網(wǎng)格移動(dòng)得越劇烈,間斷附近的分辨率就越高,但可能會(huì)犧牲光滑區(qū)域的分辨率且相應(yīng)增加了迭代步計(jì)算。所以采取哪類網(wǎng)格細(xì)化準(zhǔn)則,構(gòu)造怎樣的偽弧長控制函數(shù)(比如多物理量耦合決定、復(fù)雜多項(xiàng)式形式等)以及如何選擇可調(diào)參數(shù),從網(wǎng)格生成這方面來說是一個(gè)系統(tǒng)問題,特別在高維空間之時(shí)。