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

王雅; 蒋灵杰; 邓小龙

doi:10.11858/gywlxb.20190748

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

doi: 10.11858/gywlxb.20190748

王雅^1,,
蒋灵杰¹,
邓小龙^{1, 2,*, ,}

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

基金项目: 中国工程物理研究院院长基金（201501043）；国家自然科学基金（U1530401）

详细信息

作者简介:
王　雅（1994－），女，学士，主要从事流体力学研究. E-mail: wangya@csrc.ac.cn

通讯作者:
邓小龙（1981－），男，博士，特聘研究员，主要从事流体力学研究.E-mail: xiaolong.deng@csrc.ac.cn

中图分类号: O359; O354.5
计量
- 文章访问数: 7677
- HTML全文浏览量: 3330
- PDF下载量: 28
出版历程
- 收稿日期: 2019-03-25
- 修回日期: 2019-05-07

Numerical Study of the Interaction between High-Speed Gas and Elliptical Column Cloud

WANG Ya^1
,,
JIANG Lingjie¹,
DENG Xiaolong^{1, 2,*
, ,}

1.
Beijing Computational Science Research Center, Beijing 100193, China
2.
University of Virginia, Charlottesville 22904, VA, USA

摘要

摘要: 高速颗粒流在天文、自然灾害、工业安全、医疗工业和国防等领域有着重要应用。采用基于分层流模型的直接数值模拟方法，对平面激波与椭圆柱云的相互作用进行数值研究，重点关注椭圆柱横截面的不同长短轴之比和椭圆柱横截面长轴与来流方向所成角度对流场的影响，从气体来流方向上的速度、x轴和y轴方向上的均方根速度、动能、内能和湍动能的分布上进行分析，对能量在计算域的上游区域、椭圆柱云区域和下游区域进行定量分析。同时针对椭圆柱改进了一维体积平均模型，利用该模型拟合了由直接数值模拟得到的反射激波和透射激波位置，获得了最适配的一维体积平均模型中的人工有效阻力系数，并探讨此系数的分布规律。
- 激波 /
- 椭圆柱 /
- 直接数值模拟 /
- 一维体积平均模型
Abstract: High-speed particle-laden flow has important applications in astronomy, natural disasters, industrial safety, medical industry, and national defense. In this work, a direct numerical simulation method based on the stratified flow model is used to study the interaction between a planar shock wave and an elliptical column cloud. The influence of the aspect ratio and the tilt angle, the distributions of the flow velocity, RMS velocities along x axis, kinetic energy, internal energy, and turbulent kinetic energy are analyzed; energy values in the upstream region, the elliptical column cloud region and the downstream region of the computational domain are quantitatively analyzed. The 1-D volume-average model is improved for elliptical columns. Based on this model, the appropriate artificial effective drag coefficients are decided by fitting the positions of the reflected shock and the transmitted shock from the direct numerical simulation results, and the distribution of the optimal artificial effective drag coefficients is also discussed.
- shock wave /
- elliptical column /
- direct numerical simulation /
- 1-D volume-averaged model

HTML全文

高速气体和固体颗粒群的相互作用是一类典型的可压缩多相流问题，广泛存在于天文、自然灾害、工业安全、医疗工业和国防等领域，如超新星爆炸^[1]、火山爆发^[2]、粉尘爆炸^[3]、无针注射^[4]和炸弹爆炸^[5]等。在高速颗粒流中，颗粒体积分数α_d是一个重要参数^[6]。当α_d $\ll 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 椭圆柱横截面几何示意图

Figure 1. Illustration of the geometry for the cross-section of the elliptical cylinder

下载: 全尺寸图片幻灯片

1. 数值方法

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

