激波作用不同椭圆氦气柱过程中流动混合研究

李冬冬 王革 张斌

1)(哈尔滨工程大学航天与建筑工程学院,哈尔滨 150001)

2)(上海交通大学航空航天学院,上海 201100)

1 引 言

当激波作用于不同密度的流体界面时会发生折射和反射等现象,流体界面也会获得一定的加速,从而导致流场中产生一系列复杂流动现象,这种现象被称为Richtmyer-Meshkov不稳定性(Richtmyer-Meshkov instability,RMI)[1,2].许多自然现象和工程问题中存在这种现象,例如武器内爆、惯性约束核聚变[3,4]、超音速燃烧[5]、超新星爆发[6]等.对RMI的研究不仅会推动涡动力学、多相流、湍流等流体力学基础难题的研究,而且在工程领域和自然界中有直接应用价值.

近几十年来,国内外学者对激波与不同密度的气体界面相互作用问题中的RMI进行了大量的数值、实验和理论研究.Haas和Sturtevant[7]采用实验手段对弱激波与R22重气柱、气泡以及He轻气柱、气泡的作用过程进行了研究,对界面变形和运动特征进行了讨论和分析.Jacobs[8]首次采用平面激光诱导荧光(planar laser-induced fluorescence,PLIF)技术研究了低马赫数激波与氦气柱相互作用的混合问题,采用激波压缩后气柱面积的减少量来表征混合.Giordano和Burtschell[9]研究了较低马赫数下激波与气柱及气泡的相互作用,数值模拟结果和一维气体动力学理论分析表明,不同密度比的气泡受相同激波压缩后,体积压缩率近似是定值,其值与气泡和周围气体的比热之比有关.Ranjan等[10]通过数值模拟方法验证了Giordano和Burtschell[9]有关压缩率的结论:高马赫数激波与轻气泡作用时,Giordano等的结论失效,且随着马赫数以及密度比的增大,体积压缩率变化存在更大的振荡.Tomkins等[11]采用PLIF技术对马赫数1.2的激波与SF6重气柱的相互作用问题进行了实验研究,并以物质质量分数为基础的瞬态扩散率和总混合率表征混合效果,对SF6与空气的混合机理进行了深入的研究和讨论.Shankar等[12]采用局部人工扩散方法对Tomkins等[11]的实验进行了数值模拟,以总混合率为参考标准,研究了SF6初始分布和示踪物质丙酮的含量对结果的影响.国内沙莎等[13,14]采用五阶加权基本无振荡格式与沉浸边界法并结合大涡模拟技术的数值方法,对激波作用R22重气柱和SF6球形气泡所致的射流混合进行了研究,分析了入射激波以及反射激波在气泡界面聚焦诱导射流的过程,详细研究了不同反射距离下反射激波与重气体的作用过程及流场结构,此外对入射激波与梯形SF6重气柱相互作用过程中的复杂波系演化也进行了详细的分析[15].Bai等[16]以及廖深飞等[17]采用粒子图像测速技术对激波诱导单椭圆、双椭圆SF6重气柱的不稳定性问题进行了实验研究,获得了激波冲击下的速度场、涡量场和环量等,定量地表征了激波界面相互作用的RMI现象,并探讨了椭圆结构对界面演化的影响.Zhai等[18]采用VAS2D程序对激波作用Kr/SF6气泡过程中射流的形成进行了深入研究,结果表明,压力扰动和斜压-涡量的沉积是激波与气柱相互作用(shock bubble interaction,SBI)问题中射流形成的两个主要因素.随后采用实验方法研究了二维V形空气/SF6界面在入射激波和反射激波作用下的RMI发展规律[19].Fan等[20]同样借助VAS2D程序对激波作用于长方形、椭圆、三角形等5种不同形状的SF6气柱进行了数值模拟,分析了该过程中气柱形态、环量等的变化.Wang等[21]采用高速摄影技术和VAS2D数值模拟方法对激波作用不同形态的氦气柱进行了研究,对界面形态特征、环量等进行了比较和分析.黄熙龙等[22]采用PLIF技术对激波作用于椭圆气柱过程中的RMI问题进行了研究,分析了该过程中界面气体聚集、转移、消散等现象.

