激波透过大孔隙率结构化球阵的数值模拟

文章编号：1671-2064（2019）23-0220-04

0 引言

激波与固体颗粒群的相互作用是超声速两相流气体动力学中的一个重要研究课题[1]，激波驱动固体颗粒群理论[2]同时也被广泛运用于军事，医疗，航空航天等领域。迄今为止，国内外已经有很多学者针对激波驱动固体颗粒的流动现象，采用实验和数值模拟的方式进行研究。

Sun等[3]开展了对平面激波透过单球时的非稳态阻力的实验测量，发现雷诺数和单球直径的增大，会造成非稳态阻力的线性递增;Saito等[4]针对含激波携带气固混合物的两相流动开展一系列数值模拟，研究了激波后影响非平衡流场结构的因素主要为非稳态阻力，并发现激波马赫数的降低和颗粒质量比的升高对阻力系数的时间依赖效应的影响更深;Igra等[5]对模型球在非稳态流场中的流动进行了数值模拟，并对模型球进行了受力分析，发现球表面压力和粘性剪切应力决定着球的阻力;Tanno等[6]利用垂直激波管展开了对模型球的实验与数值模拟，发现激波刚接触到模型球时，所受到的非稳态阻力最大，衍射波停滞在模型球后，会在短时间内使阻力值下降至负;Jourdan等[7]利用水平激波管开展了激波驱动面粉颗粒群的流动相关实验，表明入射激波的初始条件和颗粒的物理因素决定颗粒喷射的数量;Shi等[8]通过激波管装置对激波与球阵的相互作用进行实验研究，表明球阵球与球之间存在不可忽视的干扰力，激波驱动单球的实验结果不可直接运用在多球实验中;Zhang等[9]对激波驱动固体颗粒群的流场研究后发现，在激波加载下的颗粒群存在一系列复杂现象，包括激波的透射、衍射、反射以及激波与激波之间的相互扰动现象。

本文对气固两相不同参数下激波与结构化球阵互相作用过程進行了非定常数值模拟，分析了可压缩性气体与固体球阵相互作用的机理，主要包括入射激波马赫数，球阵体积分数等对整体流动的影响，解析了激波遇到不同参数结构化球阵的反射和衍射现象。

1 几何模型

由于真实的固体颗粒群形状大小排列方式各不相同，为方便对不同颗粒群进行研究，把固体颗粒假设成具有周期性的固定空间球阵，从物理及数学上抽象出颗粒群的简化模型，建立其表征体，即球阵模型。为了更好的描述出球阵模型的具体特性，采用自然界常见的面心立方堆积（FCC）排列方式对球阵进行排列，其立方体单元结构阵列如图1所示。

本文中球体固定直径40mm，结合单元体几何特性，选取1%，3%，5%，7%，9%五种体积分数，计算出每种体积分数对应单元体的尺寸，选择五个单元体沿流向方向并排放置，为防止前后单元体球阵的扰动，选择中心单元体进行监测。采用参数化三维造型软件Pro/E对整个FCC球阵的流道进行建模，为了保证流动的充分发展，对FCC球阵进出口段进行了延长。表1展示了各体积分数下FCC球阵的具体几何参数。

2 数值模拟的网格划分及算法

2.1 网格划分

基于对称性，采用ANSYS Workbench软件对1/8个流道进行非结构化四面体网格划分，为了在兼顾计算精度的同时节省计算资源，对流道整体网格尺寸设置为4mm，对球面网格尺寸加密为0.4mm，在稀疏网格过渡区域采用指数增长方式进行过渡。最终得到了总体网格质量在0.6以上的高质量非结构化网格，整体网格和球面局部加密网格见图2。

2.2 控制方程

采用三维非定常雷诺平均Navier-Stokes方程作为气体（空气）流动的基本控制方程，湍流模型选用k-epsilon湍流模型。除连续性方程和能量方程外，动量方程表达形式为：

（1）

其中：U为速度;ρ为液体密度;p为压强;SM为控制方程源项;τ为分子应力张量。

标准k-epsilon湍流模型的方程表达式为：

（2）

（3）

其中：μ为动力粘度，空气粘度采用Sutherland定律进行计算;μt为湍流粘性系数;Gk为平均速度梯度形成的湍动能生成项;Gb是存在浮力的影响形成的湍动能的生成项;ε为湍动能耗散率;k为湍动能;YM是可压速湍流脉动膨胀对总的耗散率的影响;、、为模型参数。在对空气密度进行求解时采用理想气体假设，其状态方程为：

（4）

其中：R为通用气体常数，T为空气温度。

2.3 数值计算方法和边界条件

为了能够更好地捕捉精细的激波结构，选用AUSM差分格式和Roe-fds通量分裂方法，方程采用二阶迎风格式进行离散，并采用隐式算法进行瞬态迭代求解，固定时间步长为5E-7s。在进行数值计算时，工作介质选为空气，采用基于密度的求解器。计算区域的进口边界条件设定为压力入口，参数设置为波后总温和总压，具体数值根据激波管理论得到;计算域的出口边界条件压力出口，压力温度设置为环境温度和环境压力;球体表面设置为固壁无滑移、绝热边界条件。由于采用1/8个流道进行计算，其余三个面均为对称面，设置为对称边界条件。在计算开始前在球阵入口100mm处patch一个区域，将此区域温度、压力和速度均设置为波后静温、静压和波后气流速度。

3 数值模拟结果及分析

3.1 不同入射马赫数下球阵内部流动特性

