-
Based on Navier-Stokes equations, combining the fifth-order weighted essentially non-oscillatory scheme with the adaptive structured grid refinement technique, the interactions between converging shock and annular SF6 layers with different initial perturbation amplitudes and thickness are numerically investigated. The evolution mechanism of shock structure and interface are revealed in detail, and the variations of the circulation, mixing rate and turbulent kinetic energy are quantitatively analyzed. The dynamic mode decomposition method is used to analyze the dynamic characteristics of the vorticity. The results show that in the case with large initial perturbation amplitude, the transmitted shock wave forms Mach reflection structures both inside and outside of the inner interface of SF6 layer, and multiple shock focusing phenomena occur in the center. After the transmitted shock wave penetrates the outer interface, the circulation increases faster, and the “spike” and “bubble” structure on inner interface develop faster, so that the amplitude of the inner and outer interfaces and the gas mixing rate increase. As for the case with larger thickness of the gas layer, the phase of the transmitted shock wave changes inside the layer, which forms “bubble” at the crest of the inner interface and “spike” at the trough. When the thickness of the gas layer decreases, the crest of the inner interface does not move inside after being impacted, and “spike” and “bubble” are generated in the late stage. The dynamic modes show that the main structure of vorticity and the exchange of positive and negative vorticity on the main structure are determined by the modes with weak growth and low frequency, but the modes with weak growth and high frequency determine rapid exchange of positive and negative vorticity at the interface in the cases with weak coupling effect.
1. 引 言
汇聚激波与界面相互作用过程包含界面演化、激波与涡相互作用以及湍流混合等物理现象, 其中还包含Richtmyer-Meshkov (RM)不稳定性、Kelvin-Helmholtz (KH)不稳定性和Rayleigh-Taylor (RT)不稳定性. 而RM不稳定现象广泛存在于超音速燃烧[1]、水下爆炸[2]、惯性约束核聚变(inertial confinement fusion, ICF)[3]等领域中. 因此, 研究汇聚激波诱导界面不稳定性中界面和激波结构的演变、流动转捩和后期湍流混合的机理及其影响因素, 可以进一步促进或抑制过程中的不稳定性. 平面激波诱导界面RM不稳定性的研究较为广泛, 理论方面主要研究界面扰动的演化规律[4], 表征不同阶段的扰动增长规律的模型[5]. 实验和数值模拟方面则主要研究不同条件(如激波强度[6]、Atwood数[7]以及气柱形状[8,9]等)下激波冲击气泡或气柱的界面演变过程和激波结构.
近几年汇聚界面不稳定性的研究受到广泛关注. 理论研究方面, Mikaelian[10]建立了球形的分层壳体的积分模型, 并扩展到N层不同的流体的积分模型, 得到了汇聚RT和RM不稳定性中湍流混合层的宽度, 且模型成立条件为幅长比小于0.1. Lombardini和Pullin[11]将Mikaelian的理论推广到三维情况, 研究汇聚激波冲击柱体模型的扰动界面的扰动振幅变化规律, 并提出三维扰动增长的线性模型. 在实验方面, Si等[12]对竖直同轴汇聚激波管进行改进, 研究了正多边形SF6重气柱界面在汇聚激波作用下的演化过程, 发现在前期界面上斜压涡量是界面演化的主要原因, 其中正三角形气柱涡量最高, 且可压缩性在整个作用过程都起到重要作用. Ding等[13]实验研究了汇聚激波冲击带有正弦扰动的单层SF6重气柱界面, 透射激波速度在离中心较远时维持在一个定值, 说明在远离中心处收敛效应可以忽略. 近几年汇聚激波研究着重于更复杂的双层界面, 其中不仅包含多次激波透射与反射, 还需要考虑到内外界面耦合效应. Ding等[14]使用竖直同轴汇聚激波管研究了外界面带有正弦初始扰动界面重气柱的汇聚RM不稳定性, 实验还展示外界面扰动幅值随时间变化过程, 同时引用Mikaelian[15]提出的耦合角来对界面耦合效应强度进行量化. Sun等[16]采用内界面有正弦初始扰动的界面进行实验, 还展示了内界面扰动幅值和气泡尖钉结构幅值随时间变化过程, 并根据扰动幅值变化将整个过程分为三个阶段. Li等[17]则将重气柱换成轻气柱进行研究, 与重质气体相比, 初始激波接触外界面后产生透射激波和稀疏波, 透射激波也产生了与外表面一致的扰动, 但是扰动在不断衰减, 到达内表面时已经十分微弱, 因此对内表面影响很小. 而数值研究方面, 计算方法上徐建于和黄生洪[18]则基于Harten Lax van Leer (HLL)黎曼求解器的SPH(smoothed particle hydrodynamics)算法对二维汇聚激波冲击四边形轻/重气界面的RM不稳定性问题进行数值模拟, 与已有实验结果对比证明了该方法的可靠性. 梁煜等[19]对比汇聚激波及平面激波冲击下SF6球形气泡演化规律, 发现汇聚激波条件产生的涡量幅度大于平面激波条件, 界面混合也更快, 当反射激波再次作用于界面后产生了正负涡量. Zhou等[20]研究了双模界面的演化过程, 发现由于耦合效应的影响, 第一模态的幅值增长取决第二模态与第一模态的波数比, 为偶数和奇数时分别对扰动增长产生抑制和促进作用. 而Tang等[21]采用不同气体组研究了Atwood数对RM不稳定性的影响. 何慧琴等[22]则研究偏心对汇聚激波冲击不同形状的二维气柱界面RM不稳定性的影响, 偏心导致界面环量分布不均, 从而使得扰动结构发生扭曲, 且尖钉两侧的环量相差较多, 结构扭曲更严重. Fu等[23]数值模拟研究了三维圆柱和球体汇聚RM不稳定性过程中的能量传输过程, 计算了两种界面RM不稳定性湍动能运输方程中各项(湍流扩散项、黏性扩散项、外力做功项、耗散项和可压缩项)沿径向的分布.
先前汇聚激波的研究主要针对不同初始条件下的界面演变过程, 而对于汇聚RM不稳定性的界面形态、尤其激波结构的演化过程等细节仍然不明确. 而且, 目前对后期界面上由于RM不稳定性出现的涡量分布以及湍流混合的机理研究也较少. 因此, 本文主要研究汇聚激波冲击具有正弦扰动双层SF6重气柱界面过程, 并改变初始扰动幅值与气层厚度, 分析界面形态与激波结构演化过程, 并对界面扰动幅值、混合率、环量及湍动能进行定量分析. 此外, 利用动模态分解(dynamic mode decomposition, DMD)方法对涡量的演化过程进行分析, 以揭示涡量输运诱导的混合机理.
2. 数值方法与验证
基于有限体积法和可压缩多组分Navier-Stokes方程[24], 结合5阶WENO (weighted essentially non-oscillatory)格式[25]以及结构化自适应网格加密技术, 对汇聚激波冲击双层重气柱界面过程进行数值模拟, 为确保对流场中复杂激波与涡结构的准确描述, 时间项采用3阶精度的Runge-Kutta法[26]. 为了减少计算域边界反射波对流场的干扰, 对计算域四周边界均采用无反射边界条件.
为验证本文数值方法的准确性, 图1是汇聚激波与双层圆形界面作用过程计算结果与Ding等[14]实验结果的对比, 计算与实验中的双层圆形重气柱模型具有相同的尺寸、气体参数及边界条件. 由图1可见, 数值结果清晰描述了汇聚激波与双层圆形界面作用过程中外界面(OI)、内界面(II)与激波(shock)位置的演变过程, 其中内外界面在受到汇聚激波冲击后开始向中心处移动, 而透射激波在中心处的汇聚反射则使内外界面减速甚至反向运动, 此过程中它们的运动轨迹与实验结果相吻合, 虽然后期界面向外的移动速度与实验值相比较小, 但误差仍可接受, 这表明本文数值方法准确可靠.
此外, 为验证本数值模拟网格无关性, 采取了1500×1500, 2000×2000和2500×2500三种不同网格密度进行模拟, 图2为此三种网格密度下t = 0.05 ms时气体密度与初始静止空气密度比ρ/ρair沿径向分布情况, 在保证准确性的同时为减少计算时间, 本文选用2000×2000这一网格密度进行数值模拟.
为研究外界面初始扰动幅值与气层厚度的影响, 设置五种不同工况进行模拟, 几何参数如表1所列, 其中R0表示外界面平均半径, r0为内界面半径, a0为初始扰动幅值, n表示初始波数, λ为初始扰动波长. 计算域边长取L = 80 mm. 初始激波马赫数为1.22. 环形气层内为SF6, 其余为空气, 比热比γ、摩尔质量M和密度ρ这三个气体参数如表2所列. 气柱内外压力均设为101.49 kPa, 初始温度均为288 K, 初始时刻气体均静止. 图3为计算模型示意图, 其中OI为外界面, II为内界面, is为初始汇聚激波, 且外界面初始扰动表示为
$ R = {R_0} + {{a}_0}\cos (n\theta ) $ .Case R0/mm r0/mm a0/mm n λ/mm a0/λ Case 1 20 10 1 6 20.94 0.048 Case 2 20 10 0.5 6 20.94 0.024 Case 3 20 10 2 6 20.94 0.096 Case 4 20 5 1 6 20.94 0.048 Case 5 20 15 1 6 20.94 0.048 Gas γ M/(g·mol–1) ρ/(kg·m–3) Air 1.399 28.967 1.23 SF6 1.103 128.491 5.45 3. 结果与讨论
3.1 界面及激波结构演变
图4为case 1中界面与激波结构的演变过程. 初始汇聚激波is穿过OI, 产生向内传播的透射激波(transmitted shock, ts1)与向外传播的反射激波(reflected shock, rs1). 由于SF6气体的声阻抗大于空气的声阻抗, 激波在SF6气体中传播速度比在空气中慢, 冲击OI波峰的激波比冲击波谷的激波更早开始减速, 导致OI的正弦扰动传递给ts1. OI的扰动幅值在波峰被is冲击后开始减小, 当is完全穿透后, OI波谷在is的冲击加速和斜压涡量的作用下开始快速向中心收缩, 导致扰动幅值快速增加. t = 0.06 ms时, ts1以近似OI的正弦扰动形态开始冲击II, ts1的波谷首先与II作用生成向外传播的反射稀疏波(reflected rarefaction wave, rrw1)和向内传播的透射激波ts2. 在ts1完全穿透II后, 生成向外传播的反射激波rs2与透射激波ts3. ts2与ts3向内传播过程中相互作用形成由马赫杆(Mach stem, m)和三波点(triple point, T)组成的马赫反射结构(t = 0.073 ms). 在向中心汇聚的过程中, ts2不断缩短, m长度几乎不变, 而ts3不断增长且与相邻ts3相交, 同时向外移动的相邻rs2也相交并生成rs3, 如t = 0.08 ms所示. 之后m和ts2向中心移动, 在t = 0.083 ms时发生第一次激波聚焦(shock focusing, SF1), 生成向外传播的二次反射激波(second reflected shock, srs). 反射激波rs3向外移动的同时也向两侧传播, 相交产生rs4. 当t = 0.096 ms时srs开始冲击II, 生成向外的二次透射激波(second transmitted shock, sts)和向中心移动的反射激波rs5, rs5向内汇聚发生二次聚焦SF2, 并生成向外移动的三次反射激波(third reflected shock, trs). II受到srs冲击发展形成了“尖钉”(spike)与“气泡”(bubble)结构, 虽然trs在“尖钉”处首先发生第三次透射产生三次透射激波(third transmitted shock, tts), 但由于激波在SF6中传播速度较慢, 因此, trs主要从“气泡”处出现并向外传播冲击OI. 之后sts冲击OI波谷, 产生向外的透射激波ts4和向内的反射稀疏波rrw2. 在t = 0.152 ms时sts穿透OI波峰, 而tts则在向外传播过程中逐渐减弱, 并在与OI上“气泡”作用后转变为微弱的膨胀波.
界面演变受涡量分布影响, 图5选取五个时刻流场的涡量分布情况对界面发展机理进行说明. II在ts1穿透后形成类似正弦曲面, 其涡量方向与OI一致(t = 0.066 ms), 而rs2在II上产生的涡量与OI对应位置涡量方向相反(t = 0.07 ms), 因此OI波峰位置对应的II波峰反而开始向内凹陷. 之后rs3诱导的涡量与rs2生成的方向相反, II上产生不稳定性, 但对II发展影响较小. 在srs冲击II后, II上不稳定性增强, 波峰形成的凹陷向内移动形成了“尖钉”结构, 波谷则向外移动形成了“气泡”结构, 且由于“尖钉”处激波与界面夹角较大, 累积了更多的同向涡量(t = 0.103 ms), 导致“尖钉”发展更快, 向内移动带动界面向内收缩, 其头部两侧形成对称的涡结构, 而“气泡”形态变化不大. 之后, “尖钉”头部向内的移动速度逐渐降低, 头部两侧涡结构开始发展, 相邻两个涡结构开始互相影响, 而“气泡”在受到rrw2冲击后产生同向涡量(t = 0.155 ms), 开始快速向外发展. OI向内收缩的同时幅值快速增大, 受到sts冲击后生成反向涡量(t = 0.135 ms), 此时OI波谷出现向外的“尖钉”, 受到tts冲击后“尖钉”开始快速发展, OI扰动幅值减小最终反相.
由于cases 2—5在ts1冲击II前的界面与激波结构演变过程类似, 这一阶段示意图省略. 图6表示case 2界面与激波结构演变过程, 由于OI初始正弦扰动幅值变小, ts1与II作用后在II上的涡量累积较少, II波峰向内凹陷不明显, 直至受到srs冲击后II上“尖钉”和“气泡”才开始缓慢发展. 图7表示case 3界面与激波结构演化过程, 由于初始扰动幅值变大, ts1为具有明显正弦扰动的形态, 且在ts1的波峰位置形成由透射反射激波(transmitted reflected shock, TRS)、ts1、马赫杆m、三波点T组成的马赫反射结构(t = 0.057 ms), 其中TRS的移动方向垂直于激波移动方向. 此外, ts1与II相互作用之后, II上生成了明显的“尖钉”“气泡”结构, 使得srs冲击“尖钉”与“气泡”的时间间隔较大, 产生了两道相位相反的反射激波rs5和rs6. rs5和rs6向中心汇聚分别发生第二次与第三次激波聚焦, 先后生成向外传播的trs和四次反射激波(forth reflected shock, frs), 在穿透II后分别形成tts与四次透射激波(forth transmitted shock, fts), 并在t = 0.134 ms时冲击OI的波谷. 界面形态与case 1相比, ts1与II夹角更大, 累积更多的斜压涡量, 且反射激波产生的不稳定性对界面影响较小, 因此II上“尖钉”更早出现. 之后srs在“尖钉”产生更多的同向涡量, 头部两侧涡发展较快卷吸周围空气, 在“气泡”内发展更充分. “气泡”受到相邻TRS相交诱导其头部处产生的高压区, 使其头部的增长受到抑制, 而在接触rrw2后开始快速发展.
从cases 1—3可以看出, 初始扰动幅值较大容易使激波在重气柱内生成马赫反射结构, 激波与II夹角更大并产生更多的涡量. 多次反射激波产生的正负交替的涡影响越小, II上更早地发展出“尖钉”, 受到srs冲击后生成两个向内的反射激波, 进而生成tts和fts. 此外, 幅值越大OI波谷越靠近II, 使“气泡”更早受到rrw2冲击并向外发展.
图8为case 4界面与激波结构演变过程. 由于II半径减小, ts1在重气柱内传播时间增加, 向内传播的过程中ts1波峰位置产生了由ts1, m和TRS组成的马赫反射结构. 当t = 0.087 ms时, ts1和m一起冲击II, 生成由ts2组成的六边形透射激波, 此时II也转变为六边形形状, 同时ts2在II波峰位置内侧也生成了马赫反射结构. 在t = 0.096 ms时发生激波聚焦生成srs, 同时II上相邻TRS相交改变了rs2的形状. 由于ts1在II上产生的涡量与OI对应位置涡量方向相反, II波谷中点向内凹陷而波峰向外凸出, 并在srs冲击后逐渐分别发展为“尖钉”和“气泡”. “尖钉”两侧相邻的涡结构十分接近, 受到空间限制无法向中心移动, 因此涡结构在“气泡”中互相混合快速发展, 而“气泡”由于没有OI的限制, 在sts和tts穿透后快速向外发展. 从图9可以看出, case 5的激波结构演变与case 1相似, 而界面形态有明显不同. II受到冲击后在多次反射激波作用下累积了正负交替的小涡, 同时II由于较强的耦合效应随OI一起向内收缩, 因此波峰没有向内凹陷. 当tts冲击OI后, OI上逐渐生成“尖钉”, II随着向外发展形成“气泡”, 而“尖钉”逐渐生成但几乎不向内发展.
对比cases 1, 4, 5可以看出, 气层厚度对整个过程的影响主要表现在两个方面. 一是激波在重气柱内的运动时间随着气层厚度增大而增加, 重气柱内马赫反射结构有更长的发展时间, 使得II受到冲击后与case 1相位相差π, 从而对II“尖钉”的生成产生影响. 二是气层厚度减小使耦合效应增强, 前期II与OI一同向内快速收缩, 后期II“气泡”向外发展受到OI限制.
3.2 径向统计量分析
图10为case 1条件下压强沿径向分布图. 在第一次激波聚焦前(t = 0.08 ms)压力峰值出现在透射激波处, 第一次激波聚焦后(t = 0.09 ms)中心区域压力突增, 峰值出现在srs处. 当t = 0.108 ms时发生第二次集激波聚焦, 中心处压力达到最大值. 而tts穿透II后, II内的压力都较高. 当t = 0.2 ms时可以看到中心处压强较之前有所降低, 且沿径向变化较小, 说明流场后期压力影响较低.
为定量地描述内外界面形态变化, 将计算域转换为直角坐标, 以case 1某时刻SF6组分分布为例, 结果如图11所示. 首先计算了不同情况下内外界面幅值随时间变化过程, 内外界面幅值分别定义为
$$ {{{a}}_{\rm i}} = ({r_{\rm o}}-{r_{\rm i}})/2 , $$ (1) $$ {{{a}}_{\rm o}} = ({R_{\rm o}}-{R_{\rm i}})/2 . $$ (2) II幅值随时间的变化如图12(a)所示. 对比cases 1—3可以看到, ai在ts1完全穿透II后先减小再缓慢增大, 而且初始幅值a0越大ai则减小越多. 这是因为ts1波谷比波峰更早冲击II, 因此II上形成波谷并向内移动, 在ts1波峰冲击II前ai达到最小值. ts1穿透II后, ai开始逐渐增大, II凹陷形成后增长为正值, 这一阶段的扰动也随着a0的增大而更加明显. 在srs冲击II后ai开始进入到线性增长阶段, 且a0越大, 这一阶段ai的增长速度越大. 在线性阶段结束后, a0小的情况更早进入非线性阶段. 在case 4中, 由于II相位与之前相差π, ai在ts1穿过II后为正, 在srs冲击II后开始快速增长. 但由于“尖钉”距离中心更近, 向内发展受到限制, 因此ai增长速度较低. 对于case 5, 由于ts1穿过II后, 波峰没有向内凹陷, 而II由于强耦合效应和OI波谷一起向中心移动, 因此ai此时为负值. 在OI“气泡”生成后, II开始发展, ai开始逐渐增加.
图12(b)给出了OI幅值随时间的变化, is冲击OI波峰后ao先快速减小, 当is完全穿透OI后ao再逐渐增大. 初期阶段ao在界面上斜压涡量作用下线性增长, 在向内收缩的过程中增长速度缓慢增大, 当rrw1传播到OI后ao开始快速增大. 但到达在峰值前由于RT效应有一个减速阶段, 在sts冲击OI波谷后ao开始减小. 当tts冲击OI波谷后, ao减小速度加快, 并随着波谷“尖钉”向外发展最终减小到负值, OI上发生反相. 在此过程中, 当sts冲击并穿透OI波峰时, 幅值减小速度有所放缓. 对比cases 1—3可以看到, 初始扰动幅值较大的OI在线性阶段有更大的增长速度, 能够更快到达快速增长阶段, 且非线性阶段增长幅度更大. 而对比case 1, case 4和case 5可以看出, 气层厚度主要影响了ao从线性增长到非线性增长的时间, 而对ao最大值影响较小. 在case 5中耦合效应增强, 非线性阶段ao快速增加.
激波与界面作用生成的涡强度可以通过界面上环量的变化来表征, 而由于对称性, 界面总环量始终为零, 因此定义环量绝对值如下:
$$ \left| \varGamma \right| = \int_D {\left| \omega \right|} {\rm d}s = \sum\nolimits_D {\left| {{\omega _i}} \right|{A_i}} , $$ (3) 其中区域D包含了所有含有SF6气体的流体微元, 即SF6质量分数 f > 0, 而ωi和Ai分别表示区域D中第i个网格内的涡量与面积, 结果如图13(a)所示, 可以看到几种情况下环量变化过程比较相似. 在is冲击OI时
$ \left|{ \varGamma }\right| $ 突增, 在ts1向II移动的过程中加速增大, 增速与初始扰动幅值正相关. 在ts1冲击II时,$ \left|{ \varGamma }\right| $ 会有一定降低. 当ts1完全穿透后$ \left|{ \varGamma }\right| $ 会小幅增大, 且在rs2和rs4生成时,$ \left|{ \varGamma }\right| $ 增速较大, 这与反射激波与II作用交替生成的同向与反向涡量有关. 而之后sts, tts冲击OI时,$ \left|{ \varGamma }\right| $ 均会突然减小再增加, 而增速相较之前有所放缓, 原因是激波强度较低且与“气泡”夹角较小. 而气层厚度的影响不仅体现在ts1冲击II前$ \left|{ \varGamma }\right| $ 增长的时间, 而且还通过耦合效应的强弱影响环量的增长. 在case 4中,$ \left|{ \varGamma }\right| $ 在ts1冲击II前增长时间较长且增速较快, 且耦合效应弱使得II形态充分发展, 激波冲击在界面上产生的斜压涡量更多. 但由于rrw2在重气柱内运动时间较长与“气泡”作用时强度降低, 后期流场涡量在黏性作用下逐渐耗散, 涡量强度有明显下降. 而在case 5中, sts冲击II后耦合效应使得II形态变化小, 导致界面上斜压涡量少, 环量增长缓慢, tts冲击OI后逐渐振荡减小.图13(b)表示不同情况混合率随时间变化图, 混合率定义为
$$\begin{split} \xi &\; = 1 - \dfrac{{\displaystyle\int_D {f(x,y,t)} \rho (x,y,t){\rm d}s}}{{\displaystyle\int_D {\rho (x,y,t){\rm d}s} }} \\ &\;= 1 - \frac{{\displaystyle\sum\nolimits_D {{f_i}{\rho _i}{A_i}} }}{{\displaystyle\sum\nolimits_D {{\rho _i}{A_i}} }} , \end{split}$$ (4) 其中
$\xi$ 表示空气与界面内SF6的混合程度;${{f}}_{{i}}$ 和${{\rho}}_{{i}}$ 分别为区域D中第i个网格中SF6所占的质量分数和密度. 最开始边界扩散使得混合率快速上升, 其中气层厚度薄、初始扰动大的情况由于两气体接触面积大而扩散更快. 之后从ts1生成到冲击II前和从ts2生成到srs冲击II的过程中, 混合率均呈现增长趋势, 且前者的增速略高. 这一阶段初始扰动对混合率增长影响较小. 在sts冲击OI后混合率快速增大, 之后增速慢慢减小, 且初始扰动幅值越大增速越大, 越早开始减速. 而case 4中由于内外界面距离较大, 在tts冲击II后, 较长时间没有激波与界面作用, 因此混合率在sts冲击OI前增速几乎为零.由涡动力学方程可知二维条件下涡拉伸项为零, 但为表示小尺度能量分布, 定义湍动能为[27]
$$ TKE = \int_D {\frac{1}{2}\rho \left[ {{{\left( {u - \overline u } \right)}^2} + {{\left( {v - \overline v } \right)}^2}} \right]} {\rm d}s , $$ (5) 其中
$ \bar{u} $ 和$ \bar v $ 为坐标变换后流向与展向平均速度. 图14(a)为case 1湍动能分布图. 在is冲击OI后, 湍动能主要集中在ts1波峰及此处之后的激波上, 当ts1完全穿透II后, 湍动能主要分布在II与rrw1之间, 这是由于产生的rrw1与rs2,3使得这个区域激波结构复杂, 而最值出现在II上. 同时, 随着ts2向中心移动与相邻ts3的交叉, 湍动能逐渐增大, 且在激波聚焦时在中心处达到峰值. 在srs冲击II后, 湍动能主要分布在“尖钉”头部, 而“气泡”上很小. 之后“尖钉”头部两侧小涡的发展使头部湍动能减小而小涡处湍动能增大, 湍动能分布变得比较平均. 之后湍动能主要分布在tts冲击OI波谷后产生的“尖钉”头部到II“气泡”之间的区域内, 且最值出现在“尖钉”根部. 从图14(b)可以看到, case 2湍动能分布与case 1比较类似, 区别主要体现在tts穿透OI波谷后, “尖钉”到II上的“气泡”之间的湍动能较小, II上“尖钉”湍动能也较小. 从图14(c)可以看到, 与case 1相比较, case 3的内外界面“尖钉”头部的湍动能更大, II“尖钉”两端小涡发展更早更快. 由于srs先后冲击“尖钉”“气泡”, 因此t = 0.094 ms时存在一个靠近中心的转折点.由图14(d)可知, case 4中ts1完全穿透II后, 相邻TRS的相交使得激波结构更加复杂, II上湍动能比case 1高出近50%. II在srs穿透后开始向中心收缩, “尖钉”头部减速向中心移动, 一段时间后开始向外移动且在发展一段时间后分布变得平均. 由于内OI间距增加, sts与tts冲击OI时间间隔较小, 因此两部分湍动能较大的区域间隔时间更短. 由图14(e)发现, case 5中ts1更早地冲击II, 产生的湍动能较低. 界面之间耦合效应使得后期OI处湍动能较低, 仅在tts冲击OI时出现峰值.
3.3 DMD分析
通过DMD方法对不同时刻涡量场提取快照序列进行分析可以得到界面上涡量相干结构变化, 并且取前4个特征值进行分析. 图15(a)—(e)为cases 1—5情况下频谱分布情况, 其中实部ωr > 0时说明模态为增长模态, ωr < 0则为衰减模态, 而ωi的绝对值越大频率越大.
图16展示了cases 1—5情况下对涡量分布进行DMD分解所得模态中, 选定的4个动模态实数部分, 由于虚部相当于将实部相位平移90°因此不做考虑, 其中正值代表逆时针旋转涡, 负值则代表顺时针旋转涡. 结合频谱分布图可以得知, cases 1—3的第一动模态DM1表示的流动结构振荡频率低且随时间增长, 代表了流场中的主要相干结构. 在中心处有六个主干结构, 且每个结构均有反对称的两个分支, 代表的是随着时间逐渐增强的内界面“尖钉”结构, 反映了后期流动的整体结构, 而且初始扰动越大, 涡量越集中在“尖钉”的两侧涡结构上, 解释了“尖钉”和“气泡”发展程度不同的原因. 而DM2中, 主干的每个分支上有两个反向旋转的涡量, 反映“尖钉”两侧涡结构向外发展的过程. 同时外界面涡量方向与DM1相反, 反映的是sts冲击外界面后界面涡量反向增长的过程. 而DM3和DM4表示的流动结构振荡频率较高增长较慢, 分布在“尖钉”“气泡”上的快速交替的正负涡量反映相邻“尖钉”两侧涡结构相互混合及内界面受到srs冲击时界面涡量反向增长的过程, 且振荡频率越高, 涡量交替越密集. 在case 4中, DM1表示的是振荡频率低且随时间衰减的流动结构, 同样有六个主干结构, 且每个结构有两个分支, 代表的是随着时间逐渐衰减的内界面“尖钉”两侧涡结构. 同时也有涡量分布在外界面上, 反映了流场后期结构. DM2中每个分支上有三个部分, 为正负交替的涡量, 与DM1相比衰减程度更大, 反映了流场后期内界面上“尖钉”发展较慢. 而DM3和DM4表示的流动结构振荡频率较高, 反映的是内界面相邻“尖钉”两侧涡结构在“气泡”中相互混合及外界面涡量反向增长的过程. 而case 5结构与之前有所不同, DM1表示的流动结构振荡频率低且随时间增长, 涡量主要分布在外界面上, 反映了流场后期涡量集中在外界面“尖钉”并随时间逐渐增强. 而DM2—4则均表示内外界面振荡频率低且随时间衰减的流动结构, 主要反映了流场后期内外界面之间密集的正负涡量交替.
4. 结 论
本文基于有限体积法和可压缩多组分Navier-Stokes方程, 并结合5阶WENO格式以及结构化自适应网格加密技术, 数值研究了汇聚激波冲击双层重气柱界面过程, 并研究了初始扰动幅值及气层厚度对界面不稳定性的影响. 得到结论如下:
1) 初始扰动幅值越大, 透射激波在重气柱内更容易生成马赫反射结构, 内界面受到冲击时与激波夹角更大, 反射激波与界面作用产生的反向涡量使其向内凹陷更加明显, 导致“尖钉”更早发展, 二次反射激波先后冲击“尖钉”“气泡”产生两道反射激波, 从而在中心发生更多次激波聚焦. 同时二次透射激波更早地冲击外界面波谷, 使得外界面“尖钉”更早形成, 内界面“气泡”也因此更早向外快速发展. 因此界面环量、混合率的增速更大, 流场后期湍动能也越高.
2) 当气层厚度较大时, 透射激波在向内传播的过程中马赫反射结构发展时间更长, 内界面受其冲击后与case 1中相位相差π. 生成的内界面“尖钉”间距很小, 中心的高压与空间限制使“尖钉”头部无法向内发展, 两侧涡结构只能在“气泡”内部向外发展. 此时耦合效应小, 外界面对内界面发展影响较小. 而气层厚度较小时, 激波结构与之前没有明显区别, 但由于耦合效应较强, 内界面随外界面波谷一同向内移动, 直到外界面“尖钉”生成后, 内界面“尖钉”“气泡”才开始发展. 因此界面环量始终较少且增速小, 湍动能较低, 而由于两气体接触面积较大, 混合率增长较大.
3) 通过对涡量进行DMD分析, cases 1—4流场后期流动结构主要有三个部分, 首先是中心处有六个振荡频率低, 且随时间增长的主干结构代表的是随着时间逐渐增强的内界面“尖钉”, 第二个是分支上正负交替涡量, 代表的内界面相邻“尖钉”两侧涡结构相互混合, 第三个是界面上振荡频率较高增长较慢的流动结构, 反映了内外界面受到向外激波冲击后涡量反向增长的过程. 而case 5的结构有所不同, 一是外界面上有振荡频率低且随时间增长的流动结构, 反映了流场后期涡量集中在外界面“尖钉”结构并随时间逐渐增强, 二是内外界面振荡频率较高且随时间衰减的正负交替涡量, 反映的是内外界面的相互影响与涡量交换.
[1] Yang J, Kubota T, Zukoski E E 1993 AIAA J. 31 854Google Scholar
[2] Cao L, Fei W L, Grosshans H, Cao N 2017 Appl. Sci. 7 880Google Scholar
[3] Lindl J D, McCrory R L, Campbell E M 1992 Phys. Today 45 32Google Scholar
[4] Richtmyer R D 1960 Commun. Pure Appl. Math. 13 297Google Scholar
[5] Fraley G 1986 Phys. Fluids 29 376Google Scholar
[6] Haehn N, Ranjan D, Weber C, Oakley J G, Anderson M H, Bonazza R 2010 Phys. Scr. T142 014067Google Scholar
[7] Haehn N, Weber C, Oakley J, Anderson M, Ranjan D, Bonazza R 2012 Shock Waves 22 47Google Scholar
[8] Luo X S, Wang M H, Si T, Zhai Z G 2015 J. Fluid Mech. 773 366Google Scholar
[9] 沙莎, 陈志华, 薛大文, 张辉 2014 物理学报 63 085205Google Scholar
Sha S, Chen Z H, Xue D W, Zhang H 2014 Acta Phys. Sin. 63 085205Google Scholar
[10] Mikaelian K O 1990 Phys. Rev. A 42 3400Google Scholar
[11] Lombardini M, Pullin D I 2009 Phys. Fluids 21 114103Google Scholar
[12] Si T, Long T, Zhai Z G, Luo X S 2015 J. Fluid Mech. 784 225Google Scholar
[13] Ding J C, Si T, Yang J M, Lu X Y, Zhai Z G, Luo X S 2017 Phys. Rev. Lett. 119 014501Google Scholar
[14] Ding J C, Li J M, Sun R, Zhai Z G, Luo X S 2019 J. Fluid Mech. 878 277Google Scholar
[15] Mikaelian K O 1995 Phys. Fluids 7 888Google Scholar
[16] Sun R, Ding J C, Zhai Z G, Si T, Luo X S 2020 J. Fluid Mech. 902 A3Google Scholar
[17] Li J M, Ding J C, Si T, Luo X S 2020 J. Fluid Mech. 884 R2Google Scholar
[18] 徐建于, 黄生洪 2019 力学学报 51 998Google Scholar
Xu J Y, Huang S H 2019 Chin. J. Theor. Appl. Mech. 51 998Google Scholar
[19] 梁煜, 关奔, 翟志刚, 罗喜胜 2017 物理学报 66 064701Google Scholar
Liang Y, Guan B, Zhai Z G, Luo X S 2017 Acta Phys. Sin. 66 064701Google Scholar
[20] Zhou Z B, Ding J C, Zhai Z G, Cheng W, Luo X S 2020 Acta Mech. Sin. 36 356Google Scholar
[21] Tang J G, Zhang F, Luo X S, Zhai Z G 2020 Acta Mech. Sin. 37 434Google Scholar
[22] 何惠琴, 翟志刚, 司廷, 罗喜胜 2016 计算物理 33 66Google Scholar
He H Q, Zhai Z G, Si T, Luo X S 2016 Chin. J. Comput. Phys. 33 66Google Scholar
[23] Fu Y W, Yu C P, Li X L 2020 AIP Adv. 10 105302Google Scholar
[24] Lombardini M, Hill D J, Pullin D I, Meiron D I 2011 J. Fluid Mech. 670 439Google Scholar
[25] Hill D J, Pullin D I 2004 J. Comput. Phys. 194 435Google Scholar
[26] Pantano C, Deiterding R, Hill D J, Pullin D I 2007 J. Comput. Phys. 221 63Google Scholar
[27] Henry D, Movahed P, Johnsen E 2015 Shock Waves 25 329Google Scholar
-
图 4 case 1的界面与激波结构演变过程示意图(ts, 透射激波; rs, 反射激波; rrw, 反射稀疏波; m, 马赫杆; T, 三波点; SF, 激波聚焦; srs, 二次反射激波; spike, “尖钉”结构; bubble, “气泡”结构; sts, 二次透射激波; trs, 三次反射激波; tts, 三次透射激波; 下文符号含义相同)
Figure 4. Evolution of the interface and shock wave structures of case 1 (ts, transmitted shock; rs, reflected shock; rrw, reflected rarefaction wave; m, Mach stem; T, triple point; SF, shock focusing; srs, the second reflected shock; spike, “spike” structure; bubble; “bubble” structure; sts, the second transmitted shock; trs, the third reflected shock; tts, the third transmitted shock. The meaning of these abbreviations is similar hereinafter).
表 1 不同双层重气柱几何参数表
Table 1. Structural parameters of cylinder of different cases.
Case R0/mm r0/mm a0/mm n λ/mm a0/λ Case 1 20 10 1 6 20.94 0.048 Case 2 20 10 0.5 6 20.94 0.024 Case 3 20 10 2 6 20.94 0.096 Case 4 20 5 1 6 20.94 0.048 Case 5 20 15 1 6 20.94 0.048 表 2 气体参数表
Table 2. Parameters of gases.
Gas γ M/(g·mol–1) ρ/(kg·m–3) Air 1.399 28.967 1.23 SF6 1.103 128.491 5.45 -
[1] Yang J, Kubota T, Zukoski E E 1993 AIAA J. 31 854Google Scholar
[2] Cao L, Fei W L, Grosshans H, Cao N 2017 Appl. Sci. 7 880Google Scholar
[3] Lindl J D, McCrory R L, Campbell E M 1992 Phys. Today 45 32Google Scholar
[4] Richtmyer R D 1960 Commun. Pure Appl. Math. 13 297Google Scholar
[5] Fraley G 1986 Phys. Fluids 29 376Google Scholar
[6] Haehn N, Ranjan D, Weber C, Oakley J G, Anderson M H, Bonazza R 2010 Phys. Scr. T142 014067Google Scholar
[7] Haehn N, Weber C, Oakley J, Anderson M, Ranjan D, Bonazza R 2012 Shock Waves 22 47Google Scholar
[8] Luo X S, Wang M H, Si T, Zhai Z G 2015 J. Fluid Mech. 773 366Google Scholar
[9] 沙莎, 陈志华, 薛大文, 张辉 2014 物理学报 63 085205Google Scholar
Sha S, Chen Z H, Xue D W, Zhang H 2014 Acta Phys. Sin. 63 085205Google Scholar
[10] Mikaelian K O 1990 Phys. Rev. A 42 3400Google Scholar
[11] Lombardini M, Pullin D I 2009 Phys. Fluids 21 114103Google Scholar
[12] Si T, Long T, Zhai Z G, Luo X S 2015 J. Fluid Mech. 784 225Google Scholar
[13] Ding J C, Si T, Yang J M, Lu X Y, Zhai Z G, Luo X S 2017 Phys. Rev. Lett. 119 014501Google Scholar
[14] Ding J C, Li J M, Sun R, Zhai Z G, Luo X S 2019 J. Fluid Mech. 878 277Google Scholar
[15] Mikaelian K O 1995 Phys. Fluids 7 888Google Scholar
[16] Sun R, Ding J C, Zhai Z G, Si T, Luo X S 2020 J. Fluid Mech. 902 A3Google Scholar
[17] Li J M, Ding J C, Si T, Luo X S 2020 J. Fluid Mech. 884 R2Google Scholar
[18] 徐建于, 黄生洪 2019 力学学报 51 998Google Scholar
Xu J Y, Huang S H 2019 Chin. J. Theor. Appl. Mech. 51 998Google Scholar
[19] 梁煜, 关奔, 翟志刚, 罗喜胜 2017 物理学报 66 064701Google Scholar
Liang Y, Guan B, Zhai Z G, Luo X S 2017 Acta Phys. Sin. 66 064701Google Scholar
[20] Zhou Z B, Ding J C, Zhai Z G, Cheng W, Luo X S 2020 Acta Mech. Sin. 36 356Google Scholar
[21] Tang J G, Zhang F, Luo X S, Zhai Z G 2020 Acta Mech. Sin. 37 434Google Scholar
[22] 何惠琴, 翟志刚, 司廷, 罗喜胜 2016 计算物理 33 66Google Scholar
He H Q, Zhai Z G, Si T, Luo X S 2016 Chin. J. Comput. Phys. 33 66Google Scholar
[23] Fu Y W, Yu C P, Li X L 2020 AIP Adv. 10 105302Google Scholar
[24] Lombardini M, Hill D J, Pullin D I, Meiron D I 2011 J. Fluid Mech. 670 439Google Scholar
[25] Hill D J, Pullin D I 2004 J. Comput. Phys. 194 435Google Scholar
[26] Pantano C, Deiterding R, Hill D J, Pullin D I 2007 J. Comput. Phys. 221 63Google Scholar
[27] Henry D, Movahed P, Johnsen E 2015 Shock Waves 25 329Google Scholar
Catalog
Metrics
- Abstract views: 3384
- PDF Downloads: 36