激波与气柱以及液滴的相互作用可加速其与外部气体的混合,从而提高特定条件下的燃烧性能,相关研究对提高超燃发动机的性能具有极其重要的作用,因而对激波加速气柱与外部气体混合机理的研究具有重要意义.本文基于双通量模型[23],结合五阶WENO方法求解二维多组分Navier-Stokes(N-S)方程,对激波与不同结构的椭圆氦气柱作用过程进行数值模拟,研究了椭圆氦气柱与周围介质之间的混合情况及几何构型对流动和混合的影响,以期为进一步研究混合燃烧效率提供一定的理论基础.

2 数值方法与计算模型

2.1 数值方法

考虑可压缩多组分的二维N-S方程组,其形式为

其中U为守恒参数组成的向量,F和G为对流通量向量,Fv和Gv为黏性及扩散通量向量,具体表达式为

这里ρ,p和E分别代表混合物密度、压力以及单位质量的总能量;u和v是混合物速度矢量在x和y方向上的分量;Yi是组分i的质量分数;Ng为物质种类总数;Jx,i和Jy,i分别为组分i的扩散在x和y方向的分量;qx和qy分别为热扩散量在x和y方向的分量;σ为偏应力张量.多组分流动的理想气体状态方程为

其中Mwi为物质i的摩尔分子量,Ru为通用气体常数,T为混合物温度.比总能量的表达式为

其中CPi为物质i的定压比热,采用美国国家航空航天局(National Aeronautics and Space Administration,NASA)提供的拟合多项式进行计算;h0fi为物质i的标准焓;H为物质的总焓;T0为参考温度.混合物的定压比热CP采用各物质的质量分数Yi和定压比热CPi加权平均得到.音速c是气相混合物的一个重要属性,定义为

对于理想气体音速c可以简化为

其中γ为混合物的比热,定义为

物质i的质量扩散通量Ji为

其中是物质i的扩散速度,其计算公式为

其中Xi为物质的摩尔分数,Di,mix为物质i在混合物中的混合平均扩散系数.

忽略Dufour效应,热扩散矢量q计算表达式如下:

其中hi为物质i的焓值,同样采用由NASA提供的多项式计算;λ为混合物的导热系数.

流体间正应力和剪应力分别为

其中µ为混合物的剪切黏性,κ为体积黏性.

本文采用简化的混合方法来计算混合物的输运性质参数 (ξ=µ,λ,κ)[24]:

对于物质剪切黏性µ、体积黏性κ和热传导系数λ,n分别为6,4/3和4.

物质的混合平均扩散系数Di,mix采用下式计算:

其中Dij为物质i在物质j中的二元扩散系数,Xj为物质j的摩尔分数.

单物质的剪切黏性、热传导系数、二元扩散系数基于Chapman-Enskog理论[24,25]和Lennard-Jones参数[26]计算获得,纯物质的体积黏性κ依据Ern和Giovangigli[24]的简化方法获得.

考虑到控制方程(1)各部分的物理性质,采用时间分裂的方法,将控制方程分裂成对流项和黏性项依次求解,分裂方式如下:

对流项的求解在空间上采用五阶WENO格式并结合双通量模型消减多物质比热不同引起的界面参数振荡的影响,时间上采用三阶Runge-Kutta格式[27];黏性项的求解中空间项采用四阶中心差分,时间上采用二阶Runge-Kutta-Chebyshev方法[28]进行积分完成.

以一维欧拉方程为例对对流项的求解方法(五阶WENO格式结合双通量模型)进行必要的描述.一维欧拉方程形式如下:

其半离散格式为

其中Ui为网格节点上的守恒变量,Fi±1/2是网格左右界面上的对流通量.定义离散算子~(U):

方程(15)可以写成如下形式:

文中采用了Abgrall和 Karni[23]提出的双通量模型消除由于变比热带来的压力和速度振荡.在双通量模型中,网格界面上的对流通量计算两次,分别标记为如图1所示,其中采用五阶WENO格式获得.当所有网格界面的两次通量计算完成后,如图2所示,空间算子~(U)可以改写为

图1 双通量模型在高阶格式中的应用Fig.1.Application of double flux model in high order scheme.

图2 网格界面对流通量Fig.2.Convective flux at cell face.

2.2 计算模型

文中采用如图3所示的计算模型,计算条件参照Haas和Sturtevant[7]的实验条件给定,气柱内轻气体为氦气,周围介质为空气,诱导激波马赫数为1.22,位于气泡右侧,波前气体静止,温度和压力分别为293 K和1 atm(1 atm=1.01325×105Pa),波后气体参数由Rankine-Hugoniot条件获得.由于对称性,选取模型的上半部分进行计算,上边界采用固体反射边界,对称轴采用对称边界条件,左右边界上采用无反射边界条件,各参数梯度为0.

以图3所示的计算模型为基础,并在网格无关性检验和算例验证可靠性的前提下,针对面积相同、几何构型不同(圆形,两个激波沿椭圆长轴作用于气柱,两个激波沿椭圆短轴作用于气柱)的氦气柱与激波的相互作用过程,分析界面变形、波系演化、界面结构参数(界面高度和长度)、气柱体积压缩率、总混合率、环量等的变化,探究气泡内介质和周围环境气体的混合机理和不同几何构型下气体混合的优劣.以半径2.5 cm的圆形截面气柱面积19.63 cm2为标准来设计气柱截面形态,如图4所示.为了方便表示,采用统一的标准椭圆方程来约束界面形状:

图3 SBI计算模型示意图Fig.3.A schematic of computational domain for SBI.

图4 气柱几何构型 (a)圆形;(b)激波沿椭圆长轴作用于气柱;(c)激波沿椭圆短轴作用于气柱Fig.4.Helium cylinders geometry:(a)Circle;(b)shock hitting along major axis of ellipse;(c)shock hitting along minor axis of ellipse.

不同几何构型的椭圆气柱对应的参数如表1所列,其中Geometry 1和Geometry 2中激波沿椭圆长轴作用气柱,Geometry 3为圆形气柱,Geometry 4和Geometry 5中激波沿椭圆短轴作用气柱.

表1 气柱边界方程参数Table 1.Parameters of gas cylinder boundary equation.

为防止虚假涡量的产生,采用有限厚度扩散层[29]获得更加光滑的初始气泡界面,(20)式为氦气质量分数初始分布函数,

其中fin,fout分别表示气柱内外氦气质量分数;d为椭圆上一点到椭圆中心的距离;r为计算域内任意一点到椭圆中心的距离;r0为椭圆中心位置;Cr为常数,值越大界面越光滑,所有算例中Cr=10000.在笛卡尔坐标系中,假定椭圆长轴或者短轴位于计算域的对称轴上,几何中心位于坐标原点,计算域内有一点坐标为(x0,y0),则r和d的表达式为

3 结果与讨论

3.1 网格无关性检验与算例验证

为了考察计算方法的可靠性,选择与参考文献[7]中激波作用氦气柱的实验结果进行对比.首先,以Geometry 3圆形气柱为例,采用三种不同尺寸的网格进行网格无关性检验,网格数量分别为粗网格(600×83,Grid 1)、中等网格(1200×166,Grid 2)、细网格(2400×330,Grid 3),相应的气柱内网格数量分别为80,160和320.

