高速气体与椭圆柱云相互作用的数值研究

王 雅，蒋灵杰，邓小龙,2

（1北京计算科学研究中心，北京 100193；2弗吉尼亚大学，美国弗吉尼亚州夏洛茨维尔 22904）

高速气体和固体颗粒群的相互作用是一类典型的可压缩多相流问题，广泛存在于天文、自然灾害、工业安全、医疗工业和国防等领域，如超新星爆炸[1]、火山爆发[2]、粉尘爆炸[3]、无针注射[4]和炸弹爆炸[5]等。在高速颗粒流中，颗粒体积分数αd是一个重要参数[6]。当αd≪ 1 时，颗粒之间彼此远离，颗粒间的碰撞效应可以忽略不计[7]；当αd≥ 0.5时，颗粒之间彼此靠近，颗粒间的碰撞是其运动的主要机制，流体对固体颗粒的作用可以忽略不计[8]。本研究的固体体积分数为15%，属于0.001 ＜αd＜ 0.5范围[8]，此时流体与颗粒的相互作用以及颗粒之间的相互作用变得尤为重要[9-10]。

在物理实验方面，Rogue等[11]在垂直激波管中进行了激波与水平颗粒床相互作用的实验。Wagner等[12]使用多相流激波管（Multiphase shock tube,MST）对激波与颗粒帘的相互作用进行实验研究，观察到反射激波和透射激波的产生和传播，以及颗粒帘宽度扩张和传播过程。Wagner等[13]利用X射线测量技术改进了激波管，以观察在整个实验过程中颗粒帘内部体积分数的分布变化。Theofanous等[14-15]使用ASOS激波管研究激波与颗粒帘的相互作用，发现了此类实验中存在时间标度（Time scaling）现象，进而总结出颗粒帘宽度在实验前期加速扩张、实验后期匀速扩张的规律。

在数值模拟方面，研究颗粒流的常用方法有点颗粒模型法[16-17]和直接数值模拟（Direct numerical simulation,DNS）方法，DNS方法又包括基于网格重构的拉格朗日-欧拉移动网格法（Arbitray Lagrangian-Eulerian,ALE）[18]和浸入边界法（Immersed boundary method,IBM）[19-21]等。Zhu等[22]使用界面解析的DNS方法研究了不同形状（椭球和圆球）的颗粒在湍流通道中的流动。Zastawny等[23]使用改进的镜像浸入边界法进行DNS，推导了流动中4种非球形颗粒的阻力和升力系数以及扭矩系数。邹立勇等[24]实验研究了在马赫数为1.18的平面激波冲击作用下，双椭圆界面R-M不稳定性演化的动力学过程。2018年Jiang等[25]基于分层流模型[26-27]在笛卡尔背景网格下，进行了高速气体与二维圆柱云相互作用前期的系统性数值模拟，并对高速气体与三维圆球云的数值模拟进行了初步探究。Deng等[28]进一步探讨了高速气体与三维圆球云相互作用后期的相关规律，发现了与Theofanous等[14-15]实验中获得的时间标度类似的现象。

本研究基于分层流模型[26-27]，在前人[25,28]的基础上拓展研究，对平面激波与椭圆柱云的相互作用进行DNS，重点关注椭圆柱横截面不同长短轴之比（λ）以及椭圆柱横截面长轴与来流方向成不同角度（θ）时对流场的影响程度。图1为椭圆柱横截面几何示意图，其中a为椭圆柱横截面长轴长，b为椭圆柱横截面短轴长，λ=a/b表示长短轴之比，θ表示x轴（流场来流方向）与长轴之间的夹角。

图1 椭圆柱横截面几何示意图Fig.1 Illustration of the geometry for the cross-section of the elliptical cylinder

1 数值方法

本研究采用的数值方法是基于Chang和Liou[27]提出的分层流模型，其控制方程如下

