0 引　言

圆球绕流曳力系数（drag coefficient，C_D）是流体力学研究的重要基础问题，其可作为量化指标广泛应用于航空航天、医药加工、石油化工及能源工业等领域中。

在通常的不可压缩连续流体中，圆球绕流曳力系数C_D被认为是雷诺数（Reynolds number，Re）的函数，一百多年来的国内外学者也对此进行了详尽研究。而在飞行器涉及的高速高空环境中，往往伴有高马赫数的可压缩性以及高努森数条件下的稀薄效应，C_D不仅与Re相关，还与马赫数（Mach number，Ma）或者努森数（Knudsen number，Kn）相关。

高马赫数条件下圆球绕流曳力系数C_D的研究主要有地面实验^[1-13]和数值模拟^[14-19]等方法。实验方法又可分为直接测量和间接计算两种方式。直接测量是指在风洞实验中通过测量细线悬挂球体偏转角度或金属支撑球体受力计算球体的C_D^{[1, 2, 4, 8-11, 13]}。1961年Peter等^[10]在干燥的空气中用细线悬挂球体，测量偏转角度来计算圆球受流体作用力，获得了3.8 < Ma < 4.3，50 < Re < 1000范围内圆球绕流的 C_D值。1962年Aroesty^[9]采用相似的方法在伯克利低密度风洞中测量得到了2 < Ma < 6，10 < Re < 10000范围内的 C_D。间接计算是指通过弹道的形式发射物体，获得位置与时间的关系，从而求解获得C_D^{[3, 5-7, 12]}。1967年Short等^[7]通过弹道实验测量获得了0.4×10⁶ < Re < 1.6×10 ⁶，0.4 < Ma < 14.5范围内圆球绕流的 C_D。1968年Crowe等^[5]通过弹道实验获得了0.01 < Re < 5.1，0.036 < Ma < 1.77范围内圆球绕流的 C_D。1972年Bailey^[3]发射泡沫塑料球至干燥空气中，记录不同时刻球的位置，计算得到球体的C_D。

风洞及弹道实验获得的C_D结果较为可信，但其成本高昂，测量过程易受其他因素影响，误差较大且难以获得流场的细节特征。随着计算机技术和流体力学计算方法的发展，数值模拟作为一种可信的方法被用来研究圆球绕流。

Tao等^[14]用格子玻尔兹曼方法（lattice Boltzmann method, LBM）模拟研究了0.1 < Re < 3.5和0.1 < Kn < 1.1条件下的圆球绕流。Nagata等 ^[15]采用计算流体力学（computational fluid dynamics，CFD）方法，模拟研究圆柱区域中50 < Re < 300和0.3 < Ma < 2条件下的圆球绕流，并获得了丰富的绕流尾迹特征。Dogra等 ^[19]使用直接数值模拟蒙特卡洛方法（direct simulation Monte Carlo, DSMC）研究了Re > 100和 Kn < 0.1条件下的圆球绕流。

我们整理了文献中不同Re和Ma下的C_D数据，绘制于图1中。可以发现，文献中C_D数据较为分散，不同文献中的变化规律有时还有所差异；另外，在1≤Re≤20和0.1≤Ma≤3这一趋势变化剧烈的范围内，特别在1≤Re≤10和0.5≤Ma≤3高马赫数低雷诺数区间内，文献中C_D数据不全，难以准确揭示该范围内的C_D与马赫数的变化规律。

在高马赫数低雷诺数条件下，努森数Kn较大，气体流动为过渡流（0.1≤Kn≤10）或者分子流（Kn≥10）。该状态下实验测量困难，而基于连续流体的CFD方法难以适用，LBM方法的玻尔兹曼方程适合描述Kn≥10的流体，分子动力学（molecular dynamics, MD）方法计算精确但计算需求极大。

葛蔚等^[20-21]提出的拟颗粒模型（pseudo-particle model，PPM）将软球模型的时驱算法结合到硬球模型（hard sphere，HS）中，具有MD计算精确和DSMC计算高效的优点，同时PPM采用时间驱动，提高了计算可扩展性。在PPM中，两个拟颗粒在接触碰撞时距离会稍小于粒子直径，将碰撞时距离的统计平均作为粒子有效直径，此时拟颗粒性质就能与硬球颗粒等效，即拟颗粒模型可以视作一种可变径硬球模型。沈国飞和葛蔚^[22]将HS与 PPM 耦合建立了硬球-拟颗粒方法（HS-PPM）：在各子区域内部采用HS模拟流场，而用PPM模拟边界区域，在保证碰撞精确性的同时解决了传统HS不易扩展的问题，便于大规模并行计算。HS-PPM为高马赫数低雷诺数条件下的圆球绕流研究提供了合适的模拟计算方法。

