入射激波强度对抛物形凹腔内激波聚焦影响的数值研究

陈鑫， 张强， 何立明， 王育虔， 荣康

(1.空军工程大学 航空航天工程学院, 陕西 西安 710038; 2.先进航空发动机协同创新中心, 北京 100191)



入射激波强度对抛物形凹腔内激波聚焦影响的数值研究

陈鑫1,2， 张强1， 何立明1， 王育虔1， 荣康1

(1.空军工程大学 航空航天工程学院, 陕西 西安 710038; 2.先进航空发动机协同创新中心, 北京 100191)

为得到入射激波强度对抛物形凹腔内激波反射聚焦过程的影响，采用结构自适应加密网格和波传播算法，计算马赫数分别为1.1、1.3、1.5、2.0、3.0和4.0的平面入射激波在二维抛物形凹腔内的反射和聚焦过程。计算结果表明：入射激波在凹腔壁面上的反射依次为直接马赫反射、反转马赫反射和过渡规则反射；三波点及马赫干在对称面上的碰撞是形成高能区域的原因；随激波强度增加，在凹腔壁面反射的激波由在聚焦前后都在对称面上反射变为只在聚焦后反射，再到最后不会反射；三波点在壁面上的反射和在对称面上的反射随入射激波马赫数增大，其距离和时间间隔缩小。

流体力学； 结构自适应加密网格； 波传播算法； 马赫反射； 激波聚焦

0　引言

自从1878年马赫发现激波的反射现象以来，经过一个多世纪的研究，平面激波在楔面上的反射现象研究取得了很大的进展[1-2]，尤其是近年来发现了新的激波结构很好地解释了困扰人们很久的von Neumann疑题。激波在曲面上的反射比在楔面上的反射更为复杂[3-4]，Ben-Dor对其进行了详细的总结[5]。

对平面激波在半封闭凹腔内反射现象的研究最早可追溯到1942年由Guderley得到的收敛激波的半径与时间的自相似函数关系，其结果在以后的实验中得到了验证[6]。最早的实验研究是1976年Sturtevant[7]开展的弱激波聚焦实验研究。此后，文献[8-17]对平面激波在无反应气体和可燃混合气体中聚焦起爆的现象进行了一系列的实验和数值模拟研究。Izumi等[9]总结了平面激波在二维抛物形凹腔内的激波聚焦过程中出现的3种不同类型。Shugaev等[11]主要分析了不同因素对聚焦诱导射流的影响。Gelfand等[14]对激波聚焦在可燃气体中诱导爆震进行了实验研究，发现激波聚焦诱导爆震有直接起爆和燃烧转爆轰(DDT)起爆两种不同的形式。

但这些研究主要关注于不同因素对激波聚焦效果的影响，而没有关注入射激波在凹腔壁面上反射过程中激波结构的演变。此外，在已开展的数值计算中，网格尺寸较大，或者采用了缩比模型，计算结果的激波分辨率不是很高，未能清晰地展现激波聚焦过程中入射激波及反射激波结构的演变过程。

随着计算机硬件技术的发展和各种高精度算法的出现，可以实现对间断面精确捕捉和光滑流场的高精度计算。一般来说，提高计算精度需要减小网格的尺寸，而在非定常计算中为满足稳定性条件要求网格尺寸和时间步长同步减小。所以，在多维非定常计算中计算量增长得很快。自从1982年Berger[18]提出结构自适应加密算法(SAMR)及在求解双曲方程中应用[19]以来，SAMR广泛运用于各种多维多尺度高分辨率非定常问题的计算中。SAMR可以分块分层次加密网格和缩小时间步长，能够显著地减少计算量。波传播算法(WPM)是Leveque[20]提出的满足二阶总变差衰减(TVD)条件，能用于求解一维、二维和三维非线性双曲方程的数值算法，并在Clawpack软件包中得到实现。

为研究马赫数对激波在抛物形凹腔内的反射过程的影响，采用SAMR和波传播算法[21]对不同强度的入射激波(马赫数Ma分别为1.1、1.3、1.5、2.0、3.0、4.0)在x=y2的抛物形凹腔内的反射和聚焦过程进行了计算和分析。

1　计算方法和物理模型

1.1控制方程

忽略传热和粘性作用，计算采用二维欧拉方程和理想气体模型：

(1)

