一種三階有限體積法及其在欠膨脹射流激波結(jié)構(gòu)數(shù)值模擬中的應(yīng)用*

謝 政,謝 建,李 良

(第二炮兵工程大學(xué)發(fā)射理論與技術(shù)軍隊重點(diǎn)實驗室,陜西西安710025)

一種三階有限體積法及其在欠膨脹射流激波結(jié)構(gòu)數(shù)值模擬中的應(yīng)用*

謝 政,謝 建,李 良

(第二炮兵工程大學(xué)發(fā)射理論與技術(shù)軍隊重點(diǎn)實驗室,陜西西安710025)

以噴管出口欠膨脹射流為研究對象,在Lagrange坐標(biāo)系下建立欠膨脹射流二維積分形式的流動方程。通過在單元交接面處進(jìn)行三階ENO(essentially nonoscillatory)格式插值,構(gòu)造得到一種適用于求解該方程的三階ENO有限體積法。采用該格式對一維Sod激波管算例和噴管出口欠膨脹射流進(jìn)行數(shù)值計算。計算結(jié)果表明,該方法具有高精度、基本無振蕩的特點(diǎn),能很好地捕捉包含激波、滑移線以及三波交點(diǎn)等復(fù)雜流場波系結(jié)構(gòu)。計算得到的波系結(jié)構(gòu)中馬赫盤的位置與實驗結(jié)果吻合很好,相對誤差小于1.1%。

激波結(jié)構(gòu);基本無振蕩;有限體積法;欠膨脹射流

在槍炮武器、火箭和航空航天等工程技術(shù)領(lǐng)域都涉及燃?xì)馍淞鳑_擊的問題,射流波系結(jié)構(gòu)的影響因素有很多[1]。為了減少燃?xì)馍淞鲗鼒鲈O(shè)備和工程設(shè)施的沖擊危害,需要了解射流流場波系結(jié)構(gòu)和流場的流動狀態(tài)。由于相關(guān)的射流實驗費(fèi)用昂貴,且通過實驗不可能詳盡地研究各種因素在不同水平對射流波系結(jié)構(gòu)的影響。然而,采用高精度的數(shù)值方法能夠有效地捕捉燃?xì)馍淞鲌龅募げúㄏ到Y(jié)構(gòu),并且與相同工況下的實驗結(jié)果能很好地吻合[2-3]。

近十幾年來,基于非結(jié)構(gòu)/結(jié)構(gòu)網(wǎng)格的高精度算法發(fā)展迅速,主要有TVD、DGM、ENO、k-exact有限體積方法等[4-6]。從A.Harten等提出ENO格式以后,許多學(xué)者從不同的思路出發(fā),提出了多種形式的ENO格式,如WENO、CENO[7-8]。徐文燦等[9]采用高分辨率的ENO格式對噴管流動問題進(jìn)行了數(shù)值計算,得到了推力隨擺角變化的規(guī)律,與實驗結(jié)果基本一致,證明了采用ENO格式能夠很好地捕獲復(fù)雜波系,反映激波和邊界層之間的相互干擾。陸霄露等[10]也采用ENO格式計算了一維非定常進(jìn)排氣流動問題,計算結(jié)果表明發(fā)動機(jī)的主要性能參數(shù)都和實驗結(jié)果吻合很好。

為了清晰地捕捉到流場中的間斷區(qū)域(激波結(jié)構(gòu)),研究一種有效求解Lagrange坐標(biāo)系中積分形式歐拉方程的三階ENO有限體積法,并且采用具有精確解的激波管算例驗證該方法的有效性。最后,采用該方法求解典型的欠膨脹射流的流動問題。

1 控制方程

在Lagrange坐標(biāo)系中,二維軸對稱非穩(wěn)態(tài)可壓縮氣體流動的積分形式歐拉方程[11]為:

式中:Ω(t)為由邊界Γ(t)包圍的運(yùn)動控制體,守恒變量矢量U=[ρ,M,E]T,通量矢量F=[0,p n,p u·n]T。其中,ρ為密度,u為速度矢量,M=(Mx,My)=ρu為動量,E為總能量,n=(nx,ny)為邊界Γ(t)的外法線向量。理想氣體的狀態(tài)方程壓力p=(γ-1)ρe,e為氣體的質(zhì)量內(nèi)能。

2 數(shù)值計算方法