本文计划应用HS-PPM方法模拟1≤Re≤20和0.1≤Ma≤3条件下圆球绕流过程，获得的曳力系数数据用于初步分析C_D(Re, Ma)。

1 圆球绕流HS-PPM模拟

HS-PPM耦合方法模拟圆球绕流示意于图2中，其中U表示来流气速，L为模拟区域长度，H为高度，W为宽度，L₁为圆球离入口距离。

本文分别将气体粒子的质量、直径作为质量、长度单位，拟颗粒更新时间步长为1个时间单位。

1.1 HS-PPM方法

在HS-PPM耦合方法中，HS用于精确计算粒子碰撞，而PPM有效扩展计算规模。两个粒子是否发生碰撞的判据表述为两粒子间距离|r₁–r₂|小于它们的半径之和以及r₁–r₂与v₁–v₂的点积为负值，则它们将以光滑硬球方式碰撞。

$ {{v}}_{1}'={{v}}_{1}-\frac{(1+e){m}_{2}}{{m}_{1}+{m}_{2}}\frac{({{v}}_{1}-{{v}}_{2})\cdot ({{r}}_{1}-{{r}}_{2})}{{|{{r}}_{1}-{{r}}_{2}|}^{2}}({{r}}_{1}-{{r}}_{2}) $

(1)

$ {{v}}_{2}'={{v}}_{2}+\frac{(1+e){m}_{1}}{{m}_{1}+{m}_{2}}\frac{({{v}}_{1}-{{v}}_{2})\cdot ({{r}}_{1}-{{r}}_{2})}{{|{{r}}_{1}-{{r}}_{2}|}^{2}}({{r}}_{1}-{{r}}_{2}) $

(2)

其中 ${{v}}_{1}' $ 、 ${{v}}_{2}' $ 和v₁、v₂分别是两粒子碰撞后和碰撞前的速度，m₁和m₂是它们的质量，e是拟颗粒的回复系数，一般为1。

根据模拟需要将区域划分为N_x×N_y×N_z个子区域，每个子区域内部采用HS，边界区采用PPM连接各HS模拟区域，示意于图2中。每个进程负责一个子区域内计算，相邻进程间构建一层虚拟缓冲区由消息传递接口（message passing interface，MPI）协议实现信息交换。

1.2 气体-圆球作用

运动的气体粒子与圆球发生接触时其路径会发生变化，涉及的气体-圆球间作用通常采用反弹或漫反射方式处理。反弹模型中入射粒子速度法向分量沿法线反转，切向分量保持。漫反射模型中出射粒子速度遵循给定温度下的偏麦克斯韦分布，与入射速度无关，具有完全热适应能力。在高马赫数流动中，实际气体粒子与圆球的作用同时结合了上述两种方式，引入热适应系数α（0≤α≤1）来描述^[23]，即粒子的出射速度中α部分来自于偏麦克斯韦分布的漫反射，（1–α）部分采用反弹模型。

$ {v}'=\alpha {v}_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{e}}+\left(1-\alpha \right){v}_{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r}}{} $

(3)

1.3 边界设置

如图2所示，气体从左边以速度U均匀流入模拟区域，从右边出去的气体粒子将在左边控制区内应用速度重新标定控温法^[20-21]被赋予初速U重新生成，而上下和前后方向采用周期性边界条件，保证气体粒子数守恒。控制区厚度为1，沿y方向和z方向被分割成4N_y×4N_z个子区。这样处理后，气体的流动长时间演化达到稳态状态。

1.4 曳力系数的计算

气体粒子与圆球之间的每次碰撞均会对圆球产生一个冲量，则在一定时间段内统计这些冲量得到圆球受到气体的平均作用力F_D为：

$ {F}_{\mathrm{D}}=\frac{{P}_{2}-{P}_{1}}{{T}_{2}-{T}_{1}}{} $

