一种三阶有限体积法及其在欠膨胀射流激波结构数值模拟中的应用*

谢 政,谢 建,李 良

(第二炮兵工程大学发射理论与技术军队重点实验室,陕西西安710025)

一种三阶有限体积法及其在欠膨胀射流激波结构数值模拟中的应用*

谢 政,谢 建,李 良

(第二炮兵工程大学发射理论与技术军队重点实验室,陕西西安710025)

以喷管出口欠膨胀射流为研究对象,在Lagrange坐标系下建立欠膨胀射流二维积分形式的流动方程。通过在单元交接面处进行三阶ENO(essentially nonoscillatory)格式插值,构造得到一种适用于求解该方程的三阶ENO有限体积法。采用该格式对一维Sod激波管算例和喷管出口欠膨胀射流进行数值计算。计算结果表明,该方法具有高精度、基本无振荡的特点,能很好地捕捉包含激波、滑移线以及三波交点等复杂流场波系结构。计算得到的波系结构中马赫盘的位置与实验结果吻合很好,相对误差小于1.1%。

激波结构;基本无振荡;有限体积法;欠膨胀射流

在枪炮武器、火箭和航空航天等工程技术领域都涉及燃气射流冲击的问题,射流波系结构的影响因素有很多[1]。为了减少燃气射流对近场设备和工程设施的冲击危害,需要了解射流流场波系结构和流场的流动状态。由于相关的射流实验费用昂贵,且通过实验不可能详尽地研究各种因素在不同水平对射流波系结构的影响。然而,采用高精度的数值方法能够有效地捕捉燃气射流场的激波波系结构,并且与相同工况下的实验结果能很好地吻合[2-3]。

近十几年来,基于非结构/结构网格的高精度算法发展迅速,主要有TVD、DGM、ENO、k-exact有限体积方法等[4-6]。从A.Harten等提出ENO格式以后,许多学者从不同的思路出发,提出了多种形式的ENO格式,如WENO、CENO[7-8]。徐文灿等[9]采用高分辨率的ENO格式对喷管流动问题进行了数值计算,得到了推力随摆角变化的规律,与实验结果基本一致,证明了采用ENO格式能够很好地捕获复杂波系,反映激波和边界层之间的相互干扰。陆霄露等[10]也采用ENO格式计算了一维非定常进排气流动问题,计算结果表明发动机的主要性能参数都和实验结果吻合很好。

为了清晰地捕捉到流场中的间断区域(激波结构),研究一种有效求解Lagrange坐标系中积分形式欧拉方程的三阶ENO有限体积法,并且采用具有精确解的激波管算例验证该方法的有效性。最后,采用该方法求解典型的欠膨胀射流的流动问题。

1 控制方程

在Lagrange坐标系中,二维轴对称非稳态可压缩气体流动的积分形式欧拉方程[11]为:

式中:Ω(t)为由边界Γ(t)包围的运动控制体,守恒变量矢量U=[ρ,M,E]T,通量矢量F=[0,p n,p u·n]T。其中,ρ为密度,u为速度矢量,M=(Mx,My)=ρu为动量,E为总能量,n=(nx,ny)为边界Γ(t)的外法线向量。理想气体的状态方程压力p=(γ-1)ρe,e为气体的质量内能。

2 数值计算方法

二维空间域划分为m×n个计算单元。Ωi+1/2,j+1/2为一个四边形的计算单元,它的顶点分别为(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)进行离散,得到轴对称Euler方程组的半离散守恒方程形式:

式中:k为计算单元4个顶点的序号,按逆时针方向排序;pk,k+1为计算单元k顶点和k+1顶点边界上的压强。由于采用有限体积方法得到的是网格平均值由可以构造高阶插值多项式p(x,y),但不能保证p(x,y)在网格边界上是连续的,因此pk,k+1不能直接从构造的插值多项式p(x,y)求得。这里需要解一个Riemann问题,通过求精确Riemann解,可以得到pk,k+1的值。而法向速度un和网格侧面积的乘积取为:

式中:uk为计算单元Ωi+1/2,j+1/2中序号为k的顶点在x方向的速度。同理,vk为y方向的速度,(xk,yk)为该顶点的坐标值。Δt时间后网格中心新速度和能量^E可由式(4)、(5)联立求得。计算单元顶点的新速度采用格点型格式直接计算,构造顶点处速度的控制体为图1中虚线边界包围的控制体在二维空间上进行三阶ENO重构,则有插值多项式:

图1 速度的控制体Fig.1 Control volume of velocity

式中:q0、qx、qy、qxx、qxy、qyy为待定系数。将边界的线段用参数方程x=x1+(x2-x1)t,y=y1+(y2-y1)t,根据文献[10],计算单元顶点运动的位置根据每个顶点的新位置求得计算单元新的面积进而得到新的密度为计算单元内流体质量。

联立式(2)～(4),在结构网格下采用高阶ENO-FVM格式,方程(1)能够表示[4]为:

