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

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

0 引言

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

為了降低這種奇異性,提高間斷分辨率,近年來發(fā)展了許多的理論和計(jì)算方法,像譜方法、奇異攝動理論、小波分析法[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ù)值解的耗散比較嚴(yán)重,間斷的分辨率不夠清晰,其也非從根本上解決振蕩。

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

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

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

1 偽弧長算法

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

1.1 離散格式

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

(1)

式中:w為質(zhì)量、動量、能量組成的守恒變量;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。

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

(9)

1.2 網(wǎng)格自適應(yī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個頂點(diǎn),從計(jì)算空間Ωc到物理空間Ωp的一一對應(yīng)坐標(biāo)映射為(x,y)=(x(ξ,η),y(ξ,η))。多維空間要同時照顧到網(wǎng)格的移動速度、方向和尺寸問題,且期望網(wǎng)格朝物理量梯度大之處進(jìn)行移動,借助變分原理可知,網(wǎng)格的自適應(yīng)移動將受泛函E(x,y)最小值影響[12]。

(10)

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

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

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

(11)

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

(12)

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

▽·(φ▽x)=0

▽·(φ▽y)=0

(13)

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

(14)

(15)

網(wǎ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á)式見式(15)。Jacobian矩陣J的計(jì)算可借助幾何守恒條件來進(jìn)行,若能給出J的初值,J的推進(jìn)便可采取式(21)來計(jì)算:

(21)

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

(22)

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

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

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

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

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

2 精度驗(yàn)證

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

2.1 一維情況

第1章的理論基于二維情況,對一維而言需降維處理計(jì)算,這里不再贅述。針對一維情況構(gòu)造一組人為解(算例中物理量為無量綱,下同),計(jì)算域取[0,2π],計(jì)算終止時間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),對于一維、二維情況范數(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ù)時的誤差和精度(1D)Table 1 Error and accuracy of FVM and PALM(1D)

2.2 二維情況

基于兩步化學(xué)反應(yīng)流模型,針對式(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ì)算終止時間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)移動變形,從而在物理量梯度大的地方集聚了較多網(wǎng)格,減小了總體均值誤差。隨著網(wǎng)格單元數(shù)量增多,PALM的L1誤差在逐漸減小,且其收斂速度較快,收斂階更接近于2。

表2 有限體積法與偽弧長算法在不同網(wǎng)格數(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ì)算終止時間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)識。

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

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

3.2 一維爆炸波

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

(32)

計(jì)算域取[0,1],計(jì)算終止時間為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)識為Palm-2-1和Palm-2-2。由于稀疏波的作用,計(jì)算過程中可能會出現(xiàn)負(fù)密度、負(fù)壓力,從而導(dǎo)致程序終止,因此對該算法附加了保正性條件,詳細(xì)理論見文獻(xiàn)[16]。

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

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

(33)

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

圖7 二維Riemann問題在t=0.3時刻的密度云圖和網(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ǎo)致額外的局部時間迭代,但移動網(wǎng)格的計(jì)算代價比之h-型方法直接加密網(wǎng)格的代價要小,對于求解較大尺度問題,在保證精度的同時,還應(yīng)減少CPU運(yùn)行時間,合適的弧長參數(shù)下,偽弧長算法的模擬效率要優(yōu)于固定網(wǎng)格方法。

3.4 雙馬赫反射問題

該算例是一個波系結(jié)構(gòu)十分復(fù)雜的流場問題。為簡化模型,將計(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ì)算終止時間tmax=0.2,偽弧長控制函數(shù)取(a1,a2)=(1,15),該問題在320×80網(wǎng)格數(shù)下計(jì)算,如圖8所示。

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

4 結(jié)論

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

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