式中：U=[ρ,ρu,ρv,ρ(e+V2/2)]T，V2=u2+v2；F=[ρu,ρu2+p,ρuv,ρ(e+V2/2)u+pu]T；G=[ρv,ρuv,ρv2+p,ρ(e+V2/2)v+pv]T.

1.2数值计算方法

二维波传播算法在求解Riemann问题的基础上定义了修正项，使其在空间和时间上都具有二阶精度，使用波限制函数满足TVD条件，可以写成如下形式[20]：

(2)

(3)

(4)

(5)

(6)

详细的推导可参见文献[20]。本文中使用Roe算法求解Riemann问题。

1.3数值计算方法

物理模型和网格划分如图1所示，左侧为入口，其余边界为反射边界条件，工作介质为氩气。间断面初始位置为x=-0.01 m，右侧为静止气流，压力为101 325 Pa，间断面左侧根据右侧气体性质和运动激波马赫数初始化为相应的值。为在复杂几何模型中使用矩形正交网格，采用嵌入式边界条件。用水平集函数φ表示右侧曲线边界，φ>0时为计算区域，φ<0为外部区域。在本文的计算中抛物线为无量纲参数x=cy2，其中y=Y/H，x=X/H，c=H×C=1，X、Y、C为有量纲几何参数。计算区域高度H=0.1 m，长度L=0.3 m，网格总共有5层，加密倍数为2×2×2×4. 在初始网格数为32×96，64×192和128×384的网格上进行了计算，得到的结果基本没有区别。为兼顾计算效率和计算精度，所以在本文的计算中选取初始网格数为64×192的网格。图1为马赫数Ma=1.1，t=250 μs时的网格分布，可以看出在激波所在位置网格层数多，在均匀流场区域网格层数少。

图1　计算网格和边界条件Fig.1　Computational domain and boundary condition

1.4计算方法验证

为验证数值算法和模型，在与文献[9]相同的条件(C=0.5，Ma=2.0)下进行计算。实验照片与数值计算得到的密度梯度幅值云图的对比如表1所示。从图中可以看出，数值计算结果与实验照片十分吻合，捕捉到了流场中的细微结构，如滑移线和马赫反射等复杂激波结构，证明本文的计算结果有较高的可信度。

表1文献[9]中的实验照片与计算结果对比

Tab.1Comparison of experimental photos and numerical results

2　结果分析

以平面激波在楔面上的反射结构为参照，分析了不同强度(Ma分别为1.1、1.3、1.5、2.0、3.0、4.0)的入射激波在x=y2的二维抛物形凹腔内反射过程中激波结构的变化。入射激波在凹腔内的反射过程可以分为入射激波在壁面上的反射、三波点在壁面上的反射和三波点在对称面上的反射3个阶段。在以下计算结果的分析中入射激波、反射激波、马赫干、三波点和滑移线分别用i、r、m、T和s表示。为便于观察和分析，若没有特别申明，在下文中采用密度梯度的幅值云图给出聚焦过程中不同条件和时刻的激波结构。

2.1入射激波在壁面上的反射

不同强度的入射激波在凹腔壁面上的反射均为马赫反射。在Ma=1.1时，由于激波非常弱，在计算精度不够时误以为是正常反射(RR)。进一步加密网格(2×2×4×4)后可以观察到非常短的马赫干，但没有发现滑移线，如图2(a)所示。在其余情况中，激波反射都为单马赫反射(SMR)。根据文献[5]中的分类标准，可根据马赫干长度的变化对SMR进行分类。在图2(b)～图2(f)中，入射激波刚进入凹腔时，三波点远离壁面，马赫干长度增加，滑移线s向壁面的方向弯曲，称之为直接马赫反射(DiMR)。随后，由于壁面曲率增大，马赫干长度减小，三波点向壁面运动，称之为反转马赫反射(InMR)。在DiMR向InMR转变过程中，会经过马赫干长度不变的位置，入射激波在此点的反射称之为固定马赫反射(StMR)。

此外，在入射激波马赫数Ma≤2.0时滑移线s向马赫干弯曲，为稳定结构，如图2(b)～图2(d)所示；在入射激波马赫数Ma分别为3.0、4.0时滑移线尾端卷曲，为不稳定结构，如图2(e)～图2(f)所示。随着入射激波i的强度增大，反射激波r的强度也增大，入射激波与反射激波的夹角θ也随之增大，如图2(b)所示。随入射激波向前运动，由于反射角度增加，θ也会增大。