为了研究不同入射马赫数下球阵内部的流动特性，固定FCC球阵体积分数为5%，分别在入射马赫数为1.4、1.7、1.9、2.3、2.6和2.9的情况下展开非定常数值仿真研究，得到了在这些入射马赫数下激波从球阵入口运动到球阵出口时压力分布和温度分布云图分别如图3和图4所示。

图3从（a）到（e）分别展示了入射马赫数从1.4到3.2时激波经过中心球阵入口、中心和出口时整个FCC球阵的压力分布运动。从图中可以看出，在不同入射马赫数下，随着激波在球阵中的流动，激波强度不断衰减。这是由于在激波流动过程中不断与球阵发生相互作用和激波沿程损失的结果。随着入射马赫数的增大，衰减效果增强，入射激波锋面的最大压力的相对衰减幅度逐渐增大。此外，随着入射马赫数的增大，经过入射激波和球阵之间相互作用，以及单球面对激波的反射作用相互叠加，在入射激波波后区域形成了反射激波面。这种现象在随着入射激波马赫数的增大犹为明显，在马赫数为1.4时，各球面反射的反射波并不能相互作用形成较强的反射激波。而随着马赫数的增大，可以看到在波后区域的逐渐形成较为明显的反射激波区域。在入射马赫数达到1.9及以上时，马赫数增加带来的反射激波增强效应逐渐减弱，可以看出马赫数为2.3，2.6和2.9时的压力分布云图基本是相似的，这一点在温度云图上同样有所反映。因此可以认为在体积分数为5%时，使得球阵能够产生最大反射效应的入射马赫数应在1.9到2.3之间。

从图4的压力云图中可以看出，随着入射马赫數的逐渐增大，波后区域形成的反射波区域温度也随之升高，这同样是由于入射马赫数越高，FCC球阵对入射激波的反射效应更强，导致反射激波的强度更高，反射激波所在区域的温度同样更高。

3.2 不同体积分数下球阵内部流动特性分析

为了研究不同体积分数下球阵内部的流动特性，固定入射马赫数为1.4，分别在球阵体积分数为1%，3%，5%，7%和9%的情况下非定常数值仿真研究，得到了在不同体积分数下激波从中心球阵入口运动到球阵出口时压力分布和温度分布云图分别如图5和图6所示。

从图5中可以看出，在不同体积分数下，入射激波经过球阵和球阵相互作用形成的反射波效应有很大不同。随着体积分数增大，这种反射波效应在增强。这是由于较大体积分数对应球阵长度更短，激波沿程损失小，同时球阵中球体之间距离小，激波经过单球形成的单个反射波更容易叠加形成一道更强大反射激波。在体积分数为1%时，几乎不能形成很明显的反射激波。而随着体积分数的增大，各球面反射的反射波并不能相互作用形成较强的反射激波。而随着马赫数的增大，可以看到在波后区域的逐渐形成较为明显的反射激波区域。在体积分数超过5%时，从压力云图上可以看出此时球阵对入射激波的反射效应随着体积分数增大而带来提升则不是很明显，7%和9%体积分数对应的压力云图基本相似。这说明对于固定入射马赫数而言，存在一个最佳的体积分数，使得球阵对入射激波的反射效应最强。以马赫数为1.4为例，对应的体积分数应在7%左右。

4 结论与展望

本文通过对不同入射马赫数和不同体积分数下入射激波透过FCC球阵这一现象进行非定常数值模拟，得到了以下结论：

（1）不同体积分数和入射马赫数下FCC球阵对入射激波的反射效应有较大不同。（2）随着入射马赫数的增大，反射效应不断增强，当马赫数增大到某一值时，反射效应则保持不变。在体积分数为5%时，此马赫数值应在1.9到2.3之间。（3）随着体积分数增大，球阵对入射激波的反射效应先增强后基本保持不变。在入射马赫数固定为1.4时，体积分数为7%的FCC球阵具有最强的反射效应。

参考文献

[1] 杨亚政，李松年，杨嘉陵.高超音速飞行器及其关键技术简论[J].力学进展，2007（04）：537-550.

[2] M.Sommerfeld.The unsteadiness of shock waves propagating through gas-particle mixtures[J].Experiments in Fluids，1985，3（4）：197-206.

[3] M.Sun，T.Saito，K.Takayama，H.Tanno.Unsteady drag on a sphere by shock wave loading[J].Shock Waves，2005，14（1-2）.

[4] T.Saito，M.Saba，M.Sun，K.Takayama. The effect of an unsteady drag force on the structure of a non-equilibrium region behind a shock wave in a gas-particle mixture[J]. Shock Waves，2007，17（4）：255-262.

[5] Igra O，Takayama K.Shock Tube Study of the Drag Coefficient of a Sphere in a Non-Stationary Flow[J].Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，1993，42（1915）：231-247.

[6] H.Tanno，K.Itoh，T.Saito，A.Abe，K.Takayama.Interaction of a shock with a sphere suspended in a vertical shock tube[J].Shock Waves，2003，13（3）：191-200.

[7] Tanno H，Itoh K，Saito T，et al.Shock wave interaction with a sphere in a shock tube[C].Symposium on Interdisciplinary Shock Wave Research，2004：483-497.

[8] 施红辉.用激波管研究超音速气固两相流[J].应用力学学报，2003，20（04）：41-45.

[9] 章利特，施红辉，王超，等.激波与可运动颗粒群相互作用反射与透射机理实验研究[J].应用力学学报，2010，27（02）：80-86.