式中：下标“g”表示气相；ρ为密度；v为速度矢量；n为控制体边界的单位外法向量；p为压力；E和H分别为每个控制体的总能量和总焓；Sg为控制体中气体相的表面积；Vg=αgVi，其中Vg为控制体中的气相体积，αg为气体的体积分数，Vi为控制体的体积。应用理想气体状态方程来封闭式（1）。采用有限体积法（Finite volume method,FVM）离散控制方程，空间重构使用三阶TVD格式，由于不同控制体的体积分数αg不都相同，每个控制体界面可以重构为气-气、固-固和气-固3个部分，如图2所示。其中气-气界面之间的通量使用AUSM+-up近似黎曼解法器[27,29-30]计算，时间推进采用三阶龙格库塔（Runge-Kutta）方法，计算网格使用笛卡尔网格。本研究主要针对激波与椭圆柱云相互作用的前期阶段，可认为此时椭圆柱固定不动。由于流场的流速较高，因此椭圆柱所受的合外力仅通过圆柱表面对压力积分获得，而忽略流体黏性对椭圆柱受力的影响。气-固界面采用滑移边界条件。使用MPI实现并行计算，以便进行大规模数值模拟。

2 平面激波与椭圆柱云相互作用研究

2.1 方法验证

通过测试网格的收敛性检验本方法的正确性。采用4种不同分辨率的网格进行数值模拟。表1列出了不同分辨率下x和y方向上的网格数Nx和Ny，以及解析椭圆短轴所使用的网格数nb，计算区域设置见图3(a)。椭圆圆心位于原点 (0 ,0)处，入射激波位于x=-0.04处，计算域为x∈[-0.140，0.140]，y∈[-0.084，0.084]，z∈[-0.020，0.020]。为了节省计算资源，当x∈[-0.056，0.056]且y∈[-0.028，0.028]时，网格间距是均匀的，且在该区域Δx=Δy；除此之外，x和y方向采用不等间距的拉伸网格，z方向上的均匀网格个数设置为2，且该方向上两侧为周期性边界条件。图3(a)的上下边界条件为周期性边界条件，左边界为入口边界条件，右边界为出口边界条件，入口边界和出口边界都按照气体的初始条件设置为固定值。气体初始条件为：波前 (p,T,u)R=(8.234 9×104Pa,294.9 K,0.0 m/s)，波后 (p,T,u)L=(2.548 9×105Pa,425.2 K,309.1 m/s)。

表1 网格收敛性分析实验中使用的4种网格Table 1 Four meshes used in the convergence analysis experiment

图3 网格收敛性分析实验示意图Fig.3 Illustration of the convergence analysis experiment

在此初始条件下可以产生马赫数Ma=1.67的运动激波。图3(b)显示了网格收敛性的数值模拟结果，主要考察了x方向上椭圆柱所受的外力Fx随时间的演化规律。从图3(b)中可以看出：随着nb的增加，数值模拟结果显示出很好的收敛性，且当nb=32时，数值模拟结果与nb=64时的结果吻合很好。因此，以下均使用nb=32的网格进行DNS。

2.2 直接数值模拟

设入射激波马赫数Ma=1.67，气体初始条件与2.1节中的初始条件相同，网格分辨率nb取为32。初始流场设置见图4，入射激波位于x=-2.5，计算域为x∈[-3.0，4.0]，y∈[-0.5，0.5]，z∈[-0.1，0.0]。为了节省计算资源，当x∈[-0.65，0.65]时，采用等间距网格，除此之外，x方向采用拉伸网格，y方向采用均匀网格，z方向上的网格个数设置为1。椭圆柱云位于x∈[-0.5，0.5]区域，其宽度L=1。椭圆柱个数Np=440，椭圆柱的横截面积均相同，其排布参照Jiang等[25]的排布方案。表2列出了λ为2、3和4时网格的设置，a为椭圆柱横截面长轴长，b为椭圆柱横截面短轴长，Δx为均匀网格区域网格的宽度，Nx和Ny分别表示x和y方向上的网格总数，N为整个计算区域内的网格总数。

图4 x-y平面计算区域设置示意图（右图为初始椭圆柱云分布图，蓝色表示低压区域，红色表示高压区域）Fig.4 Illustration of the computational domain setting in the x-y plane (The right plot shows the initial distribution of the elliptical cylinder cloudThe red and blue regions represent the high-pressure and low-pressure regions,respectively.)

表2 平面激波与椭圆柱云相互作用数值模拟使用的网格设置Table 2 Mesh settings in numerical simulation of the interaction between plane shock and elliptical column cloud