$\left\{ {\begin{aligned} & {\frac{\partial }{{\partial t}}\int_{{V_{\rm g}}} {{\rho _{\rm g}}} {\rm d}{V_{\rm g}} + \oint_{{S_{\rm g}}} {\left( {{\rho _{\rm g}}{{ v}_{\rm g}}} \right) \cdot { n}{\rm d}{S_{\rm g}} = 0} }\\ & {\frac{\partial }{{\partial t}}\int_{{V_{\rm g}}} {{\rho _{\rm g}}}{{{ v}_{\rm g}}} {\rm d}{V_{\rm g}} + \oint_{{S_{\rm g}}} {\left( {{\rho _{\rm g}}{{ v}_{\rm g}}{{ v}_{\rm g}}} \right)} \cdot { n}{\rm d}{S_{\rm g}} + \oint_{{S_{\rm g}}} p { n}{\rm d}{S_{\rm g}} = 0}\\ & {\frac{\partial }{{\partial t}}\int_{{V_{\rm g}}} {{\rho _{\rm g}}} {E_{\rm g}}{\rm d}{V_{\rm g}} + \int_V p \frac{{\partial {\alpha _{\rm g}}}}{{\partial t}}{\rm d}V + \oint_{{S_{\rm g}}} {\left( {{\rho _{\rm g}}{H_{\rm g}}{{ v}_{\rm g}}} \right)} \cdot { n}{\rm d}{S_{\rm g}} = 0} \end{aligned}} \right.$

(1)

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

图 2 分层流模型示意图^[25]（i–1、i和i + 1表示网格索引号，界面

$\overline {ab}$ 和

$\overline {ef}$ 为气-气界面，

$\overline {bc}$ 和

$\overline {fg}$ 为气-固界面，

$\overline {cd}$ 和

$\overline {gh}$ 为固-固界面）

Figure 2. Illustration of the stratified flow model^[25] (Where i –1, i and i + 1 are the indexes of the cell.

$\overline {ab}$ and

$\overline {ef}$ are the interfaces between the gas phases,

$\overline {bc}$ and

$\overline {fg}$ are the interfaces between the gas phase and solid phase,

$\overline {cd}$ and

$\overline {gh}$ are the interfaces between the solid phases.)

下载: 全尺寸图片幻灯片

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

2.1 方法验证

通过测试网格的收敛性检验本方法的正确性。采用4种不同分辨率的网格进行数值模拟。表1列出了不同分辨率下x和y方向上的网格数N_x和N_y，以及解析椭圆短轴所使用的网格数n_b，计算区域设置见图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 × 10⁴ Pa, 294.9 K, 0.0 m/s)，波后(p, T, u)_L = (2.548 9 × 10⁵ Pa, 425.2 K, 309.1 m/s)。

表 1 网格收敛性分析实验中使用的4种网格

Table 1. Four meshes used in the convergence analysis experiment

Mesh	n_b	N_x	N_y
1	8	112	64
2	16	224	128
3	32	448	256
4	64	896	512

下载: 导出CSV

| 显示表格

图 3 网格收敛性分析实验示意图

Figure 3. Illustration of the convergence analysis experiment

下载: 全尺寸图片幻灯片

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

2.2 直接数值模拟

设入射激波马赫数Ma = 1.67，气体初始条件与2.1节中的初始条件相同，网格分辨率n_b取为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。椭圆柱个数N_p = 440，椭圆柱的横截面积均相同，其排布参照Jiang等^[25]的排布方案。表2列出了λ为2、3和4时网格的设置，a为椭圆柱横截面长轴长，b为椭圆柱横截面短轴长，Δx为均匀网格区域网格的宽度，N_x和N_y分别表示x和y方向上的网格总数，N为整个计算区域内的网格总数。

表 2 平面激波与椭圆柱云相互作用数值模拟使用的网格设置

Table 2. Mesh settings in numerical simulation of the interaction between plane shock and elliptical column cloud

λ	N_p	a	b	Δx/10^–4	N_x	N_y	N/10⁶
2	440	0.029 44	0.014 72	4.59	3 894	2 176	8.5
3	440	0.036 04	0.012 02	3.75	4 763	2 666	12.7
4	440	0.041 64	0.010 41	3.25	5 500	3 078	16.9

下载: 导出CSV

| 显示表格

图 4 x-y平面计算区域设置示意图（右图为初始椭圆柱云分布图，蓝色表示低压区域，红色表示高压区域）

Figure 4. Illustration of the computational domain setting in the x-y plane (The right plot shows the initial distribution of the elliptical cylinder cloud. The red and blue regions represent the high-pressure and low-pressure regions, respectively.)

下载: 全尺寸图片幻灯片

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