图2　入射激波在凹腔壁面上反射过程Fig.2　Reflection of incident shock wave from cavity wall

2.2三波点在壁面的反射

入射激波在壁面上的反射由最初的DiMR转变为InMR，最后三波点在壁面反射，形成过渡规则反射(TRR)。在入射激波较弱时，TRR结构清晰可见，如图3所示。在Ma=1.1时，三波点在壁面反射过后形成类似TRR的结构(见图3(d))，但没有观察到滑移线。这种结构在文献[9]中作为一个未解决的问题有所提及。根据计算结果判断，这应该是一个不完整的TRR结构。在实验中由于分辨率限制和粘性作用难以观察到微弱的滑移线。而在计算中，数值耗散也使得微弱的滑移线难以持续存在。

图3　Ma分别为1.1、1.3和1.5时三波点在壁面上的反射Fig.3　Reflection of triple points from cavity wall with incident Mach numbers of 1.1, 1.3 and 1.5

在Ma=1.3时，三波点在反射前贴于壁面，如图3(b)所示。然后反射形成RR和SMR组合的激波结构，即TRR，如图3(e)所示。滑移线s的转折点大致为三波点在壁面碰撞的位置。在Ma=1.5时，TRR结构中的马赫干很短，通过局部放大才能辨认出反射为TRR结构，如图3(f)所示。从图3(a)～图3(c)中可以看出，随着入射激波马赫数增大，三波点在壁面反射的位置下移(向凹腔顶点靠拢)，TRR结构空间尺度缩小，变得难以分辨。

入射激波强度进一步增大，三波点在壁面反射位置的高度继续降低，最后在Ma=4.0时观察不到三波点在壁面的反射(见图4)。对比图2和图3、图4中可以发现，入射激波与反射激波的夹角θ随激波运动和马赫数增大而增大。在入射马赫数Ma=3.0和Ma=4.0时，滑移线的不稳定进一步加剧，形成Kelvin-Helmholtz不稳定结构(KH结构)，如图4(b)和图5(c)所示。

图4　Ma分别为2.0、3.0和4.0时三波点在壁面上的反射Fig.4　Reflection of triple points from cavity wall with incident Mach numbers of 2.0, 3.0 and 4.0

图5　Ma分别为1.1和1.3时三波点在对称面上的反射Fig.5　Reflection of triple points from symmetric plane with incident Mach numbers of 1.1 and 1.3

2.3三波点及马赫干在对称面上的反射

三波点在壁面反射后，形成TRR激波结构(见图5(a)和图5(e))，在激波强度较弱时很容易捕捉到波系结构的发展，如图5所示。TRR结构中的入射激波i在壁面上的反射激波r2，与反射激波r和TRR中的马赫干m形成一个新的三波结构，同时，弧形的反射激波在对称面上规则反射，形成两个反射点R1和R2，如图5(b)和图5(f)所示。马赫干m继续向对称面运动，然后在对称面上反射，如图5(c)和图5(g)所示。由于马赫干m向对称面弯曲，马赫干的反射先于三波点的反射。马赫干及三波点在对称面上的反射形成了一道弧形的反射激波r3，并向远离凹腔顶点的方向运动，如图5(d)和图5(h)所示。由于θ较小，反射激波的形状为凹形曲线，所以反射激波r在聚焦前后都在对称面上反射，形成R1和R2两个反射点。从壁面反射的激波在对称面上的规则反射和聚焦时形成的弧形激波在对称面上的马赫反射组成的激波形态类似于TRR结构，如图5(d)和图5(h)所示。反射激波r在对称面上的两个相交点中靠近凹腔顶点的一个交点R1与聚焦形成的反射激波r3和m合并成一个三波点并向远离顶点的方向运动，在Ma=1.3时可以观察到明显的滑移线，见图5(h)。

在入射激波马赫数Ma=1.5和Ma=2.0时(见图6)，激波聚焦过程与图4所示的激波聚焦过程类似同，但有其不同之处。随着Ma增大，TRR结构在时间和空间上的尺度缩小，从而使得入射激波在凹腔顶点反射后形成的马赫干m在对称面碰撞时长度减小，如图6(c)和图6(g)所示。这种尺度上的缩小使得聚焦的强度增大。此外，反射激波r只在聚焦后在对称面上反射(聚焦前没有在对称面上反射)，如图6(d)和图6(h)所示。图6(h)中的滑移线在对称面上有一个转折点，说明三波点在对称面上发生了反射。正是因为这次聚焦后形成的三波点的在对称面上的反射引起反射激波r在对称面上的再一次反射。