在数值模拟过程中，可以观察到类似于高速气体与圆柱云相互作用的物理现象[25]。图5显示了当λ=2、θ=0°时，无量纲时间t为1.3、1.5、1.8、2.4和3.5时流场中的压强分布，其中压强采用初始时刻波前压强（8.234 9×104Pa）进行了无量纲化处理。从图5中可以看出：在激波与椭圆柱云相互作用的过程中产生了一道向流场上游运动的反射激波和一道向流场下游运动的透射激波；当无量纲时间为2.4和3.5时，在椭圆柱云内部和流场下游部分区域流场扰动较大，原因是此区域存在激波的反射以及激波之间的相互作用。

图5 当λ=2、θ=0°时，不同无量纲时间下流场的无量纲压强分布Fig.5 Distributions of the dimensionless pressure at different dimensionless time when λ=2 and θ=0°

为了定量分析流场，先给出本研究所涉及的部分变量的定义[8]，变量 φ的体积平均可以定义为

式中：V为连续相和离散相的总体积。连续相的相平均（或雷诺平均）定义为

式中：Vc为V中连续相的体积。质量平均定义为

式中：下标i表示速度的3个方向表示速度的质量平均表示的平方。连续相的体积平均后的内能、动能和湍动能的定义[8]分别为

为了便于描述，将计算域分为3部分，即上游区域、椭圆柱云区域和下游区域，如图6(e)所示，其中椭圆柱云的边界分为椭圆柱云上游边界（UFC）和椭圆柱云下游边界（DFC）。图6展示了在t=3.5，λ为2、3、4时，椭圆柱横截面长轴与来流方向呈不同角度时的流场速度（u）、流场内能和流场动能在计算域中沿着x方向的分布。从图6中可以看出：随着θ从0°增大到135°，入射激波与椭圆柱云正面冲击的有效面积先增大后减小，椭圆柱云对入射激波的反射效果也先增强后减弱；当θ达到90°时，入射激波与椭圆柱云正面冲击的有效面积最大，椭圆柱云对入射激波的反射效果最强。从图6中来流速度u的分布可以看出，当θ=90°时，反射激波位置离UFC最远，而透射激波位置离DFC最近。从内能分布图中也可以看出，与其他角度相比，当θ=90°时，反射激波与UFC之间的区域中的内能相对于初始时刻的内能增加更大，而下游区域中内能增加最小。从动能分布中可以看出，当θ=90°时，入射激波对下游的影响相比其他角度来说更小，说明流场中的能量更多以内能形式保留在从反射激波面到UFC之间的区域中，而较少输入流场的下游。由于椭圆柱云的分布具有沿中心轴上下近似对称的性质，如图4所示，因此当θ=45°和135°时，流场中主流速度(u)、内能和动能的分布也具有相似性。

图6 t=3.5时不同λ下θ分别为0°、45°、90°、135°时的流场速度、流场内能和流场动能分布（灰色矩形区域表示椭圆柱云，RS、TS、UFC、DFC分别表示反射激波、透射激波、椭圆柱云上游边界、椭圆柱云下游边界）Fig.6 Distributions of the fluid velocity,fluid internal energy and fluid kinetic energy with different λ when θ equals to 0°,45°,90°,135° at dimensionless time t=3.5 (The gray rectangular regions stand for the elliptical cylinder cloud.Hereafter,RS,TS,UFC and DFC mean reflected shock,transmitted shock,the upstream front of elliptical column cloud,and the downstream front of elliptical column cloud,respectively.)

图7显示了在无量纲时间t=3.5且θ不变，λ分别为2、3、4时的主流速度、流场内能和流场动能在x方向上的分布。当θ=0°时，随着λ从2增加到4，相同面积的椭圆逐渐变得细长，入射激波与椭圆柱云正面冲击的有效面积减小，椭圆柱云对入射激波的反射效果减弱。从图7中来流速度u的分布来看，当θ=0°、λ=2时，反射激波间断面离UFC最远，透射激波的间断面离DFC最近。从内能分布图中也可以看出，随着λ的增大，反射激波与UFC之间的区域中的内能相对于初始时刻内能上升的幅度减小，传递到下游的内能相对于初始时刻内能逐渐增加，流场下游区域的内能占流场总内能的比例逐渐增加。从动能分布来看，随着λ的增加，动能整体增加。而当θ=90°时，随着λ的增大，反射激波与UFC之间的区域中的内能相对于初始时刻内能上升的幅度增加，传递到下游的内能相对于初始时刻内能逐渐减少，流场上游区域中的内能占流场总内能的比例逐渐增加，且此时流场中的动能随着λ的增加而减小，但差别不大。当θ=45°时，不同的λ下，入射激波与椭圆柱云正面冲击的有效面积差异较小，因此λ的不同对主流速度、流场内能和流场动能的影响较小。

