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

吴博文 章利特 余秋李 刘天程 冯子龙

摘 要：为研究激波透过大孔隙率结构化球阵的气体流动特性，采用Fluent软件三维雷诺平均Navier-Stokes方程和k-epsilon湍流模型，对不同颗粒（或固相）体积分数和入射马赫数时激波透过面心立方球阵的气相流场进行了数值模拟，重点分析了球阵内部的流动特性和反射波强度的参数影响。研究结果表明，不同颗粒体积分数和入射马赫数下FCC球阵对入射激波的反射效应有较大不同。随着入射马赫数的增大，反射效应不断增强，当马赫数增大到某一确定值后，反射效应则趋于稳定。在颗粒体积分数为5%时，此确定马赫数值应介于1.9与2.3之间。而随着体积分数增大，球阵对入射激波的反射效应先增强后基本趋于稳定。当入射马赫数固定为1.4时，体积分数为7%的FCC球阵具有最强的反射效应。

关键词：FCC球阵;激波;反射波;体积分数

中图分类号：TH311 文献标识码：A 文章编号：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一个区域，将此区域温度、压力和速度均设置为波后静温、静压和波后气流速度。