图5是以氦气质量分数显示的不同网格尺度下的界面变形情况.结果表明在三种网格尺度下界面变形过程基本相似,较细的网格体现出更多的界面细节特征.图6为不同网格尺度下界面三个特征点(特征点位置如图6中右下角所示)在激波作用He气泡过程中的位置变化(界面以氦气质量分数0.01作为标准),三种不同网格下所获得的结果基本相同.在后续的计算分析中为了获得较为清晰的流动细节均采用细网格(2400×330).

图5 网格无关性检验(左,600×83;中,1200×166;右,2400×330) (a)62µs;(b)240µs;(c)427µs;(d)674µsFig.5. Grid re finement test(left,600×83;middle,1200×166;right,2400×330):(a)62µs;(b)240µs;(c)427µs;(d)674µs.

图6 不同网格尺寸下特征点位置的变化Fig.6.Movement of characteristic points for different grids.

图7给出了激波作用于氦气泡前期,特征点位置的变化及与其他计算结果[30,31]的对比,图中曲线数据符合较好,表明本文所采用的数值方法具有较高的精度,可以用于激波作用氦气柱的研究中.

为进一步检验方法的准确性,将在细网格下获得的数值纹影结果与文献[7]实验采集到的纹影结果进行定性的对比,如图8所示,数值纹影图清晰地表明了作用过程中界面变形和波系演化过程,进一步说明所采用数值方法的可靠性和精度.

图7 激波作用于氦气泡前期,特征点位置的变化及与其他计算结果[30,31]的比较Fig.7.Space-time diagram for the interaction of a shock wave with a helium bubble at initial stage;comparisons with the results of Ref.[30,31].

图8 数值纹影与实验结果[7]的比较 (a)62µs;(b)102µs;(c)467µs;(d)674µsFig.8.Comparison between numerical schlieren and experimental results[7]:(a)62 µs;(b)102 µs;(c)467 µs;(d)674µs.

3.2 不同几何构型下氦气混合过程分析

3.2.1 界面形态和特征

1)界面变形和波系演化

激波作用于不同几何构型的椭圆氦气柱中界面形态及波系演化分别如图9—图13所示.在所有的算例中,激波由右向左作用于气柱,将激波接触界面最右端的时刻定义为0µs时刻.

图9为激波沿椭圆长轴作用于半长轴a=4 cm、半短轴b=1.5625 cm的氦气柱过程中界面和波系演化情况.当诱导激波与气柱最右端接触并沿椭圆气柱边界向下游运动的过程中,在界面处产生了不规则反射现象,如图9(a)和图14(a)所示.由于气柱内氦气的声阻抗较空气低,进入气柱内的透射激波(transmitted shock wave)传播速度要快于在空气中传播的诱导激波(incident shock wave),透射激波在界面处发生折射,在空气中形成一道新的激波(free-precursor shock wave),与诱导激波相遇后形成三叉激波结构.当透射激波运动到下游界面后,在空气中形成一道向下游传播的二次透射激波,同时在气柱内形成一道向上游界面传播的反射激波,如图9(b)所示.其后波系经过复杂的反射、折射及相互作用,变得极其复杂,从图9(c)和图9(d)可以清晰地观察到这一现象.当激波接触界面后,界面顶端会形成一个小的空气射流(air jet)(图9(b)),随着流动的发展,空气射流不断增长并向下游渗透,该过程中伴随着由RMI所致的主涡及小的速度剪切引起的次级涡的出现和发展(图9(c)—(f)),初始的小空气射流最终导致流动分离成两个独立发展的涡团(图9(g)和图9(h)).

图9 几何构型1数值纹影图(激波沿椭圆长轴作用气柱;a=4.0 cm,b=1.5625 cm)Fig.9.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 1(shock hitting along major axis of ellipse;a=4.0 cm,b=1.5625 cm).