图 5 当λ = 2、θ = 0°时，不同无量纲时间下流场的无量纲压强分布

Figure 5. Distributions of the dimensionless pressure at different dimensionless time when λ = 2 and θ = 0°

下载: 全尺寸图片幻灯片

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

$\bar \phi \equiv \frac{1}{V}\mathop \int\nolimits_V \phi {\rm{d}}V$

(2)

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

$\left\langle \phi \right\rangle \equiv \frac{1}{{{V_{\rm{c}}}}}\int\nolimits_{{V_{\rm{c}}}} \phi {\rm{d}}V$

(3)

式中：V_c为V中连续相的体积。质量平均 ${\tilde \phi }$ 定义为

$\tilde \phi \equiv \frac{{\overline {\rho \phi } }}{{\bar \rho }} = \frac{{\left\langle {\rho \phi } \right\rangle }}{{\left\langle \rho \right\rangle }}$

(4)

均方根（RMS）速度 ${u''_i}$ 定义为

${u''_i} \equiv \sqrt {\widetilde {u_i^2} - \tilde u_i^2}$

(5)

式中：下标i表示速度的3个方向， $\widetilde {u_{{i}}^2}$ 表示速度 ${u_i^2}$ 的质量平均， $\tilde u_i^2$ 表示 $\tilde u_i$ 的平方。连续相的体积平均后的内能、动能和湍动能的定义^[8]分别为

${E_{\rm i}} = \left\langle {{\rho _{\rm c}}} \right\rangle \tilde e = \frac{{\left\langle p \right\rangle }}{{\gamma - 1}}$

(6)

${E_{\rm{k}}} = \frac{{\left\langle {{\rho _c}} \right\rangle {{\tilde u}_i}{{\tilde u}_i}}}{2}$

(7)

$k = \frac{{{{\left\langle {{\rho _{\rm{c}}}{{u''_i}}u''_i} \right\rangle }}}}{2}$

(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分别表示反射激波、透射激波、椭圆柱云上游边界、椭圆柱云下游边界）

Figure 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时，流场速度、流场内能和流场动能的分布（灰色矩形区域表示椭圆柱云）

Figure 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方向的分布

Figure 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 \in \left[ {-3.0, 4.0} \right]$ 分为3部分，分别是计算域上游区域（ $x \in \left[ {-3.0, -0.5} \right]$ ）、椭圆柱云区域（ $x \in \left[ {-0.5, 0.5} \right]$ ）和计算域下游区域（ $x \in \left[ {0.5, 4.0} \right]$ ），在每个计算区域内对内能、动能和湍动能进行积分，如图9所示，图中标出了每个区域具体积分值以便分析。

图 9 t = 3.5时不同θ和λ下流场内能、流场动能和流场湍动能在计算域上游区域（

$x \in \left[ { - 3.0, - 0.5} \right]$ ）、椭圆柱云区域（

$x \in \left[ { - 0.5,0.5} \right]$ ）和计算域下游区域（

$x \in \left[ {0.5, 4.0} \right]$ ）分布

Figure 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 \in \left[ { - 3.0, - 0.5} \right]$ , elliptical column cloud area

$x \in \left[ { - 0.5,0.5} \right]$ , the downstream area of the domain

$x \in \left[ {0.5, 4.0} \right]$ at dimensionless time t = 3.5

下载: 全尺寸图片幻灯片

从内能分布来看，对于不同的λ和θ，在上游区域中的内能占整个流场中内能的比例最大。随着λ的增大，主要差异体现在：当θ = 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°时，湍动能主要分布在椭圆柱云区域。

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

图 10 不同角度θ和不同λ下流场内能、流场动能和流场湍动能随无量纲时间t的变化

Figure 10. Variations of the fluid internal energy, fluid kinetic energy and fluid turbulent kinetic energy with dimensionless time t at different θ and λ

下载: 全尺寸图片幻灯片

2.3 一维体积平均模型

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

$\frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{c}}}\left\langle \rho \right\rangle } \right) + \frac{\partial }{{\partial x}}\left( {{\alpha _{\rm{c}}}\left\langle \rho \right\rangle \tilde u} \right) = 0$

(9)