图7 t=3.5时不同θ下λ分别为2、3、4时，流场速度、流场内能和流场动能的分布（灰色矩形区域表示椭圆柱云）Fig.7 Distributions of the fluid velocity,fluid internal energy and fluid kinetic energy with different θ,when λ equals to 2,3,4,at dimensionless time t=3.5,where the gray rectangular regions stand for the elliptical cylinder cloud

图8显示了在无量纲时间t=3.5时，相同θ下，当λ从2增加到4，沿x和y方向上的RMS速度（u′′，v′′）和湍动能的分布，虚线之间的区域表示椭圆柱云。可以观察到，当θ为0°、45°和135°时，λ的变化对u′′、v′′和湍动能分布的影响较小，此时u′′、v′′和湍动能数值较大的区域分布在DFC附近。当θ=90°时，随着λ的增加，u′′、v′′和湍动能的分布区域逐渐收窄，此时湍动能数值较大的区域分布在UFC附近。

图8 t=3.5，θ=0°,45°,90°,135°时流场RMS速度u′′、v′′以及湍动能k在不同 λ下沿x方向的分布Fig.8 Distributions of the fluid RMS velocity u′′,v′′ and turbulent kinetic energy k in the x-direction at different λ,when θ is equal to 0°,45°,90°,135° at dimensionless time t=3.5

为了探究更详细的规律，需要对流场内能、流场动能和流场湍动能进行定量分析，将计算域x∈[-3.0,4.0]分为3部分，分别是计算域上游区域（x∈[-3.0,-0.5]）、椭圆柱云区域（x∈[-0.5,0.5]）和计算域下游区域（x∈[0.5,4.0]），在每个计算区域内对内能、动能和湍动能进行积分，如图9所示，图中标出了每个区域具体积分值以便分析。

从内能分布来看，对于不同的λ和θ，在上游区域中的内能占整个流场中内能的比例最大。随着λ的增大，主要差异体现在：当θ=0°时，上游区域内能从23.8减小到20.6；而当θ=90°时，变化趋势则相反，内能从29.4增加到32.9，同时内能在此区域中的值比其他角度在此区域中的值都大，在椭圆柱云区域，内能从4.2减小到1.9，在下游区域，内能从7.8减小到6.8；当θ=45°，135°时，在不同λ下，每个区域的内能对λ的变化不敏感。

从动能分布来看，当θ=0°时，随着λ的增大，动能在上游区域从1.18增加到1.54，在椭圆柱云区域从0.10增加到0.23，在下游区域从0.83增加到1.53；当θ=90°时，趋势则相反，随着λ的增大，动能在上游区域、椭圆柱云区域和下游区域都逐渐减小；当θ=45°，135°时，随着λ的增大，动能在3个区域的值缓慢增加。

从湍动能分布来看，当θ处于不同角度下，湍动能在上游的值均很小，主要分布在椭圆柱云区域和下游区域。当θ=0°,45°,135°时，随着λ的增加，湍动能在椭圆柱云区域和下游区域的值都逐渐减小；当θ=90°时，湍动能主要分布在椭圆柱云区域。

图9 t=3.5时不同θ和λ下流场内能、流场动能和流场湍动能在计算域上游区域（x∈[-3.0,-0.5]）、椭圆柱云区域（x∈[-0.5,0.5]）和计算域下游区域（x ∈[0.5,4.0]）分布Fig.9 Distributions of the fluid internal energy,fluid kinetic energy and fluid turbulent kinetic energy at different θ and λ in three different regions,that is the upstream area of the domain x ∈[-3.0,-0.5],elliptical column cloud area x∈[-0.5,0.5],the downstream area of the domain x ∈[0.5,4.0]at dimensionless time t=3.5

图10 不同角度θ和不同λ下流场内能、流场动能和流场湍动能随无量纲时间t的变化Fig.10 Variations of the fluid internal energy,fluid kinetic energy and fluid turbulent kinetic energy with dimensionless time t at different θ and λ