图6　Ma分别为1.5和2.0时三波点在对称面上的反射Fig.6　Reflection of triple points from symmetric plane with incident Mach numbers of 1.5 and 2.0

从Ma=2.0时，激波聚焦前后凹腔顶点附近的温度和压力云图分布可以看出(见图7)，三波点以及马赫干在对称面上的反射后，形成了高温高压区域.

图7　Ma为2.0时聚焦前后顶点附近温度和压力分布Fig.7　Temperature and pressure distribution at the vicinity of cavity apex before and after focusing

图8　Ma分别为3.0和4.0时三波点在对称面上的反射Fig.8　Reflection of triple points from symmetric plane with incident Mach numbers of 3.0 and 4.0

在入射激波马赫数Ma=3.0和Ma=4.0时，激波强度增加，三波点在壁面和对称面上的反射几乎同时进行，如图8所示。由于激波强度大，运动速度快，两个反射发生的间隔也很短。这种激波在狭小空间的快速多次碰撞比三波点的一次碰撞形成的高能区能量密度更大，反射激波也更强。在入射激波较强时，反射激波r在聚焦前后不会在对称面上相交。在Ma=3.0时，聚焦反射激波远离凹腔顶点运动过程中马赫干m的长度先增加后减小，而后又增加。这与Ma=2.0的情况类似，但马赫干在此过程中没有消失，所以三波点没有碰撞，反射激波r也不会相交。

3　结论

采用结构自适应加密网格和波传播算法，对不同马赫数的平面入射激波在型面为x=y2的凹腔内反射过程进行了数值计算，定性描述和分析了在反射和聚焦过程中激波结构的运动和发展。得到了以下结论：

1)不同强度的入射激波在半封闭凹腔壁面上的反射过程都包括DiMR、InMR及TRR 3个过程。

2)三波点以及马赫干在对称面上的碰撞是激波聚焦形成高温高压区的根本原因。

3)随着入射激波变强，反射激波由在聚焦前后都在对称面上反射变为只在聚焦后反射，到最后不会反射。反射激波与入射激波的夹角和反射激波的形状是决定反射激波能否再次反射的主要因素。

4)随着入射激波强度的增大，三波点在壁面的反射点和对称面上的反射点之间的距离变小，时间间隔缩短，激波聚焦的效果进一步加强。

References)

[1]杨旸，姜宗林，胡宗民. 激波反射现象的研究进展[J]. 力学进展, 2012, 42:142-160.

YANG Yang, JIANG Zong-lin, HU Zong-min. Advances in shock wave reflection phenomena[J]. Advances in Mechanics, 2012, 42:142-160. (in Chinese)

[2]Ben-Dor G. A state-of-the-knowledge review on pseudo-steady shock-wave reflections and their translation criteria[J]. Shock Waves, 2006, 15(3/4):277-294.

[3]Sebastien Gruber. Weak shock wave reflections from concave curved surfaces[G]. Johannesburg, South Africa: University of the Witwatersrand, 2012.

[4]Skews B W, Kleine H, Bode C, et al. Shock wave reflection from curved surfaces[C]∥22nd International Congress on Theoretical and Applied Mechanics. Adelaide, Australia: IUTAM, 2008.

[5]Ben-Dor G. Shock wave reflection phenomena[M]. Berlin: Springer-Verlag, 2007:247-306.

[6]Veronica Eliasson. On focusing of shock waves[R]. Stockholm, Sweden: Royal Institute of Technology KTH Mechanics,2007.

[7]Sturtevant. The focusing of weak shock waves [J]. Fluid Mech, 1976, 73:651-671.

[8]Kishige H, Teshima K, Nishida M. Focusing of shock waves reflected from and axisymmetric parabolic wall[J]. Shock Waves, 1991, 1:341-345.

[9]Izumi K, Aso S, Nishida M. Experimental and computational studies focusing process of shock waves reflected from parabolic reflectors[J]. Shock Waves, 1994, 3: 213-222.