$\frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{c}}}\left\langle \rho \right\rangle \tilde u} \right) + \frac{\partial }{{\partial x}}\left[ {{\alpha _{\rm{c}}}\left( {\left\langle \rho \right\rangle {{\tilde u}^2} + \left\langle p \right\rangle } \right)} \right] = \left\langle p \right\rangle \frac{{\partial {\alpha _{\rm{c}}}}}{{\partial x}} - \frac{{2\sqrt {{\lambda ^2}{{\sin }^2}\theta + {{\cos }^2}\theta } }}{{{\text{π}}a}}{C_{\rm{d}}}{\alpha _{\rm{d}}}\left\langle \rho \right\rangle \left| {\tilde u} \right|\tilde u$

(10)

$\frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{c}}}\left\langle \rho \right\rangle {{\tilde e}_T}} \right) + \frac{\partial }{{\partial x}}\left( {{\alpha _{\rm{c}}}\left\langle \rho \right\rangle \tilde u{{\tilde h}_T}} \right) = 0$

(11)

$\left\langle \rho \right\rangle {{\tilde e}_T} = \left\langle \rho \right\rangle \tilde e + \frac{1}{2}\left\langle \rho \right\rangle {{\tilde u}^2}$

(12)

$\left\langle \rho \right\rangle {{\tilde h}_T} = \left\langle \rho \right\rangle \tilde e + \frac{1}{2}\left\langle \rho \right\rangle {{\tilde u}^2} + \left\langle p \right\rangle$

(13)

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

表 3 人工有效阻力系数C_d的最优取值

Table 3. Optimal values of artificial effective drag coefficient C_d

λ	C_d
λ	θ = 0°	θ = 15°	θ = 30°	θ = 45°	θ = 60°	θ = 90°	θ = 135°
2	1.19	1.06	1.06	1.42	2.22	6.00	1.68
3	0.86	0.67	0.69	1.48	2.80	14.00	1.29
4	0.57	0.58	0.47	1.30	3.00	36.00	1.15

下载: 导出CSV

| 显示表格

图 11 t = 3.5时不同的λ和θ下一维体积平均模型与DNS拟合结果

Figure 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中的数据得到如图12所示的C_d分布。由图12可以看出： $\theta \in \left[ {{0^\circ },{{45}^\circ }} \right]$ 时，λ越大，C_d越小；当 $\theta \in \left[ {{60^\circ },{{120}^\circ }} \right]$ 时，λ越大，C_d越大；当θ = 90°时，C_d比相同λ、其他角度下的C_d都要大。在目前拟合的 $\theta$ 和λ范围内，C_d最小值出现在图12中A处，最大值出现在图12中B处。

图 12 人工有效阻力系数C_d的最优取值分布

Figure 12. Distribution for the optimal value of artificial effective drag coefficient C_d

下载: 全尺寸图片幻灯片

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 \in \left[ {1.2,{\rm{}}3.5} \right]$ 时，流场总内能随时间逐渐增大，相同λ下，当θ = 90°时，流场总内能增长最快，θ = 0°时，流场总内能增长最慢。当θ = 0°、λ = 3或者θ = 0°、λ = 4时，计算域中的总动能随时间的增加而逐渐增大，而其他情况则随时间的增加而逐渐减小。当θ = 90°时，计算域中的总动能减小较快。湍动能随时间的增加而逐渐增大，相同θ下，λ越小，湍动能越大。在相同λ的条件下，当θ = 90°时，动能转换成内能的效率最高。

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

图 1 椭圆柱横截面几何示意图

Figure 1. Illustration of the geometry for the cross-section of the elliptical cylinder

下载: 全尺寸图片幻灯片

图 2 分层流模型示意图^[25]（i–1、i和i + 1表示网格索引号，界面 $\overline {ab}$ 和 $\overline {ef}$ 为气-气界面， $\overline {bc}$ 和 $\overline {fg}$ 为气-固界面， $\overline {cd}$ 和 $\overline {gh}$ 为固-固界面）

Figure 2. Illustration of the stratified flow model^[25] (Where i –1, i and i + 1 are the indexes of the cell. $\overline {ab}$ and $\overline {ef}$ are the interfaces between the gas phases, $\overline {bc}$ and $\overline {fg}$ are the interfaces between the gas phase and solid phase, $\overline {cd}$ and $\overline {gh}$ are the interfaces between the solid phases.)