二維空間域劃分為m×n個計算單元。Ωi+1/2,j+1/2為一個四邊形的計算單元,它的頂點(diǎn)分別為(xi,j,yi,j)、(xi+1,j,yi+1,j)、(xi+1,j+1,yi+1,j+1)、(xi,j+1,yi,j+1)。Si+1/2,j+1/2為計算單元的面積。采用格心型有限體積法,所有的物理量都存儲在計算單元的中心。因此,采用均值形式的密度動量以及總能量分別如下式所示:

式中:Mx、My分別為x方向和y方向的動量。

采用有限體積法,對控制方程式(1)進(jìn)行離散,得到軸對稱Euler方程組的半離散守恒方程形式:

式中:k為計算單元4個頂點(diǎn)的序號,按逆時針方向排序;pk,k+1為計算單元k頂點(diǎn)和k+1頂點(diǎn)邊界上的壓強(qiáng)。由于采用有限體積方法得到的是網(wǎng)格平均值由可以構(gòu)造高階插值多項式p(x,y),但不能保證p(x,y)在網(wǎng)格邊界上是連續(xù)的,因此pk,k+1不能直接從構(gòu)造的插值多項式p(x,y)求得。這里需要解一個Riemann問題,通過求精確Riemann解,可以得到pk,k+1的值。而法向速度un和網(wǎng)格側(cè)面積的乘積取為:

式中:uk為計算單元Ωi+1/2,j+1/2中序號為k的頂點(diǎn)在x方向的速度。同理,vk為y方向的速度,(xk,yk)為該頂點(diǎn)的坐標(biāo)值。Δt時間后網(wǎng)格中心新速度和能量^E可由式(4)、(5)聯(lián)立求得。計算單元頂點(diǎn)的新速度采用格點(diǎn)型格式直接計算,構(gòu)造頂點(diǎn)處速度的控制體為圖1中虛線邊界包圍的控制體在二維空間上進(jìn)行三階ENO重構(gòu),則有插值多項式:

圖1 速度的控制體Fig.1 Control volume of velocity

式中:q0、qx、qy、qxx、qxy、qyy為待定系數(shù)。將邊界的線段用參數(shù)方程x=x1+(x2-x1)t,y=y1+(y2-y1)t,根據(jù)文獻(xiàn)[10],計算單元頂點(diǎn)運(yùn)動的位置根據(jù)每個頂點(diǎn)的新位置求得計算單元新的面積進(jìn)而得到新的密度為計算單元內(nèi)流體質(zhì)量。

聯(lián)立式(2)～(4),在結(jié)構(gòu)網(wǎng)格下采用高階ENO-FVM格式,方程(1)能夠表示[4]為:

求解過程忽略3階以上高階項的影響,方程(1)的求解轉(zhuǎn)化為常微分方程的求解問題。常微分方程的求解采用具有TVD性質(zhì)的三階Runge-Kutta方法,如下式[10]所示:

式中:時間步長Δt根據(jù)式來確定,其中,r表示第r個計算步;邊界中最短邊的長度為計算單元Ωi+1/2,j+1/2內(nèi)的聲速;計算中Courant數(shù)λ取為0.5。

3 計算結(jié)果與分析

3.1 Sod激波管算例

計算域為[0,2],所采用的網(wǎng)格數(shù)為100,初始條件為:

采用Newmann邊界條件,取CFL系數(shù)為0.5,計算推進(jìn)到t=0.5 s時終止計算。此時,流場中包含一個激波、一個接觸間斷和一個膨脹波。圖2給出了采用數(shù)值方法得到的密度分布曲線。從圖2可以看出,三階ENO有限體積法的激波分辨率很高,能準(zhǔn)確捕捉到激波結(jié)構(gòu),將激波被抹平的厚度限制在一個網(wǎng)格單元尺度。表1給出了不同網(wǎng)格尺度下,三階ENO有限體積法的密度誤差L1范數(shù)及其收斂精度[6]。從表1可以看出,隨著網(wǎng)格尺度的減小,計算方法體現(xiàn)出了網(wǎng)格收斂性,并且收斂的精度大致相當(dāng),計算時間以指數(shù)倍增長。

圖2 Sod激波管密度分布Fig.2 Density profiles for the Sod problem

3.2 噴管出口欠膨脹射流問題

根據(jù)文獻(xiàn)[1],計算了噴管出口欠膨脹射流例子,計算域為噴管子午面截得的一半?yún)^(qū)域,如圖3所示。圖中AB表示噴管的出口半徑de,AC表示噴管的中心軸線,上游邊界BD和外邊界DE為靜止大氣條件,下游邊界CE采用外推。采用無量綱格式計算,參考長度de=10 mm,參考壓力p0= 101.325 k Pa,參考溫度T=300 K。計算區(qū)域為[0,9]×[0,4]計算區(qū)域內(nèi)布置270×160個網(wǎng)格,整個計算域初始時為靜止大氣條件,在t=0時刻射流突然噴出。計算了兩種工況,工況1輕度欠膨脹,噴管出口壓力比為2.06;工況2重度欠膨脹,噴管出口壓力比為26.4。