图10展示了流场总内能、流场总动能和流场湍动能在t∈[1.2,3.5]的演化过程。因为t∈[0,1.2]时，流场中的激波未与椭圆柱云相互作用，因此该过程中的能量变化规律并不重要，可以忽略。这里主要讨论当x∈[1.2,3.5]时3个量的变化。3个量中流场总内能总是随着时间的增加而逐渐增大。当θ=90°时，流场总内能增长得最快；λ越大，θ越小，流场中总内能增长得越慢；到θ=0°时，流场总内能增长得最慢。总动能方面，当θ=0°、λ=3或者θ=0°、λ=4时，计算域中的总动能随时间的增加而逐渐增大，而其他情况下则随时间的增加而逐渐减小。当θ=90°时，计算域中的总动能减小得较快。就湍动能而言，总体上都随着时间的增加而逐渐增大，而且在相同θ下，λ越小，湍动能越大。综合来看，在相同λ的条件下，当θ=90°时，入射激波与椭圆柱云正面冲击的有效截面积最大，整个流场中总内能增长得最快，湍动能增长得最慢，而总动能减少得最快。

2.3 一维体积平均模型

本节讨论改进Jiang等[25]的一维体积平均模型与DNS的结果对比。改进后的一维体积平均模型的具体形式如下

式中：Cd为人工有效阻力系数，αd为固体体积分数。表3给出了在t=3.5时不同λ和θ下最优Cd的取值，其中θ=90°时以整体拟合趋势为准，其他角度则以一维体积平均模型的反射激波和透射激波位置与DNS结果的激波位置拟合最优为准。图11展示了λ=2，3，4，且θ=0°，45°，90°和135°时的速度分布。从图11中可以看出，一维体积平均模型与当前DNS结果的拟合效果较好。

图11 t=3.5时不同的λ和θ下一维体积平均模型与DNS拟合结果Fig.11 Fitting results of the 1-D volume-averaged model and DNS at different λ and θ when the dimensionless time is equal to 3.5

表3 人工有效阻力系数Cd的最优取值Table 3 Optimal values of artificial effective drag coefficient Cd

利用表3中的数据得到如图12所示的Cd分布。由图12可以看出：θ∈[0°,45°]时，λ越大，Cd越小；当θ∈[60°,120°]时，λ越大，Cd越大；当θ=90°时，Cd比相同λ、其他角度下的Cd都要大。在目前拟合的θ和λ范围内，Cd最小值出现在图12中A处，最大值出现在图12中B处。

图12 人工有效阻力系数Cd的最优取值分布Fig.12 Distribution for the optimal value of artificial effective drag coefficient Cd

3 结 论

研究了不同椭圆柱横截面长短轴之比λ和不同椭圆柱横截面长轴与来流方向所成角度θ对流场的影响，得到如下结论。

（1）在t=3.5时，保持λ不变，θ从0°增大到90°，激波与椭圆柱云正面冲击的有效横截面积逐渐增大，反射激波的位置逐渐远离椭圆柱云，透射激波的位置逐渐靠近椭圆柱云，此时流场中的能量更多以内能形式保留在从反射激波面到椭圆柱云上边界之间的区域中。

（2）在t=3.5时，保持θ不变，当θ=0°时，随着λ从2增大到4，反射激波的位置逐渐靠近椭圆柱云，透射激波的位置逐渐远离椭圆柱云，流场下游区域的内能占流场总内能的比例逐渐增加。θ=90°时，λ从2增加到4，此时流场上游区域中的内能占流场总内能的比例逐渐增加。θ=45°的内能分布和动能分布对λ的变化不敏感。

（3）当t=3.5、θ=90°时，随着λ的增加，u′′、v′′和湍动能的分布区域逐渐收窄，此时湍动能数值较大的区域分布在椭圆柱云上边界附近。而其他角度下，u′′、v′′和湍动能的分布对λ的变化不敏感，且u′′、v′′和湍动能数值较大区域分布在椭圆柱云下边界附近。

（4）当t∈[1.2,3.5]时，流场总内能随时间逐渐增大，相同λ下，当θ=90°时，流场总内能增长最快，θ=0°时，流场总内能增长最慢。当θ=0°、λ=3或者θ=0°、λ=4时，计算域中的总动能随时间的增加而逐渐增大，而其他情况则随时间的增加而逐渐减小。当θ=90°时，计算域中的总动能减小较快。湍动能随时间的增加而逐渐增大，相同θ下，λ越小，湍动能越大。在相同λ的条件下，当θ=90°时，动能转换成内能的效率最高。

（5）一维体积平均模型在采用合适的有效人工阻力系数Cd时，可以较好地拟合当前DNS结果。