求解过程忽略3阶以上高阶项的影响,方程(1)的求解转化为常微分方程的求解问题。常微分方程的求解采用具有TVD性质的三阶Runge-Kutta方法,如下式[10]所示:

式中:时间步长Δt根据式来确定,其中,r表示第r个计算步;边界中最短边的长度为计算单元Ωi+1/2,j+1/2内的声速;计算中Courant数λ取为0.5。

3 计算结果与分析

3.1 Sod激波管算例

计算域为[0,2],所采用的网格数为100,初始条件为:

采用Newmann边界条件,取CFL系数为0.5,计算推进到t=0.5 s时终止计算。此时,流场中包含一个激波、一个接触间断和一个膨胀波。图2给出了采用数值方法得到的密度分布曲线。从图2可以看出,三阶ENO有限体积法的激波分辨率很高,能准确捕捉到激波结构,将激波被抹平的厚度限制在一个网格单元尺度。表1给出了不同网格尺度下,三阶ENO有限体积法的密度误差L1范数及其收敛精度[6]。从表1可以看出,随着网格尺度的减小,计算方法体现出了网格收敛性,并且收敛的精度大致相当,计算时间以指数倍增长。

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

3.2 喷管出口欠膨胀射流问题

根据文献[1],计算了喷管出口欠膨胀射流例子,计算域为喷管子午面截得的一半区域,如图3所示。图中AB表示喷管的出口半径de,AC表示喷管的中心轴线,上游边界BD和外边界DE为静止大气条件,下游边界CE采用外推。采用无量纲格式计算,参考长度de=10 mm,参考压力p0= 101.325 k Pa,参考温度T=300 K。计算区域为[0,9]×[0,4]计算区域内布置270×160个网格,整个计算域初始时为静止大气条件,在t=0时刻射流突然喷出。计算了两种工况,工况1轻度欠膨胀,喷管出口压力比为2.06;工况2重度欠膨胀,喷管出口压力比为26.4。

表1 三阶ENO有限体积法密度误差Table 1 Density error for third-order ENO finite volume method

图3 计算区域Fig.3 Computational domain

图4 工况1下计算结果等值线图和纹影照片Fig.4 Contour maps and schlieren photograph in case 1

图5 工况2下计算结果等值线图和纹影照片Fig.5 Contour maps and schlieren photograph in case 2

图4(a)～(c)分别给出了工况1在第1 000个计算时间步时刻的马赫数、压力、密度等值线分布。从图4可以看出,射流的流场结构类似钻石状,马赫盘略有呈现,入射激波、反射激波、马赫盘等构成的波胞交替产生。在下游边界处,由于射流与静止大气之间的对流作用,流场中产生了间断面。由于Euler方程可以看作是高Reynolds数下的近似,因此边界层区域附近将出现Kelvin-Helmholtz不稳定性,越往下游不稳定性越明显。计算结果与正格式数值计算的结果[12]和实验纹影照片图4(d)反映的流场结构吻合较好。

图5(a)～(c)分别为工况2下马赫数、压力和密度等值线分布。与工况1相比,由于喷管出口压力比增大,图5中有明显的马赫盘结构,在马赫盘的边缘与入射激波和反射激波交汇形成三叉激波,出现三波交点和滑移线,马赫盘上游射流的结构较稳定,下游流场出现了明显的Kelvin-Helmholtz不稳定性。由于出口压力比增大,计算域内只能形成一个波胞,并且马赫盘的过渡区间很窄,说明本文中采用的方法具有较强的捕捉激波的能力。图5(d)为相同流动条件下的实验纹影照片[13],实验结果中第一个马赫盘距喷管出口平面的距离为4.56de。采用数值方法得到的马赫盘的位置为4.51de,文献[12]中采用正格式的结果为4.68de,两种算法相比较,采用本文方法得到的马赫盘位置更精确。与实验得到的马赫盘位置比较,本文计算结果的误差小于1.1%,且与实验纹影照片所反映的流场结构吻合较好,表明该方法具有较高的精度,能精确捕捉激波位置,并且在激波面上不会产生振荡或抹平间断现象。

4 结 论

以ENO格式为基础,通过在单元交界面处进行三阶ENO格式插值,构造了一类求解Lagrange坐标系中积分形式欧拉方程的三阶ENO有限体积法。数值计算结果表明,该方法具有较高的数值精度,适用于非稳态流场的数值计算,并且具有较强的激波捕捉能力,能够清晰地模拟出复杂流场的射流结构,与相同工况下实验结果吻合较好。

[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] 刘小军,傅德彬,牛青林,等.燃气射流冲击传热特性的数值模拟[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] 薛晓春,余永刚,张琦.双股燃气射流在充液室内扩展特性的实验研究[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].国防科技大学学报,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格式在射流数值模拟中的应用[J].空气动力学学报,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].内燃机工程,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] 朱孙科,陈二云,马大为,等.燃气自由射流的正格式数值模拟[J].空气动力学报,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国标学科代码:13025

:A

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

(责任编辑 张凌云)

2015-10-26;

:2016-03-31

国家自然科学基金项目(51475462)

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