(4)

其中P₂–P₁表示T₂–T₁时间内动量改变量的累加值。圆球绕流的曳力系数C_D计算公式为：

$ {C}_{{\rm{D}}}=\frac{{F}_{\mathrm{D}}}{A\rho {U}^{2}/2}=\frac{8{F}_{\mathrm{D}}}{{\rho U}^{2}{\text{π}} {D}^{2}}{} $

(5)

其中D为圆球直径。

在本文涉及的条件下，体系运行1×10⁶时间后流动达到稳定，随后开始统计计算曳力系数，统计持续时间为4×10⁶（Re<15）或1×10⁶（Re≥15）。

2 结果与讨论

圆球绕流HS-PPM模拟得到的C_D不仅与Re和Ma相关，还受到模拟设置的影响，主要包括模拟区域设置和热适应系数α。因此本部分首先考察模拟区域设置和热适应系数等模拟设置对C_D的影响规律，从而获得合理的模拟设置用于1≤Re≤20和0.1≤Ma≤3条件下圆球绕流C_D数据的补充和分析。

2.1 曳力系数模拟计算

以Re = 19.4和Ma = 1.195条件下的圆球绕流模拟为例，说明曳力系数的计算过程和误差分析。在L×W×H = 276.8×276.8×276.8区域内用3.24×10⁶气体粒子，模拟了Re = 19.4和Ma = 1.195条件下的圆球绕流过程，记录不同时间段内的气体粒子与圆球间碰撞信息，用于曳力系数计算，结果如图3所示。从图3中可以看到，模拟计算的曳力系数在1×10⁶时间后趋于定值， ${{\overline C}_{\mathrm{D}}}=2.468$ 。在计算曳力系数的1×10⁶到2×10⁶时间内，气体粒子与圆球之间发生1.42×10⁶次碰撞，则曳力系数的统计误差 $\Bigg({E}_{{C}_{\mathrm{D}}}\propto \dfrac{1}{Ma\sqrt{N}} \Bigg)$ 在0.1%以下。

圆球绕流曳力系数的计算精度还与流场是否充分发展相关，因此比较了三个不同条件下的圆球绕流流场模拟结果，见图4。从中可以发现：随着Re增大，尾流在圆球后部逐渐分离，形成对涡（图4（b））；当Ma > 1时，在圆球前面形成激波，此时需要更宽的模拟区域用于激波面往两侧发展（ 图4（c））。这些发现与文献[22，24]中的结果较为相符。

2.2 模拟区域的影响

圆球绕流HS-PPM模拟在一个有限空间内进行，获得的结果通常会与无限大流场内圆球绕流有差异，因此有必要讨论如何设置模拟区域使偏差在合理范围内。

在一定的Re和Ma条件下，分别改变模拟区域长度、宽度（高度）以及圆球位置，作了一系列圆球绕流HS-PPM模拟，其中长度或宽度范围为5D～40D，圆球位置为长度的0.25处、0.38处和0.5处。

图5表示Re = 19.4和Ma = 1.195条件下，C_D模拟结果随长度或宽度的变化趋势。从图中可以看到，随着长度或宽度的增大，C_D逐渐趋于稳定值C_D∞，且在20D后C_D与C_D∞偏差较小。图5还显示了该条件下C_D实验值处于热适应系数α = 0及α = 1下的两个C_D∞之间，从另一个角度论证了HS-PPM模拟获得圆球绕流C_D的合理性。

进一步模拟研究发现：Ma < 1和 Re < 10时，C_D与模拟区域的相关性较小；Ma < 1和 Re > 10时，随着 Re的增大，圆球绕流逐渐不稳定，需要适当增大区域长度满足尾迹的发展；Ma > 1和 Re < 10时，圆球前方流体出现激波面，此时需要增大区域宽度（高度）用于激波面的发展；Ma > 1和 Re>10时，模拟区域长度和宽度（高度）均需要增大。基于上述认识，我们给出了不同条件下模拟区域的建议设置，列于表1中。并按照建议的区域设置模拟了不同Re和Ma下的圆球绕流，得到的C_D结果与C_D∞偏差在2%以内（见表2），表明有限区域的模拟结果可以近似表达无限区域内圆球绕流。

2.3 热适应系数的影响