图10和图11分别为激波沿椭圆长轴作用于半长轴a=3.25 cm、半短轴b=1.9231 cm的氦气柱与半径为2.5 cm的圆形氦气柱过程中界面和波系演化情况.其中波系演化模式及界面变形情况与几何构型1基本相似,但不难发现由于长轴的变短和短轴长度的增加(离心率变小),激波入射角度、压力梯度和密度梯度之间的方向发生了改变,一方面导致了界面前端空气射流形成时间增加,射流初期的形态也在增大,由其主导的流动不稳定性的发展有所滞后.不规则反射形成的激波结构中三叉点的位置及自由前体激波更贴近界面(图14(b)和图14(c)).

图12和图13分别为激波沿椭圆短轴作用于半长轴b=3.25 cm、半短轴a=1.9231 cm的氦气柱与半长轴b=4.0 cm、半短轴a=1.5625 cm的氦气柱过程中界面和波系演化情况.由于椭圆界面曲率的变化,诱导激波角度较小,初始压力梯度和密度梯度的夹角也较小,因而此时界面变形和波系演化与上述的三种几何构型表现不同.在波系演化上,激波作用后产生的是规则反射现象(图12(a)和图13(a)),图14(d)和图14(e)给出了详细的波系情况,激波作用于重/轻界面后,在氦气柱内产生了一道向下游传播的透射激波,同时反射出一道在空气中向上游传播的稀疏波.在界面变形方面,激波作用后,界面的前端并未迅速产生空气射流,相反地在激波压缩作用下,界面上游较大区域界面曲率变小,产生了类似平面结构的状态(图12(b)和图12(c),图13(b)—图13(e)),其后随着涡在这一平面结构的末端产生,界面弯曲并向下游渗透,流动继续发展,最终形成与上述三种情况类似的两个独立涡团结构(图12(d)—(h)、图13(f)—(h)).

图10 几何构型2数值纹影图(激波沿椭圆长轴作用于气柱;a=3.25 cm,b=1.9231 cm)Fig.10.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 2(shock hitting along major axis of ellipse;a=3.25 cm,b=1.9231 cm).

图11 几何构型3数值纹影图(激波作用于半径为2.5 cm的圆形气柱)Fig.11.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 3(shock hitting circular helium cylinder;a=b=2.5 cm).

图12 几何构型4数值纹影图(激波沿椭圆短轴作用于气柱;a=1.9231 cm,b=3.25 cm)Fig.12.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 4(shock hitting along minor axis of ellipse;a=1.9231 cm,b=3.25 cm).

图13 几何构型5数值纹影图(激波沿椭圆短轴作用于气柱;a=1.5625 cm,b=4.0 cm)Fig.13.Sequence of numerical schlieren images of shock-accelerated elliptic helium cylinder for Geometry 5(shock hitting along minor axis of ellipse;a=1.5625 cm,b=4.0 cm).

图14 不同形状下激波作用气柱20µs后的典型波系结构 (a)几何构型1;(b)几何构型2;(c)几何构型3;(d)几何构型4;(e)几何构型5Fig.14.Typical wave systems derived from numerical simulation for shock-accelerated elliptic cylinder at 20µs:(a)Geometry 1;(b)Geometry 2;(c)Geometry 3;(d)Geometry 4;(e)Geometry 5.

2)界面特征点位置变化

为了定量化地研究界面形态的变化,依据数值模拟得到的参数分布,以氦气质量分数0.01作为物质界面,得到界面三个特征点的位置(界面上游位置upstream、界面下游位置downstream、射流位置jet)变化如图15所示,图中同样给出了特征点的位置示意图.

图15 界面特征点位置及界面长度和高度的变化 (a)几何构型1;(b)几何构型2;(c)几何构型3;(d)几何构型4;(e)几何构型5Fig.15.Movement of distorted upstream,jet and downstream interfaces,and the interface length and height:(a)Geometry 1;(b)Geometry 2;(c)Geometry 3;(d)Geometry 4;(e)Geometry 5.