[10]Liang S M, Wu L N, Hsu R L. Numerical investigation of axisymmetric shock wave focusing over paraboloidal reflectors[J]. Shock Waves, 1999, 9:367-379.

[11]Shugaev F V, Serov A O, Shtemenko L S. Formation of a jet and vortices behind a shock wave reflected from a concave body[J]. Shock Waves, 1999, 9:31-35.

[12]Grönig H. Past, present and future of shock focusing research[C]∥Proceedings of International Workship on Shock Wave Focusing. Sendai, Japan: Institute of High Speed Mechanics, 1989: 1-37.

[13]Chan C K, Lau D, Thibault P A, et al. Ignition and detonation initiation by shock focusing[J]. AIP Conference Proceedings, 1990,208: 161-166.

[14]Gelfand B E, Khomik S V, Bartenev A M, et al. Detonation and deflagration initiation at the focusing of shock waves in combustible gaseous mixture[J]. Shock Waves, 2000, 10:197-204.

[15]王春, 司徒明, 韩肇元. 用于爆震引燃的激波聚焦无反应流场数值模拟[J]. 推进技术, 2003, 24:156-159.

WANG Chun, SITU Ming, HAN Zhao-yuan. Numerical investigations on cold flowfields of shock focusing for ignition of pulse detonation[J]. Journal of Propulsion Technology, 2003, 24:156-159.(in Chinese)

[16]董刚, 唐敖, 叶经方, 等. 激波聚焦诱导点火和爆轰的数值研究[J]. 爆炸与冲击, 2005, 25(5): 437-444.

DONG Gang, TANG Ao, YE Jing-fang, et al. Numerical studies on initiation and detonation induced by shock wave focusing[J]. Explosion and Shockwaves, 2005, 25(5):437-444. (in Chinese)

[17]滕宏辉, 张德良, 李辉煌, 等. 用环形激波聚焦实现爆轰波直接起爆的数值模拟[J]. 爆炸与冲击, 2005, 25:512-518.

TENG Hong-hui, ZHANG De-liang, LI Hui-huang, et al. Numerical investigation of detonation direct initiation induced by toroidal shock wave focusing[J]. Explosion and Shockwaves, 2005, 25:512-518. (in Chinese)

[18]Berger M. Adaptive mesh refinement for hyperbolic differential equations[D]. CA: Stanford University, 1982.

[19]Berger M, Colella P. Local adaptive mesh refinement for shock hydrodynamics[J]. Journal of Computational Physics, 1988, 82:64-84.

[20]Leveque R J. Finite volume method for hyperbolic problem[M]. Oxford City: Cambridge University Press, 2002.

[21]Ralf Deiterding. Block structured adaptive mesh refinement, theory, implementation and application[J]. ESAIM Proceeding, 2011, 34:97-150.

Numerical Research on the Effect of Incident Shock Strength on Focusing of Shock Wave in a Parabolic Cavity

CHEN Xin1,2, ZHANG Qiang1, HE Li-ming1, WANG Yu-qian1, RONG Kang1

(1.Aeronautics and Astronautics Engineering College, Air Force Engineering University, Xi’an 710038, Shaanxi, China; 2.Collaborative Innovation Center of Advanced Aero-Engine,Beijing 100191, China)

In order to obtain the influence of incident shock strength on the focusing of shock wave in a parabolic cavity, the adaptive refined mesh and wave propagation method are used to simulate the reflection and focusing of planar shock wave in a 2D parabolic cavity with incident shock Mach numbers of 1.1, 1.3, 1.5, 2.0, 3.0 and 4.0. Numerical results indicate that the reflection of incident shock from the cavity wall is direct Mach reflection, invert Mach reflection and transition regular reflection in order. High energy zone is resulted from the collision of triple points as well as Mach stems. The shock wave reflecting from the cavity wall before and after focusing only reflects from symmetry plane after focusing and finally does not reflect as the incident shock strength increases. The distance and interval between triple points reflecting from wall and symmetry plans decrease with the increase in the Mach number of incident shock wave.

fluid mechanics; structure adaptive mesh refinement; wave propagation method; mach reflection; shock wave focusing

2014-04-14

国家自然科学基金青年科学基金项目(51106178)

陈鑫(1976—), 男, 副教授。E-mail:chenxin7605@aliyun.com

V231.2+2

A

1000-1093(2015)04-0631-06

10.3969/j.issn.1000-1093.2015.04.009