表1 三階ENO有限體積法密度誤差Table 1 Density error for third-order ENO finite volume method

圖3 計算區(qū)域Fig.3 Computational domain

圖4 工況1下計算結(jié)果等值線圖和紋影照片F(xiàn)ig.4 Contour maps and schlieren photograph in case 1

圖5 工況2下計算結(jié)果等值線圖和紋影照片F(xiàn)ig.5 Contour maps and schlieren photograph in case 2

圖4(a)～(c)分別給出了工況1在第1 000個計算時間步時刻的馬赫數(shù)、壓力、密度等值線分布。從圖4可以看出,射流的流場結(jié)構(gòu)類似鉆石狀,馬赫盤略有呈現(xiàn),入射激波、反射激波、馬赫盤等構(gòu)成的波胞交替產(chǎn)生。在下游邊界處,由于射流與靜止大氣之間的對流作用,流場中產(chǎn)生了間斷面。由于Euler方程可以看作是高Reynolds數(shù)下的近似,因此邊界層區(qū)域附近將出現(xiàn)Kelvin-Helmholtz不穩(wěn)定性,越往下游不穩(wěn)定性越明顯。計算結(jié)果與正格式數(shù)值計算的結(jié)果[12]和實驗紋影照片圖4(d)反映的流場結(jié)構(gòu)吻合較好。

圖5(a)～(c)分別為工況2下馬赫數(shù)、壓力和密度等值線分布。與工況1相比,由于噴管出口壓力比增大,圖5中有明顯的馬赫盤結(jié)構(gòu),在馬赫盤的邊緣與入射激波和反射激波交匯形成三叉激波,出現(xiàn)三波交點(diǎn)和滑移線,馬赫盤上游射流的結(jié)構(gòu)較穩(wěn)定,下游流場出現(xiàn)了明顯的Kelvin-Helmholtz不穩(wěn)定性。由于出口壓力比增大,計算域內(nèi)只能形成一個波胞,并且馬赫盤的過渡區(qū)間很窄,說明本文中采用的方法具有較強(qiáng)的捕捉激波的能力。圖5(d)為相同流動條件下的實驗紋影照片[13],實驗結(jié)果中第一個馬赫盤距噴管出口平面的距離為4.56de。采用數(shù)值方法得到的馬赫盤的位置為4.51de,文獻(xiàn)[12]中采用正格式的結(jié)果為4.68de,兩種算法相比較,采用本文方法得到的馬赫盤位置更精確。與實驗得到的馬赫盤位置比較,本文計算結(jié)果的誤差小于1.1%,且與實驗紋影照片所反映的流場結(jié)構(gòu)吻合較好,表明該方法具有較高的精度,能精確捕捉激波位置,并且在激波面上不會產(chǎn)生振蕩或抹平間斷現(xiàn)象。

4 結(jié) 論

以ENO格式為基礎(chǔ),通過在單元交界面處進(jìn)行三階ENO格式插值,構(gòu)造了一類求解Lagrange坐標(biāo)系中積分形式歐拉方程的三階ENO有限體積法。數(shù)值計算結(jié)果表明,該方法具有較高的數(shù)值精度,適用于非穩(wěn)態(tài)流場的數(shù)值計算,并且具有較強(qiáng)的激波捕捉能力,能夠清晰地模擬出復(fù)雜流場的射流結(jié)構(gòu),與相同工況下實驗結(jié)果吻合較好。

[1] Matsuda T,Umeda Y,Ishii R,et al.Numerical and experimental studies on chocked under-expanded jets[C]∥19th AIAA,Fluid Dynamics,Plasma Dynamics,and Lasers Conference.Honolulu,HI,USA,1987,7:87-1378-281.

[2] 劉小軍,傅德彬,牛青林,等.燃?xì)馍淞鳑_擊傳熱特性的數(shù)值模擬[J].爆炸與沖擊,2015,35(2):229-235. Liu Xiaojun,Fu Debin,Niu Qinlin,et al.Numerical simulation of heat transfer for exhausted gas jet impinging [J].Explosion and Shock Waves,2015,35(2):229-235.