对于所有五种几何构型界面,上游界面位置upstream随时间的变化可以分为两个阶段:流动的初期,上游界面在诱导激波的作用下获得一个较大的速度ui;其后,随着透射激波在气柱下游界面的反射激波作用于上游界面后,上游界面重新获得一个新的速度uf. 在两个阶段上,几何构型1、几何构型2、几何构型3、几何构型4和几何构型5的平均速度ui/uf分别为196.3/134.2 m/s,193.4/120.2 m/s,181.3/105.2 m/s,153.1/93.5 m/s,136.9/83.3 m/s.下游界面在透射激波作用前基本保持静止,随后同样以一个基本恒定的速度沿流动方向发展,其平均运动速度分别为150.5,148.4,138.8,128.6,121.8 m/s,后期涡环(射流处)的速度分别为184.4,159.1,143.2,143.3,144.0 m/s.界面长度(length)定义为界面上下游之间的距离,在流动发展初期,由于激波的压缩作用,界面长度首先呈现变小的趋势,流动过程中存在一个长度基本不变的平台期,主要是因为此时界面上下游速度相近,其后由于界面下游速度高于上游速度,界面长度又有所上升,平台期的长度随着椭圆参数a值的减小而减小.界面高度(height)在透射激波作用到气柱上下缘之前基本不变,而后在透射激波的作用下获得一个向外的速度,使得气柱高度有所增加,增长速率随着椭圆参数a的减小而增加,其后界面高度变化主要取决于该过程中空气射流或者涡旋的发展,特别是对于激波沿椭圆短轴作用的几何构型3和4,后期由于涡旋在内侧产生,发展过程中不断弯曲内侧的界面,在独立涡团形成之间基本不对界面高度产生影响(图15(d)和图15(e)).

3)气柱体积压缩率

气泡体积变化可以很好地反映流动过程中不同介质的混合情况,因而常被作为SBI问题中的重要分析参数.Giordano和Burtschell[9]定义了以体积压缩率γc表征气泡体积随时间的变化,

图16 不同椭圆界面下气柱体积压缩率变化Fig.16.Volume compressibility for different elliptic interfaces.

上述的数值结果与分析表明,界面初始形状对SBI过程中界面变形和波系演化有很大的影响.界面变形和波系演化特点主要由椭圆结构自身曲率的变化和激波作用方向共同决定.激波强度一定时,激波在界面处的反射类型依赖于激波诱导角度θ(incident angleθ)的大小.由于椭圆界面各处曲率不同,针对上述五种椭圆构型及激波作用方向,采用一种简单的近似方法估计激波诱导角度,如图17所示.在这种近似方法下得到的诱导角度θ分别为68.6◦,59.4◦,45.0◦,30.6◦和21.3◦.分析结果显示在激波诱导角度θ为68.6◦,59.4◦和45.0◦时界面发生的是不规则反射现象,如图14(a)—(c)所示,而30.6◦和21.3◦的诱导角度下发生简单的规则反射现象.尽管研究的界面结构形式不同,Zhai等[32]在对正方形、等腰三角形和菱形气柱的研究中,同样指出在诱导角度为60.0◦,60.3◦和90.0◦时,激波在重/轻界面处的反射类型为不规则反射.

同样地,若将该角度用于衡量压力梯度和密度梯度方向的不一致性,分析界面形态的变化.激波沿长轴作用于气柱(此时不一致性较大,平均夹角分别为68.6◦,59.4◦),气柱前端产生空气射流,向氦气柱内浸入,进而主导流动和混合的发展;反之,当激波沿短轴作用于气柱(此时不一致性较小,平均夹角分别为30.6◦和21.3◦),由于激波的压缩作用,在上游界面形成平行于诱导激波的近平面结构,随后涡旋在该平面结构的末端产生并主导流动和混合的发展过程.界面特征点速度及结构参数也因此有所差异.

图17 激波诱导角度θ的近似方法Fig.17.Approximate method of shock incident angle θ.

