基于坐標變換的強間斷問題偽弧長算法

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

(北京理工大學 爆炸科學與技術國家重點實驗室, 北京 100081)

0 引言

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

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

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

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

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

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

1 偽弧長算法

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

1.1 離散格式

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

(1)

式中:w為質量、動量、能量組成的守恒變量;F(w),G(w),S(w)為w的函數。

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

(2)

式中:dσ為面積微元。

(3)

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

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

(4)

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

(5)

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

(6)

圖1 網格單元計算示意圖Fig.1 Schematic diagram of grid cell calculation

(7)

式中:ψ為非線性限制器函數。在計算中ψ取[4]:

(8)

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

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

(9)

1.2 網格自適應移動與計算空間的轉換

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

(10)

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

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

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

(11)

它對應于泛函臨界點。

(12)

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

▽·(φ▽x)=0

▽·(φ▽y)=0

(13)

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

(14)

(15)

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

圖2 二維新舊網格轉換示意圖Fig.2 Schematic diagram of old grid transforming to new grid in 2D

設Dk為經過一次迭代物理空間中網格變化導致邊?kKi, j所掃過的面積,則有:

(16)

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

(17)

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

(18)

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

(19)

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

xξyη-xηyξ

(20)

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

(21)

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

(22)

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

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

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

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

第6步:若求解達到終點時間,程序終止,否則轉到第3步。

2 精度驗證

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

2.1 一維情況

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

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

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

(23)

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

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

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

(24)

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

(25)

驗證結果如表1所示?？梢钥闯鰝位￠L算法并沒有因網格變形而損失太多的精度,隨著計算網格數量的增多,其L1誤差在逐漸減小,計算結果逐漸接近于精確解,并且,基于MUSCL的偽弧長算法收斂階也為二階。將偽弧長算法的誤差和收斂階同有限體積法的進行比較,能看到相同網格數目下,PALM的L1誤差一般比固定網格的要小。

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

2.2 二維情況

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

(26)

式中:ωα,ωβ為化學反應變化速率,且有:

(27)

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

狀態(tài)方程:

(28)

式中:γ為比熱比。計算域取[0,2π]×[0,2π],計算終止時間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)

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

ρ(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)

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

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

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

3 數值算例

下面基于偽弧長算法開展算例驗證。

3.1 Sod激波管問題

初值條件:

(31)

計算域為[-5,5],計算終止時間tmax=2.0,邊界條件為自由輸出邊界條件,γ取1.4,偽弧長控制函數為(a1,a2)=(1.7,200),網格數N=150。計算結果如圖3所示,分別以Palm、Fixed來代表偽弧長算法和固定網格下有限體積法的數值標識。

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

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

3.2 一維爆炸波

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

(32)

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

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

圖4 爆炸波問題數值結果比較Fig.4 Numerical results of blast-wave problem

圖5 不同參數下網格軌跡圖Fig.5 Grid trajectories under different parameters

3.3 二維Riemann問題

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

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

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

(33)

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

圖7 二維Riemann問題在t=0.3時刻的密度云圖和網格自適應分布圖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的固定網格算法,其模擬的精度較差,密度分布梯度線比較粗糙,而偽弧長算法能更銳利地追蹤梯度界面。雖然網格自適應移動會導致額外的局部時間迭代,但移動網格的計算代價比之h-型方法直接加密網格的代價要小,對于求解較大尺度問題,在保證精度的同時,還應減少CPU運行時間,合適的弧長參數下,偽弧長算法的模擬效率要優(yōu)于固定網格方法。

3.4 雙馬赫反射問題

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

(34)

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

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

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

4 結論

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

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