下载: 全尺寸图片幻灯片

图 3 网格收敛性分析实验示意图

Figure 3. Illustration of the convergence analysis experiment

下载: 全尺寸图片幻灯片

图 4 x-y平面计算区域设置示意图（右图为初始椭圆柱云分布图，蓝色表示低压区域，红色表示高压区域）

Figure 4. Illustration of the computational domain setting in the x-y plane (The right plot shows the initial distribution of the elliptical cylinder cloud. The red and blue regions represent the high-pressure and low-pressure regions, respectively.)

下载: 全尺寸图片幻灯片

图 5 当λ = 2、θ = 0°时，不同无量纲时间下流场的无量纲压强分布

Figure 5. Distributions of the dimensionless pressure at different dimensionless time when λ = 2 and θ = 0°

下载: 全尺寸图片幻灯片

图 6 t = 3.5时不同λ下θ分别为0°、45°、90°、135°时的流场速度、流场内能和流场动能分布（灰色矩形区域表示椭圆柱云，RS、TS、UFC、DFC分别表示反射激波、透射激波、椭圆柱云上游边界、椭圆柱云下游边界）

Figure 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时，流场速度、流场内能和流场动能的分布（灰色矩形区域表示椭圆柱云）

Figure 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，θ = 0°, 45°, 90°, 135°时流场RMS速度 ${u''}$ 、 ${v''}$ 以及湍动能k在不同λ下沿x方向的分布

Figure 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

下载: 全尺寸图片幻灯片

图 9 t = 3.5时不同θ和λ下流场内能、流场动能和流场湍动能在计算域上游区域（ $x \in \left[ { - 3.0, - 0.5} \right]$ ）、椭圆柱云区域（ $x \in \left[ { - 0.5,0.5} \right]$ ）和计算域下游区域（ $x \in \left[ {0.5, 4.0} \right]$ ）分布

Figure 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 \in \left[ { - 3.0, - 0.5} \right]$ , elliptical column cloud area $x \in \left[ { - 0.5,0.5} \right]$ , the downstream area of the domain $x \in \left[ {0.5, 4.0} \right]$ at dimensionless time t = 3.5

下载: 全尺寸图片幻灯片

图 10 不同角度θ和不同λ下流场内能、流场动能和流场湍动能随无量纲时间t的变化

Figure 10. Variations of the fluid internal energy, fluid kinetic energy and fluid turbulent kinetic energy with dimensionless time t at different θ and λ

下载: 全尺寸图片幻灯片

图 11 t = 3.5时不同的λ和θ下一维体积平均模型与DNS拟合结果

Figure 11. Fitting results of the 1-D volume-averaged model and DNS at different λ and θ when the dimensionless time is equal to 3.5

下载: 全尺寸图片幻灯片

图 12 人工有效阻力系数C_d的最优取值分布

Figure 12. Distribution for the optimal value of artificial effective drag coefficient C_d

下载: 全尺寸图片幻灯片

表 1 网格收敛性分析实验中使用的4种网格

Table 1. Four meshes used in the convergence analysis experiment

Mesh n_b N_x N_y

1 8 112 64
2 16 224 128
3 32 448 256
4 64 896 512

下载: 导出CSV

表 2 平面激波与椭圆柱云相互作用数值模拟使用的网格设置

Table 2. Mesh settings in numerical simulation of the interaction between plane shock and elliptical column cloud

λ N_p a b Δx/10^–4 N_x N_y N/10⁶

2 440 0.029 44 0.014 72 4.59 3 894 2 176 8.5
3 440 0.036 04 0.012 02 3.75 4 763 2 666 12.7
4 440 0.041 64 0.010 41 3.25 5 500 3 078 16.9

下载: 导出CSV

表 3 人工有效阻力系数C_d的最优取值

Table 3. Optimal values of artificial effective drag coefficient C_d

λ C_d
θ = 0° θ = 15° θ = 30° θ = 45° θ = 60° θ = 90° θ = 135°