对界面变形、特征点位置和特征点速度的分析表明,激波强度和气柱面积一定的情况下,激波沿椭圆长轴作用于气柱时,离心率越大,流动发展得越迅速,界面的变形及氦气与周围气体的混合越快;激波沿椭圆短轴的作用规律与此相反,并且整体流动混合也慢于前者.所有的几何构型下,体积压缩率变化趋势基本相同,最终稳定值的差异并不大,这与Giordano和Burtschell[9]关于体积压缩率的结论也基本相符(体积压缩率只与激波强度和气体比热比有关),认为气体体积压缩率并不能很好地反映气柱几何结构对流动和混合的影响.

3.2.2 环 量

激波作用于不连续气柱的过程中,由于斜压机制的作用会在物质界面处产生大量的涡旋.在这些沉积于界面的涡量的驱动下,气柱形态产生大的变形.由压力梯度和密度梯度不一致而产生的斜压涡w可以表示为

进一步对氦气柱与激波作用中产生的环量进行计算,结果如图18所示.数值模拟中环量Γ采用如下方式计算:

其中D表示计算域,w(x,y,t)表示垂直于区域D的涡量.

图18 不同椭圆界面下上半部分环量的变化 (a)正环量;(b)负环量;(c)总环量Fig.18.Circulation versus time for different elliptic interfaces:(a)Positive circulation;(b)negative circulation;(c)total circulation.

环量的计算中只考虑模型的上半部分,负环量在其中起主导作用.激波沿椭圆长轴作用于椭圆气柱和圆形气柱的过程中,正、负环量的变化过程基本相似.在作用的初期,由于激波和界面的作用,负环量迅速增加;其后随着气柱前缘空气射流的产生和发展,正、负环量以基本相似的趋势在增加,直到空气射流基本发展完全,在小的速度剪切引起的次级涡的共同作用下,此后环量的变化相对复杂.激波沿椭圆短轴作用于气柱的过程中环量的变化基本与上述过程相似,同样的负环量在激波的作用下迅速增加,但是由于椭圆自身的几何结构,此种作用方式下,激波与界面接触时间较短,因而导致负环量增加时间较短,增加值也较小,其后在其他波系的作用下缓慢增加,并且随着涡的产生,增长速率有所增加,对于几何构型4,小的速度剪切造成的次级涡在后期同样导致了环量值的下降.总体来看,在激波与氦气柱相互作用的过程中,尽管椭圆气柱结构不同,但环量的变化趋势基本相同,并且随着椭圆结构参数a的减小(从几何1到几何5),该过程中正、负环量所能达到的最大值(大小)在不断减小,总环量的稳定值也在不断减小.

结合对界面及波系的分析,正、负环量的变化主要经历三个时期:1)作用初期,环量主要由诱导激波与界面的压缩作用产生,此时正环量值缓慢增加,负环量值迅速增加;2)作用中期,环量主要受RMI产生的主涡控制,此时正、负环量值均迅速增加;3)作用后期,环量主要由RMI产生的主涡和小的速度剪切引起的次级涡共同控制,主涡倾向于增加环量值,次涡倾向于减小环量值.由于正负环量在中后期变化的一致性,总环量在整个发展过程中快速稳定在一个基本恒定的值上,其中几何构型1达到稳定期经历的时间最长,且最终稳定值大.环量的变化同样可以体现SBI过程中的激波压缩效应、RMI引起的主涡的产生发展及其对流动和混合的作用,定量地反映流动和混合机制.

3.2.3 总混合率

物质扩散率(混合率)在时间和空间上的分布可以很好地反映物质之间的混合特征.Tomkins等[11]以物质质量分数为基础,定义物质瞬时标量扩散率或混合率(instantaneous scalar dissipation rate or mixing rate)χ(x,t)为

总混合率由瞬时混合率进行空间积分得到,定义如下:

其中DHe为氦气在混合气体中的扩散系数,为初始氦气最大质量分数,YHe为氦气质量分数.

图19为不同椭圆界面下氦气总混合率随时间的变化.根据总混合率的定义可知其大小主要取决于氦气浓度梯度和界面面积的大小(二维中表现为接触界面长度).在作用初期,诱导激波一方面导致界面面积迅速减小,另一方面导致部分激波掠过的区域氦气浓度梯度增加,从总混合率迅速下降的结果可以得知在这一阶段,界面面积的减小占据着主导地位.其后随着空气射流和涡的产生,使得氦气与空气的接触界面面积增加,进而导致混合率又迅速变大.后期总混合率存在一个极大值,分析认为此后虽然界面面积在主涡的作用下依然增长,但流动中小的次级涡的出现,倾向于将氦气的浓度在空间上呈现均匀化分布的趋势,因此此时的浓度梯度有所减小,可以认为作用的后期氦气与空气处在一个充分混合的状态中,混合率的变化也因此较为复杂.总而言之,从总混合率的角度出发,同样可以将混合分为三个时期:1)激波压缩期,此阶段波系的作用占据主导地位,一方面导致界面面积急剧减小,另一方面使得部分激波掠过的区域氦气浓度梯度增加,总混合率呈下降趋势;2)RMI阶段,在此时期,RMI引起的主涡在流动过程中不断弯曲和伸展进入氦气的空气,使得接触界面面积增加,总混合率也随之迅速增加;3)充分混合阶段,此阶段以总混合率达到最大值为开始的标志,流动中出现的次级涡一方面使得氦气在核心涡区分布更加均匀,另一方面和主涡共同作用使得界面面积进一步有所增加,这一阶段混合率变化比较复杂,该阶段气体处于充分混和状态.不同几何构型的对比表明,总混合率的时间和空间分布虽然在整体上基本相似,但具体的流动和混合细节存在一定的差异,随着椭圆几何参数a的减小,流动和混合发展的越来越慢,进入充分混合阶段所需的时间越来越长,尤是激波沿椭圆短轴作用气柱的情况.

图19 不同椭圆界面下氦气的总混合率Fig.19.Total mixing rate versus time for different elliptic interfaces.

4 结 论

本文采用数值模拟方法研究了激波作用于面积相同,结构不同的椭圆型气柱过程中的流动混合情况,通过对该过程中气柱界面变形、波系演化、体积压缩率、环量和总混合率的分析,得到以下结论.

1)当激波沿椭圆长轴作用于椭圆气柱时,由于椭圆前缘附近界面的曲率特征,激波入射角度较大,密度梯度与压力梯度方向的不一致性也较大,此时诱导激波作用后界面发生不规则发射现象,界面最前端形成小的空气射流,主导其后的流动和混合过程,流动的RMI出现的也较早,小的速度剪切引起的Kelvin-Helmholtz不稳定性影响也越大,有利于氦气更大程度地与环境介质混合.并且此时离心率越大,流动发展越快.

2)当激波沿椭圆短轴作用于椭圆气柱时,此时短轴越短,激波诱导角度较小,密度梯度与压力梯度方向的不一致性也较小,诱导激波作用后界面发生规则反射现象,前缘界面被压缩成近平面结构,随后涡在平面结构的边缘产生,主导流动和混合的发展,这种作用情况下RM不稳定性发展较为缓慢.此时离心率越大,流动发展越慢.

3)对环量和总混合率的定量对比分析表明,SBI过程主要分为三个阶段:激波压缩期、RMI引起的主涡发展期、速度剪切引起的次级涡发展期.并且椭圆结构参数a越大,环量值和总混合率也越大.结合对界面变形的分析可以推断激波沿椭圆长轴作用于气柱,流动过程中的RMI产生越早发展越快,越有利于氦气与周围环境气体的混合,并且离心率越大,混合发展得越快越充分.