[3] 薛曉春,余永剛,張琦.雙股燃?xì)馍淞髟诔湟菏覂?nèi)擴(kuò)展特性的實驗研究[J].爆炸與沖擊,2013,33(5):449-455. Xue Xiaochun,Yu Yonggang,Zhang Qi.Experimental study on expansion characteristics of twin combustion-gas jets in liquid filled chamber[J].Explosion and Shock Waves,2013,33(5):449-455.

[4] Ivan L,Groth C P T.High-order central ENO finite-volume scheme with adaptive mesh refinement[C]∥18th AIAA Computational Fluid Dynamics Conference.Miami,Florida,2007.

[5] Luo H,Luo L P,Nourgaliev R,et al.A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids[J].Journal of Computational Physics,2010,229(19):6961-6978.

[6] 范進(jìn)之,李樺.高精度有限體積法與間斷有限元法的比較[J].國防科技大學(xué)學(xué)報,2014,36(5):33-38. Fan Jinzhi,Li Hua.Comparison of high-precision finite volume method and discontinuous Galerkin method[J]. Journal of National University of Defense Technology,2014,36(5):33-38.

[7] Harten A,Enquist B,Osher S,et al.Uniformly high order essentially non-oscillatory schemes[J].Journal of Computational Physics,1987,71(2):231-303.

[8] 程曉晗,封建湖,聶玉峰.求解雙曲守恒方程的WENO型熵相容格式[J].爆炸與沖擊,2014,34(4):501-507. Cheng Xiaohan,Feng Jianhu,Nie Yufeng.WENO type entropy consistent scheme for hyperbolic conservation laws [J].Explosion and Shock Waves,2014,34(4):501-507.

[9] 徐文燦,黃振宇.高精度ENO格式在射流數(shù)值模擬中的應(yīng)用[J].空氣動力學(xué)學(xué)報,2001,19(4):401-406.Xu Wencan,Huang Zhenyu.Flow field calculation with high resolution ENO[J].Acta Aerodynamica Sinica, 2001,19(4):401-406.

[10] 陸霄露,鄧康耀.進(jìn)排氣一維非定常流動的基本無振蕩有限體積法的研究[J].內(nèi)燃機(jī)工程,2013,34(2):52-61. Lu Xiaolu,Deng Kangyao.Study of essentially non-oscillatory finite method for one-dimension unsteady intake and exhaust flows[J].Chinese Internal Combustion Engine Engineering,2013,34(2):52-61.

[11] Wang Yongjian,Zhao Ning,Wang Donghong,et al.A kind essentially non-oscillatory finite volume scheme in Lagrangian coordinates[J].Journal on Numerical Methods and Computer Application,2007,28(4):250-259.

[12] 朱孫科,陳二云,馬大為,等.燃?xì)庾杂缮淞鞯恼袷綌?shù)值模擬[J].空氣動力學(xué)報,2011,29(3):380-384. Zhu Sunke,Chen Eryun,Ma Dawei,et al.Numerical simulation of gas free jet by positive schemes[J].Acta Aerodynamica Sinica,2011,29(3):380-384.

[13] Ruggles A J,Ekoto I W.Experimental investigation of nozzle aspect ratio effects on under-expanded hydrogen jet release characteristics[J].International Journal of Hydrogen Energy,2014,39(35):20331-20338.

A three-order finite volume method and its application to under-expanded jet shock wave structure simulation

Xie Zheng,Xie Jian,Li Liang
(Military Key Laboratory for Armament Launch Theory&Technology,The Second Artillery Engineering University,Xi’an710025,Shaanxi,China)

By considering the under-expanded jet flow from nozzle exit,the integral form Euler equations for unsteady compressible flow in the Lagrange coordinates of a moving control volume was developed.By using three-order essentially non-oscillatory(ENO)interpolations at cell interfaces,a three-order ENO finite volume method for the integral form Euler equations was presented.The Sod shock tube case and nozzle outlet under-expanded jet shock wave structures were used to test the proposed scheme.The numerical results demonstrate that the method is accurate and non-oscillatory, and it can capture the wave structures of jet flow fields including shock cell structure,slip lines,jet boundary and the triple point well.Meanwhile,the simulated Mach disk locations in wave structures coincide with the experimentally measured ones,especially the error of the first Mach disk locations in wave structures between the numerical results and the experimental results was less than 1.1%.

shock wave structures;essentially non-oscillatory;finite volume method;under-expanded jet

O354;V211國標(biāo)學(xué)科代碼:13025

:A

10.11883/1001-1455(2017)02-0347-06

(責(zé)任編輯 張凌云)

2015-10-26;

:2016-03-31

國家自然科學(xué)基金項目(51475462)

謝 政(1989— ),男,博士研究生,xiez19891121@163.com。