2 1.19 1.06 1.06 1.42 2.22 6.00 1.68
3 0.86 0.67 0.69 1.48 2.80 14.00 1.29
4 0.57 0.58 0.47 1.30 3.00 36.00 1.15

下载: 导出CSV

参考文献(30)

[1]	罗志全. 核心坍缩型超新星的相关物理过程及爆发机制的研究 [D]. 成都: 四川大学, 2006. LUO Z Q. Research on the physical process and explosion mechanism of core-collapse supernova [D]. Chengdu: Sichuan University, 2006.
[2]	CHOJNICKI K, CLARKE A B, PHILLIPS J C. A shock-tube investigation of the dynamics of gas-particle mixtures: implications for explosive volcanic eruptions [J]. Geophysical Research Letters, 2006, 33(15): 292–306.
[3]	张莉聪, 徐景德, 吴兵, 等. 甲烷-煤尘爆炸波与障碍物相互作用的数值研究 [J]. 中国安全科学学报, 2004(8): 85–88. ZHANG L C, XU J D, WU B, et al. Study on numerical value of reaction between barrier and explosion wave of methane-coal dust [J]. China Safety Science Journal, 2004(8): 85–88.
[4]	QUINLAN N J, KENDALL M A F, BELLHOUSE B J, et al. Investigations of gas and particle dynamics in first generation needle-free drug delivery devices [J]. Shock Waves, 2001, 10(6): 395–404. doi: 10.1007/PL00004052
[5]	张晓立, 解立峰, 洪滔, 等. 激波管驱动石英砂颗粒抛洒的数值模拟 [J]. 高压物理学报, 2014, 28(1): 97–102. doi: 10.11858/gywlxb.2014.01.016 ZHANG X L, XIE L F, HONG T, et al. Numerical simulation of quartz sand dispersion under shock tube loading [J]. Chinese Journal of High Pressure Physics, 2014, 28(1): 97–102. doi: 10.11858/gywlxb.2014.01.016
[6]	ZHANG F, FROST D L, THIBAULT P A, et al. Explosive dispersal of solid particles [J]. Shock Waves, 2001, 10(6): 431–443. doi: 10.1007/PL00004050
[7]	RUDINGER G. Fundamentals of gas-particle flow [M]. Elsevier Scientific Publishing Company, 1980.
[8]	REGELE J D, RABINOVITCH J, COLONIUS T, et al. Unsteady effects in dense, high speed, particle laden flows [J]. International Journal of Multiphase Flow, 2014, 61: 1–13. doi: 10.1016/j.ijmultiphaseflow.2013.12.007
[9]	ZAREI Z, FROST D L, TIMOFEEV E V. Numerical modelling of the entrainment of particles in inviscid supersonic flow [J]. Shock Waves, 2011, 21(4): 341–355. doi: 10.1007/s00193-011-0311-5
[10]	JACOBS G B, DON W S, DITTMANN T. High-order resolution Eulerian-Lagrangian simulations of particle dispersion in the accelerated flow behind a moving shock [J]. Theoretical and Computational Fluid Dynamics, 2012, 26(1/2/3/4): 37–50. doi: 10.1007/s00162-010-0214-6
[11]	ROGUE X, RODRIGUEZ G, HAAS J F, et al. Experimental and numerical investigation of the shock-induced fluidization of a particles bed [J]. Shock Waves, 1998, 8(1): 29–45. doi: 10.1007/s001930050096
[12]	WAGNER J L, BERESH S J, KEARNEY S P, et al. A multiphase shock tube for shock wave interactions with dense particle fields [J]. Experiments in Fluids, 2012, 52(6): 1507–1517. doi: 10.1007/s00348-012-1272-x
[13]	WAGNER J L, KEARNEY S P, BERESH S J, et al. Flash X-ray measurements on the shock-induced dispersal of a dense particle curtain [J]. Experiments in Fluids, 2015, 56(12): 213. doi: 10.1007/s00348-015-2087-3
[14]	THEOFANOUS T G, MITKIN V, CHANG C H. The dynamics of dense particle clouds subjected to shock waves. Part 1. experiments and scaling laws [J]. Journal of Fluid Mechanics, 2016, 792: 658–681. doi: 10.1017/jfm.2016.97
[15]	THEOFANOUS T G, CHANG C H. The dynamics of dense particle clouds subjected to shock waves. Part 2. modeling/numerical issues and the way forward [J]. International Journal of Multiphase Flow, 2017, 89: 177–206. doi: 10.1016/j.ijmultiphaseflow.2016.10.004
[16]	WANG L P, ROSA B, GAO H, et al. Turbulent collision of inertial particles: point-particle based, hybrid simulations and beyond [J]. International Journal of Multiphase Flow, 2009, 35(9): 854–867. doi: 10.1016/j.ijmultiphaseflow.2009.02.012
[17]	LING Y, WAGNER J L, BERESH S J, et al. Interaction of a planar shock wave with a dense particle curtain: modeling and experiments [J]. Physics of Fluids, 2012, 24(11): 113301. doi: 10.1063/1.4768815
[18]	HU H H. Direct simulation of flows of solid-liquid mixtures [J]. International Journal of Multiphase Flow, 1996, 22(2): 335–352. doi: 10.1016/0301-9322(95)00068-2
[19]	XIONG Q, LI B, ZHOU G, et al. Large-scale DNS of gas-solid flows on Mole-8.5 [J]. Chemical Engineering Science, 2012, 71: 422–430. doi: 10.1016/j.ces.2011.10.059
[20]	PICANO F, BREUGEM W P, BRANDT L. Turbulent channel flow of dense suspensions of neutrally buoyant spheres [J]. Journal of Fluid Mechanics, 2015, 764: 463–487. doi: 10.1017/jfm.2014.704
[21]	WANG S, VANELLA M, BALARAS E. A hydrodynamic stress model for simulating turbulence/particle interactions with immersed boundary methods [J]. Journal of Computational Physics, 2019, 382: 240–263. doi: 10.1016/j.jcp.2019.01.010
[22]	ZHU C, YU Z, SHAO X. Interface-resolved direct numerical simulations of the interactions between neutrally buoyant spheroidal particles and turbulent channel flows [J]. Physics of Fluids, 2018, 30(11): 115103. doi: 10.1063/1.5051592
[23]	ZASTAWNY M, MALLOUPPAS G, ZHAO F, et al. Derivation of drag and lift force and torque coefficients for non-spherical particles in flows [J]. International Journal of Multiphase Flow, 2012, 39: 227–239. doi: 10.1016/j.ijmultiphaseflow.2011.09.004
[24]	邹立勇, 廖深飞, 刘金宏, 等. 双椭圆界面Richtmyer-Meshkov流动中的相互干扰效应 [J]. 高压物理学报, 2015, 29(3): 191–198. doi: 10.11858/gywlxb.2015.03.005 ZOU L Y, LIAO S F, LIU J H, et al. Interaction effect of two ellipse Richtmyer-Meshkov flows [J]. Chinese Journal of High Pressure Physics, 2015, 29(3): 191–198. doi: 10.11858/gywlxb.2015.03.005
[25]	JIANG L J, DENG X L, TAO L. DNS study of initial-stage shock-particle curtain interaction [J]. Communications in Computational Physics, 2018, 23(4): 1202–1222.
[26]	STEWART H B, WENDROFF B. Two-phase flow: models and methods [J]. Journal of Computational Physics, 1984, 56(3): 363–409. doi: 10.1016/0021-9991(84)90103-7
[27]	CHANG C H, LIOU M S. A robust and accurate approach to computing compressible multiphase flow: stratified flow model and AUSM(+)-up scheme [J]. Journal of Computational Physics, 2008, 227(10): 5360–5360. doi: 10.1016/j.jcp.2008.01.041
[28]	DENG X, JIANG L, DING Y. Direct numerical simulation of long-term shock-particle curtain interaction [C]//2018 AIAA Aerospace Sciences Meeting. Florida: American Institute of Aeronautics and Astronautics, 2018.
[29]	LIOU M S. Ten years in the making-AUSM-family [C]//15th AIAA Computational Fluid Dynamics Conference, 2001: 2521.
[30]	LIOU MS. A sequel to AUSM, Part II: AUSM+-up for all speeds [J]. Journal of Computational Physics, 2006, 214(1): 137–170. doi: 10.1016/j.jcp.2005.09.020