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

謝 政,謝 建,李 良

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

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

謝 政,謝 建,李 良

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

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

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

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

近十幾年來(lái),基于非結(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格式對(duì)噴管流動(dòng)問(wèn)題進(jìn)行了數(shù)值計(jì)算,得到了推力隨擺角變化的規(guī)律,與實(shí)驗(yàn)結(jié)果基本一致,證明了采用ENO格式能夠很好地捕獲復(fù)雜波系,反映激波和邊界層之間的相互干擾。陸霄露等[10]也采用ENO格式計(jì)算了一維非定常進(jìn)排氣流動(dòng)問(wèn)題,計(jì)算結(jié)果表明發(fā)動(dòng)機(jī)的主要性能參數(shù)都和實(shí)驗(yàn)結(jié)果吻合很好。

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

1 控制方程

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

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

2 數(shù)值計(jì)算方法

二維空間域劃分為m×n個(gè)計(jì)算單元。Ωi+1/2,j+1/2為一個(gè)四邊形的計(jì)算單元,它的頂點(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為計(jì)算單元的面積。采用格心型有限體積法,所有的物理量都存儲(chǔ)在計(jì)算單元的中心。因此,采用均值形式的密度動(dòng)量以及總能量分別如下式所示:

式中:Mx、My分別為x方向和y方向的動(dòng)量。

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

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

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

圖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],計(jì)算單元頂點(diǎn)運(yùn)動(dòng)的位置根據(jù)每個(gè)頂點(diǎn)的新位置求得計(jì)算單元新的面積進(jìn)而得到新的密度為計(jì)算單元內(nèi)流體質(zhì)量。

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

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

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

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

3.1 Sod激波管算例

計(jì)算域?yàn)閇0,2],所采用的網(wǎng)格數(shù)為100,初始條件為:

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

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

3.2 噴管出口欠膨脹射流問(wèn)題

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

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

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

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

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

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

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

4 結(jié) 論

以ENO格式為基礎(chǔ),通過(guò)在單元交界面處進(jìn)行三階ENO格式插值,構(gòu)造了一類(lèi)求解Lagrange坐標(biāo)系中積分形式歐拉方程的三階ENO有限體積法。數(shù)值計(jì)算結(jié)果表明,該方法具有較高的數(shù)值精度,適用于非穩(wěn)態(tài)流場(chǎng)的數(shù)值計(jì)算,并且具有較強(qiáng)的激波捕捉能力,能夠清晰地模擬出復(fù)雜流場(chǎng)的射流結(jié)構(gòu),與相同工況下實(shí)驗(yàn)結(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ò)展特性的實(shí)驗(yàn)研究[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].國(guó)防科技大學(xué)學(xué)報(bào),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].空氣動(dòng)力學(xué)學(xué)報(bào),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)排氣一維非定常流動(dòng)的基本無(wú)振蕩有限體積法的研究[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].空氣動(dòng)力學(xué)報(bào),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國(guó)標(biāo)學(xué)科代碼:13025

:A

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

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

2015-10-26;

:2016-03-31

國(guó)家自然科學(xué)基金項(xiàng)目(51475462)

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