在高马赫数下稀薄效应明显，气体-圆球间通常介于无滑移条件和完全滑移条件之间，HS-PPM模拟圆球绕流引入热适应系数α来表达气体-圆球间作用中漫反射成分：α = 0时，镜面反弹滑移；α = 1时，完全漫反射。而α的不同取值会对C_D结果造成影响，佐证于图5中。赵琪和马琳博等^[23-24]已对此进行了初步分析，发现圆球绕流C_D模拟值与热适应系数α正相关。

基于上述发现，将不同α下的C_D与实验结果进行比较，获得两者偏差最小的热适应系数合理值α_r。图6表示Re = 5.1和Ma = 0.4条件下不同α模拟得到不同C_D，可以发现C_D与α的关系呈显著的线性关系，C_D=1.78α+4.16（R² = 0.998），进而得到C_D与实验值相符的α_r = 0.643。类似地，完成了Re = 19.4和Ma = 1.195、Re = 30.2和Ma = 1.896条件下的C_D与α的关系研究，获得对应的α_r，结果列于表3中。采用α_r作了验证性模拟，获得的C_D与实验值偏差在1%以内。

不同Re和Ma条件下的α_r存在一定的差异，表3也表明α_r约在0.65-0.75之间，将α_r统一设置为0.7可以平衡精度和严格验证所需的大批量模拟需求。我们采用α_r = 0.7模拟计算得到其他条件下的C_D，并与相应的实验值比较，结果列于表4中。可以发现，α_r = 0.7时的模拟值与实验值偏差均较小（5%以内），故后续模拟均采用α_r = 0.7。

2.4 1≤Re≤20和0.1≤Ma≤3的曳力系数C_D

在得到合适的区域设置条件和热适应系数后。应用HS-PPM模拟圆球绕流获得了1≤Re≤20，0.1≤Ma≤3范围内的曳力系数C_D，结果绘制于图7中。详细的模拟参数和结果列于附录A的表格中。

图7中的实心点表示HS-PPM模拟计算得到的C_D，空心点表示Ma→0条件下标准阻力系数C_D (Re)。从图7中可以看出，在1≤Re≤20范围内，C_D随着Ma的增加而逐渐下降；Re越大，C_D下降趋势趋缓。

马赫数Ma对C_D(Re)的修正可以表达为下面形式，

$ {C}_{{\rm{D}}}\left(Re,Ma\right)={C}_{{\rm{D}}}\left(Re\right)[a{e}^{-bMa}+(1-a\left)\right] $

(6)

其中[ae^–bMa + (1–a)]为曳力系数修正项，系数a和b与Re相关。

按照公式(6)拟合C_D数据得到了不同Re下的系数a和b（列于表5中），并将系数a和b分别拟合成Re的函数（如图8）。

$ a\left(Re\right)=-0.0238Re+0.845 $

(7)

$ b\left(Re\right)=0.203+6.19/\left(Re+1.18\right){} $

(8)

3 结　论

本文应用HS-PPM方法模拟了1≤Re≤20和0.1≤Ma≤3条件下的圆球绕流过程。在讨论确认合理的模拟区域设置和热适应系数α_r = 0.7后，获得不同Re和Ma下的圆球绕流曳力系数C_D数据，并基于此建立高马赫数低雷诺数条件下马赫数Ma对C_D(Re)的修正模型，

$ {C}_{{\rm{D}}}\left(Re,Ma\right)={C}_{{\rm{D}}}\left(Re\right)[a\left(Re\right){e}^{-b\left(Re\right)Ma}+(1-a\left(Re\right)\left)\right]{} $

其中，

$\begin{split}& a\left(Re\right)=-0.0238Re+0.845\\& b\left(Re\right)=0.203+6.19/ \left(Re+1.18\right) \end{split}{\text{。}}$

高马赫数下圆球绕流流动形态复杂，本文仅对1≤Re≤20和0.1≤Ma≤3范围内的圆球绕流曳力系数作了定量分析，而其他范围内的曳力系数变化趋势仍难以用统一的方式进行描述。此外，高马赫数特别是Ma > 3条件下，激波的存在和发展对圆球绕流流动形态影响巨大，而且会发生强烈的高速热效应，这些都会对曳力系数的计算产生重要影响，需要进一步深入研究。