引 言

高超声速飞行器以极高速(Ma > 10)飞行, 将在近壁区产生激波压缩、黏性摩擦、气动加热等作用, 需要设计特定的热防护系统以抵御气动热载荷. 烧蚀防热通过飞行器表面材料反应、热解等将表面热载荷带离, 是飞行器热防护的重要形式之一^[1-2]. 热防护系统常用的纤维编制类、颗粒增强类复合材料(见图1), 在高焓气动载荷下材料颗粒相或纤维束发生力学剥蚀, 对飞行器流态、气动性能、热载荷产生影响^[3-4]. 另外, 依托电弧加热类风洞设备开展材料响应评价实验时, 在极高焓电弧作用下, 材料将发生微结构变形、蠕化、剥离等现象, 形成的烧蚀颗粒进入流场, 污染电弧风洞流场, 进而干扰材料性能评价^[5-6]. 因此, 烧蚀颗粒的剥离动力学问题是高超声速领域的共性基础性科学问题, 研究气动载荷下烧蚀剥离颗粒的输运特性对防热材料体系评价和热防护系统性能评估均具有重要而又紧迫的工程价值.

针对飞行器烧蚀及颗粒剥离问题, Li等^[7]研究了碳化复合材料烧蚀面后退的非线性热解层模型, 并预测了烧蚀响应; 文献[8]利用线性稳定理论研究了防热材料的温度响应及烧蚀气体产物对高温边界层不稳定性的影响; 国义军等^[9]在前人基础上提出了氧化烧蚀的双平台理论和反应控制机理. 上述研究重点关注烧蚀后退量工程计算、烧蚀剥离的力学/热学条件、烧蚀机理等问题, 未见有针对烧蚀颗粒在近壁运动的动力学机理以及对边界层的调制等方面的研究. 另一方面, 作为一个经典的力学问题, 颗粒多相流已在风沙、化工等颗粒流动问题中得到了较为深入的研究, 并形成了颗粒解析的相关数值模拟方法和软件^[10-13]. 文献[14]综述了颗粒、液滴和气泡与湍流的相互作用和调制, 及其相关的实验和模拟技术. Yu和Shao^[15]提出基于虚拟区域方法的颗粒悬浮流直接数值模拟, 并证明了其在各种典型颗粒流动均适用并具有较好的精度和鲁棒性. 文献[16]通过理论分析得出了球形悬浮颗粒在有界流体中的侧向迁移现象. 针对气动载荷下防热材料剥离颗粒动力学问题, 亟需借鉴颗粒解析的模拟方法, 进一步开展气动载荷下防热材料剥离颗粒输运特性研究, 以期揭示烧蚀颗粒在近壁运动的动力学机理以及对边界层的调制规律, 支撑高超声速飞行器热防护系统设计.

本文针对防热材料烧蚀剥离问题, 采用近壁流动量纲分析方法, 将气动环境下材料烧蚀颗粒剥离过程模化为单个圆球惯性烧蚀颗粒在Couette流动中的动力学问题, 采用颗粒解析的直接数值模拟方法(particle resolved-direct numerical simulation, PR-DNS), 开展烧蚀颗粒剥离和输运特性原理性分析, 为烧蚀质量损失的预测和边界层调制研究提供理论基础.

1.   颗粒及流动模型

1.1   典型烧蚀颗粒模型

常用防热烧蚀材料为纤维编织/颗粒填充增强复合材料. 烧蚀过程分为热化学烧蚀和机械剥蚀两部分, 前者指烧蚀材料与高温气体组分之间的化学反应以及烧蚀材料在高温情况下的升华过程^[17-18], 后者指基体和基质密度不同造成材料烧蚀速率差异而暴露出的纤维束、黏结剂优先烧蚀而暴露出的颗粒填充物等, 在气流压力以及剪切力作用下被剥离吹走^[19-20], 如图1所示.

图  1  纤维增强防热材料在高温气流作用下的机械剥落示意图

Figure  1.  Schematic diagram of mechanical spalling of fiber reinforced ablative composite materials under high temperature air flow

下载: 全尺寸图片 幻灯片

单根纤维的直径一般1 ~ 10 μm, 纤维断裂一般是多根一起发生, 产生的断裂产物一般几微米到几十微米; 颗粒填充物反应产物的典型粒径为0.1 ~ 10 μm, 因此, 实际烧蚀颗粒的尺寸在微纳米量级0.1 ~ 10 μm^[21-22]. 当颗粒粒径大于1 μm时, 只需考虑颗粒与流体相互作用力; 当颗粒直径小于1 μm时, 还需要考虑颗粒之间的分子间作用力、颗粒的团聚作用等. 实际烧蚀颗粒不仅粒径尺寸多变, 而且具有不同的长径比, 且由于烧蚀掉的材料不同, 颗粒密度也不同. 本文以真实数据为参考, 研究单个圆球形惯性颗粒, 直径量级$d_{\rm{p}}^*$ ~ ${{O}}$(1~10) μm (本文上标*代表有量纲物理量, 其他为无量纲物理量, 除了角度量θ, α, β).

1.2   典型流动模型

高超声速飞行器飞行过程中, 来流马赫数Ma高达10 ~ 20, 自由来流的空气参数为高空20 ~ 30 km, 这是烧蚀问题发生的典型参数区间^[23-24], 密度$ \rho _1^* $ = 18 g/m³, 温度$ T_1^* $ = 226.5 K, 压强$ p_1^* $ = 1.197 kPa, 声速$ a_1^* $ = 301.7 m/s, 来流速度$ U_1^* $ = 6034 m/s, 根据“Properties of the U.S. Standard Atmosphere 1976”^[25], 在典型烧蚀空间高度范围内, 运动黏性系数$ \nu _1^* $ = 8.3 × 10⁻⁴ m²/s. 飞行器钝头典型直径D^*为0.05 ~ 0.2 m. 经过激波压缩后, 典型空气参数范围: 密度$ \rho _2^* $ ~ O (0.1) kg/m³, $ T_2^* $ ~ O (10³) K, $ U_2^* $ ~ O (10³) m/s, $ \nu _2^* $ ~ O (10⁻⁴) m²/s, 见表1. 图2给出了激波前后流体参数示例, 其中μ为动力黏性系数.

表  1  激波前后来流参数范围

Table  1.  Flow parameters beside the shock

	Before shock		After shock
density/(g·m−3)	$ \rho _1^* $	18	$ \rho _2^* $	O (0.1)
temperature/K	$ T_1^* $	226.5	$ T_2^* $	O (103)
inflow velocity/(m·s−1)	$ U_1^* $	6034	$ U_2^* $	O (103)
coefficient of kinematic viscosity/(m2·s−1)	$ \nu _1^* $	8.3 × 10−4	$ \nu _2^* $	O (10−4)

下载: 导出CSV 

| 显示表格

图  2  飞行器头部流动示意简图

Figure  2.  Flow around the aircraft nose

下载: 全尺寸图片 幻灯片

飞行器头部位置流动可模化为经典的楔形物体绕流问题, 其控制方程为Falkner-Skan方程, 是一个典型的层流边界层相似解问题. Falkner-Skan方程^[26]表述为

边界条件为

式中$ \beta $为压力梯度, 且$ \beta $与等效半锥角$ \alpha $的关系为

根据$ \beta $取不同的值, 飞行器头部流动可模化为三种不同的典型标模流动, 如图2所示.

(1) $ \beta {\text{ = }}1 $, Hiemenz流动

Hiemenz流动区域有一个向下的挤压力, 正压力梯度, 加速流动.

(2) $ 0 < \beta < 1 $, Falkner-Skan流动

Falkner-Skan流动区域无向下挤压力, 正压力梯度, 加速流动.

(3) $ \beta {\text{ = 0}} $, Blasius流动

Blasius流动区域无向下挤压力, 无压力梯度, 平行剪切流, 本文重点研究该区域.

边界层外边界处的势流速度$ U_e^* $为

其中, x^*为沿壁面位移, c为常数.

不同压力梯度的Falkner-Skan流动速度分布见图3, 本文研究的壁面烧蚀颗粒在微米量级$d_{\rm{p}}^*/{\text{m}}\sim {{O}}\left( {{{10}^{ - 6}} \sim {{10}^{ - 5}}} \right)$, 因此不妨先关注距离壁面10 ~ 100倍颗粒直径, 即y^*≤0.1 m范围内颗粒的启动和输运情况.

图  3  不同压力梯度Falkner-Skan流动速度分布图

Figure  3.  Velocity profiles of Falkner-Skan flow at different half cone angles

下载: 全尺寸图片 幻灯片

在$ \beta {\text{ = }}0 $的平行剪切流近壁区y^*≤0.1 m, 相似变量

此时, 剪切率B^*接近常数, 如图3红色虚线所示, 因此, 作为一种原理性研究, 可用剪切率为常数的Couette流动来研究微米尺度颗粒在近壁面的输运特性. 根据定义, 剪切率B^*和速度u^*的关系式为

自相似解求解过程中, 有

其中, $\varphi $^*为流函数, $\eta $为相似变量, 将式(8)代入式(7), 并经过变换得剪切率

其中

式中, D^*为飞行器钝头典型直径. 代入实际参数可以计算出剪切率B^* ~ O (10⁷) s⁻¹. 特征长度${\delta ^{\text{*}}}{\text{ = }} $$\sqrt {{{{{{\nu ^*}} \mathord{\left/ {\vphantom {{{\nu ^*}} B}} \right. } B}}^*}} \sim {{O}}\left( {{{10}^{{{ - }}5}}} \right)$m.

流体推动颗粒运动, 颗粒运动的能量由流体惯性提供, 表征颗粒周围的流体惯性的特征参数为颗粒雷诺数

表征颗粒惯性的特征参数为斯托克斯数

式中, ${\rho _r}$为颗粒/流体密度比. 由式(11)和式(12)可知, 颗粒在流体中的运动由$ R{e_{\text{p}}} $和$ St $两个无量纲参数决定, $ {U_{\text{p}}} \sim g\left( {R{e_{\text{p}}},St} \right) \sim \left( {{\rho _r},{{d_{\text{p}}^*} \mathord{\left/ {\vphantom {{d_{\text{p}}^*} {{\delta ^{\text{*}}}}}} \right. } {{\delta ^{\text{*}}}}}} \right) $. 因此, 在数值实验中, 只需要保证${\rho _r}$和${{{d_{\text{p}}^*}} \mathord{\left/ {\vphantom {{{d_{\text{p}}^*}} \delta^* }} \right. } \delta^* }$与真实物理过程一致, 即可保证数值实验模拟真实物理过程的合理性.

真实情况下, 颗粒直径与特征长度的比为

颗粒密度与气体密度比为

2.   数值方法及验证算例

传统的PP-DNS (point particle-DNS)是对颗粒动力学过程的一种近似, 其忽略了有限尺寸带来的尾流效应, 且需要的流体信息是基于未扰流动信息, 而真实PP-DNS求解器中的流场都是受扰动流场^[27-28].

而PR-DNS是基于第一性原理的数值模拟方法: 由于网格尺寸比颗粒尺寸小至少一个数量级, 使得流体求解器可以获得颗粒周围的小尺度流体动力学信息, 然后根据柯西定理, 对有限尺寸颗粒表面的流体动力学信息进行积分, 从而获得基于拉格朗日坐标下的颗粒瞬时受力, 即

其中, Σ为有限尺寸颗粒控制体的表面积.

本研究采用浙江大学开发的基于并行计算的PR-DNS代码来进行数值实验^[29-30], 研究单个惯性颗粒在流动近壁区的动力学问题, 从而揭示其剥离及输运特性.

2.1   中性颗粒在平面Couette流中的验证

验证算例选取平面Couette流中的中性浮力颗粒. 首先在一个长方体计算域中构造一个数值Couette流动, 如图4所示. 计算域为OX × OY × OZ = 16$ {r_{\text{p}}} $ × 16$ {r_{\text{p}}} $ × 8$ {r_{\text{p}}} $, 球形颗粒半径${r_{\text{p}}} = 1$(代码中所有量均为无量纲量), 则上下壁面距离H = 16, 上壁面速度U_w = 8, $\kappa = {r_{\text{p}}}/H$. 在计算域中构造流动, 参数(详见表2)如下: 特征流体剪切率B = 0.5, 特征流体密度$ \rho _{\text{f}}^{}{\text{ = 1}} $, 特征运动黏度$ \nu {\text{ = 1}} $, 特征长度为$ \delta {\text{ = }}\sqrt {\nu {\text{/}}B} {\text{ = }}\sqrt 2 $, 特征速度$ U{\text{ = }}\sqrt {B\nu } {\text{ = }}\sqrt {0.5} $, 特征时间$t = \delta /U{\text{ = }}2$, 来与理论解对比, 验证PR-DNS的准确性和合理性.

图  4  长方体计算域及其相关参数示意图

Figure  4.  Cuboid computing domain and parameters

下载: 全尺寸图片 幻灯片

表  2  数值计算参数设置

Table  2.  Numerical calculation parameters settings

	Dimensional parameters (superscript *)	Reference quantity (subscript ∞)	Dimensionless parameters	Dimensionless values
				validation	numerical experiment
shear rate	${B^*}$	${B_\infty }$	$ B = {{{B^*}} \mathord{\left/ {\vphantom {{{B^*}} {{B_\infty }}}} \right. } {{B_\infty }}} $	0.5	1
fluid density	$ \rho _{\text{f}}^{\text{*}} $	$ \rho _\infty ^{} $	$ \rho _{\text{f}}^{}{\text{ = }}{{\rho _{\text{f}}^{\text{*}}} \mathord{\left/ {\vphantom {{\rho _{\text{f}}^{\text{*}}} {\rho _\infty ^{}}}} \right. } {\rho _\infty ^{}}} $	1	1
particle density	$ \rho _{\text{p}}^{\text{*}} $	$ \rho _\infty ^{} $	$ \rho _{\text{p}}^{}{\text{ = }}{{\rho _{\text{p}}^{\text{*}}} \mathord{\left/ {\vphantom {{\rho _{\text{p}}^{\text{*}}} {\rho _\infty ^{}}}} \right. } {\rho _\infty ^{}}} $	1	10000 20000 30000
coefficient of kinematic viscosity	$ {\nu ^{\text{*}}} $	$ {\nu _\infty } $	$ \nu {\text{ = }}{{{\nu ^{\text{*}}}} \mathord{\left/ {\vphantom {{{\nu ^{\text{*}}}} {{\nu _\infty }}}} \right. } {{\nu _\infty }}} $	1	1
characteristic length	$ {\delta ^{\text{*}}} $	$ {\delta _\infty }{\text{ = }}\sqrt {{\nu _\infty }{\text{/}}{B_\infty }} $	$ \delta {\text{ = }}\sqrt {\nu {\text{/}}B} $	$\sqrt 2 $	1
characteristic velocity	${U^*}$	${U_\infty } = \sqrt {{B_\infty }{\nu _\infty }} $	$U = \sqrt {B\nu } $	$\sqrt {0.5} $	1
characteristic time	${t^*}$	${t_\infty } = {{{\delta _\infty }} \mathord{\left/ {\vphantom {{{\delta _\infty }} {{U_\infty }}}} \right. } {{U_\infty }}}$	$t = {\delta \mathord{\left/ {\vphantom {\delta U}} \right. } U}$	2	1
particle diameter	$d_{_{\text{p}}}^*$	$ {\delta _\infty }{\text{ = }}\sqrt {{\nu _\infty }{\text{/}}{B_\infty }} $	${{{{d_{\text{p}}} = d_{\text{p}}^*} \mathord{\left/ {\vphantom {{{d_{\text{p}}} = d_{\text{p}}^*} \delta }} \right. } \delta }_\infty }$	2	0.5, 0.75, 1, 1.5, 2, 4, 6, 8, 10, 12, 14

下载: 导出CSV 

| 显示表格

初始流场全场施加一个平面Couette流$ U = By,\;\;V = 0,\;\;W = 0 $, 之后每个时刻, OY施加Dirichlet边界条件, OX和$ OZ $方向边界施加周期边界条件, 均匀网格尺寸为$ \varDelta {\text{ = }}{r_{\text{p}}}{\text{/16}} $. 已验证结果与网格尺寸$ \varDelta $、时间步长dt、计算域大小OX × OY × OZ无关.

计算结果图5表明, PR-DNS数值方法预测的中性颗粒在Couette流中的法向速度与文献[31-32]理论解、Wang^[33]数值模拟结果吻合良好.

图  5  Couette流中单个中性颗粒法向速度随颗粒到壁面位移的变化

Figure  5.  The normal velocity of a single neutral particle in Couette flow varies with the particle displacement to the wall

下载: 全尺寸图片 幻灯片

2.2   颗粒解析的直接数值模拟方法

构造一个数值Couette流动, 如表2所示, 其特征流体剪切率B = 1, 特征流体密度$ \rho _{\text{f}}^{}{\text{ = 1}} $, 特征运动黏度$ \nu {\text{ = 1}} $, 则特征长度为$ \delta {\text{ = }}\sqrt {\nu {\text{/}}B} {\text{ = }}1 $, 特征速度$ U{\text{ = }}\sqrt {B\nu } {\text{ = 1}} $, 特征时间$t = \delta /U{\text{ = }}1$. 为保证$ R{e_{\text{p}}} $和St与真实物理过程一致, 数值实验中的${\rho _r} = {\rho _{\text{p}}}/{\rho _{\text{f}}}\sim {{O}}\left( {{{10}^4}} \right)$, $ {d_{\text{p}}}/\delta \sim {O} \left( 1 \right) $, 因此, ${\rho _r}$取值分别为10000, 20000, 30000, $ {d_{\text{p}}}/\delta $取值分别为0.5, 0.75, 1, 1.5, 2, 4, 6, 8, 10, 12, 14.

本文考虑数值计算域是一个长方体, 如上述图4所示, 计算域OX × OY × OZ=10r_p × 10r_p × 5r_p, $ OX \times OY \times OZ{\text{ = 20}}{r_{\text{p}}} \times {\text{20}}{r_{\text{p}}} \times {\text{10}}{r_{\text{p}}} $, 其中$ {r_{\text{p}}} $为无量纲颗粒半径. 网格尺寸为$\varDelta \text{ = }{r}_{\text{p}}/6.4,\; \varDelta \text{ = }{r}_{\text{p}}/12.8$, 时间步长为$\mathrm{d}t\text{ = 0}{.002\;5}t, \text{ 0}\text{.005}t$. 已验证PR-DNS的数值模拟结果与计算域、网格尺寸、时间步长无关, 其中一组网格无关验证如图6所示, 本文使用均匀网格. 初始流场全场施加一个平面Couette流$ U = By{,^{}}V = 0, W = 0 $, 之后每个时刻, OY施加Dirichlet边界条件, OX和$ OZ $方向边界施加周期边界条件.

图  6  网格无关性验证算例(ρ_r = 10000, r_p = 0.25, $ OX \times OY \times OZ{\text{ = 5}}{r_{\text{p}}} \times {\text{5}}{r_{\text{p}}} \times 2.{\text{5}}{r_{\text{p}}} $)

Figure  6.  Grid independent verification(ρ_r = 10000, r_p = 0.25, $ OX \times OY \times OZ{\text{ = 5}}{r_{\text{p}}} \times {\text{5}}{r_{\text{p}}} \times 2.{\text{5}}{r_{\text{p}}} $)

下载: 全尺寸图片 幻灯片

3.   颗粒剥离和输运特性模拟结果分析

3.1   密度比对颗粒输运的影响

相同直径不同密度比颗粒的输运轨迹如图7所示. 从图7可以看出, 颗粒水平位移约是法向位移的20倍, 水平位移远大于法向位移, 即大密度比球形颗粒在Couette流动近壁区以水平输运为主. 除此之外, 密度比越大的颗粒, 法向位移越大, 即越重的颗粒被吹的越高. 为了揭示个中原因, 进一步研究了颗粒输运速度的变化规律.

图  7  不同密度比颗粒输运轨迹 (d_p/δ = 1)

Figure  7.  Particle transport trajectories with different density ratios (d_p/δ = 1)

下载: 全尺寸图片 幻灯片

相同直径不同密度比颗粒的水平输运速度沿流向变化规律如图8所示. 从图8可以看出, 随着密度比$ {\rho _r} $的增大, 颗粒水平输运速度${U_{\rm{p}}}$呈现减小趋势. 而相同直径不同密度比颗粒水平滑移速度沿流向变化规律如图9所示. 由图9可知, 随着密度比$ {\rho _r} $的增大, 颗粒/流体之间的水平滑移速度$\left| {{U_{\text{f}}}{{ - }}{U_{\rm{p}}}} \right|{\text{ = }}\left| {B{y_p}{{ - }}{U_{\rm{p}}}} \right|$呈增加趋势. 图8和图9共同佐证了越重的颗粒越难以被吹远, 这是符合认知规律的, 根据St ~ $ {\rho _r} $(式(12)), 在相同体积下, 密度比$ {\rho _r} $越大, 颗粒惯性St也越大, 沿流向越难以被吹远.

图  8  不同密度比颗粒水平速度沿流向变化规律 (d_p/δ = 1)

Figure  8.  The horizontal velocity of particles with different density ratios varies along the flow direction (d_p/δ = 1)

下载: 全尺寸图片 幻灯片

图  9  不同密度比颗粒水平滑移速度沿流向变化规律 (d_p/δ = 1)

Figure  9.  The horizontal slip velocity of particles with different density ratios varies along the flow direction (d_p/δ = 1)

下载: 全尺寸图片 幻灯片

进一步考察法向运动情况. 相同直径不同密度比颗粒的法向速度随流向变化规律如图10所示. 由图可知, 颗粒密度ρ_r越大, 颗粒法向速度越小. 相同直径不同密度比颗粒法向位移随时间变化曲线如图11所示. 由图可知, 颗粒密度ρ_r越大, 颗粒法向位移越小. 上述现象说明相同体积下, 越重的颗粒, 其惯性St越大, 沿法向也越难以被抬高.

图  10  不同密度比颗粒法向速度沿流向变化规律(d_p/δ = 1)

Figure  10.  Normal velocity of particles with different density ratios varies along the flow direction (d_p/δ = 1)

下载: 全尺寸图片 幻灯片

图  11  不同密度比颗粒法向位移随时间变化规律(d_p/δ = 1)

Figure  11.  The normal displacement of particles with different density ratios varies with time (d_p/δ = 1)

下载: 全尺寸图片 幻灯片

然而, 有趣的是从上述图7可观察到大密度比的颗粒输运轨迹反而在较高的位置, 这可能是由于密度比对颗粒水平输运速度的影响远大于对法向速度的影响, 即密度比的增加对颗粒水平输运速度的抑制作用远大于对法向速度的抑制作用.

3.2   直径对颗粒输运的影响

由于颗粒直径d_p越大, 则颗粒惯性St ( ~ d_p)越大, 而大颗粒的重心离开壁面的距离也越高, 其重心处流体速度越大($ {U_{\text{f}}}{\text{ = }}B{y_{\text{p}}} $), 颗粒受到流体的推力也越大. 因而, 直径对颗粒输运的影响必然是颗粒惯性和流体惯性竞争的结果. 相同密度比、不同直径的颗粒水平输运速度如图12所示. 从图12可以看出, 直径d_p越大的颗粒, 水平输运速度$ {U_{\text{p}}} $越小, 体现了颗粒惯性作用优势.

图  12  不同直径颗粒水平速度沿流向变化规律($ {\rho _r}{\text{ = }}10\;000 $)

Figure  12.  The horizontal velocity of particles with different diameters varies along the flow direction ($ {\rho _r}{\text{ = }}10\;000 $)

下载: 全尺寸图片 幻灯片

然而, 相同密度比、不同直径的颗粒水平滑移速度沿流向变化规律, 如图13所示, 越大的颗粒, 其周围流体速度$ {U_{\text{f}}}{\text{ = }}B{y_{\text{p}}} $也越大, 这导致颗粒/流体之间的滑移速度$\left| {{U_{\text{f}}}{{ - }}{U_{\rm{p}}}} \right|$就越大, 体现出流体惯性作用优势.

图  13  不同直径颗粒水平滑移速度沿流向变化规律(${\rho _r}{\text{ = }}10\;000$)

Figure  13.  The horizontal slip velocity of particles with different diameters varies along the flow direction ($ {\rho _r}{\text{ = }}10\;000 $)

下载: 全尺寸图片 幻灯片

图14为不同直径颗粒法向速度沿流向变化规律, 表明直径越大的颗粒, 法向速度越大, 即越大的颗粒, 越容易被吹高. 图15为不同直径颗粒的输运轨迹, 从中看出, 颗粒直径越大, 法向位移越大. 这是由于剪切流中的颗粒会受到向上的Saffman升力$ {F_L}^{\text{*}} $^[34-35]

图  14  不同直径颗粒法向速度沿流向变化规律($ {\rho _r}{\text{ = }}10\;000 $)

Figure  14.  The normal velocity of particles with different diameters varies along the flow direction ($ {\rho _r}{\text{ = }}10\;000 $)

下载: 全尺寸图片 幻灯片

图  15  不同直径颗粒的运动轨迹 ($ {\rho _r}{\text{ = }}10\;000 $)

Figure  15.  Trajectories of particles of different diameters ($ {\rho _r}{\text{ = }}10\;000 $)

下载: 全尺寸图片 幻灯片

其中$R{e_G}$是剪切雷诺数. 从颗粒所受升力表达式中可以看出$F_{_L}^{\text{*}} \sim d_{\text{p}}^{\text{*}}$, 因此, 颗粒直径$d_{\text{p}}^{\text{*}}$越大, 颗粒所受升力$F_{_L}^{\text{*}}$越大. 文献[36]也给出了单个颗粒在平行剪切流中由Lamb矢量诱导的垂直于流动方向的升力: $ F_{_L}^{\text{*}}{\sim}U_{_{{\text{slip}}}}^{\text{*}} \times \omega _{_{{\text{slip}}}}^{\text{*}} $, 因此, 颗粒和流体之间的滑移速度越大(见图13), 颗粒受到升力$F_{_L}^{\text{*}}$越大, 则颗粒法向速度和法向位移越大.

由图15还可以看出, 相同密度比不同直径的颗粒, 水平位移是法向位移的10 ~ 20倍, 颗粒也是以水平输运为主. 图16进一步给出了受颗粒影响的流场压力云图, 由图16可知, 水平方向压差明显大于法向压差, 颗粒主要在较大的水平压差作用下运动, 以水平输运为主.

图  16  含颗粒的Couette流场压力云图

Figure  16.  Pressure contour of particle influenced flow field

下载: 全尺寸图片 幻灯片

3.3   归一化的颗粒输运规律

为了进一步找出颗粒输运的统一规律, 分别选取颗粒水平位移达到$ {X_{\text{p}}}{{ - }}{X_{{\text{p0}}}}{\text{ = }}{r_{\text{p}}}{,^{}}2{r_{\text{p}}}{,^{}}3{r_{\text{p}}} $为颗粒启动的阈值, 认为颗粒运动到启动阈值时的速度为颗粒启动速度U_p0, 即烧蚀颗粒从壁面的剥离速度. 研究发现无量纲颗粒启动速度${{\left( {B{r_{_{\rm{p}}}}{{ - }}{U_{{\text{p0}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {B{r_{_p}}{\text{ - }}{U_{{\text{p0}}}}} \right)} {S{t^{0.01}}}}} \right. } {S{t^{0.01}}}}$对阈值的选取不敏感. 因此, 本文不妨选取$ {X_{\text{p}}}{{ - }}{X_{{\text{p0}}}}{\text{ = }}{r_{\text{p}}} $为阈值, 获得所对应的颗粒启动速度, 进而可得颗粒启动速度随颗粒雷诺数Re_p的变化规律, 如图17所示.

图  17  不同密度比颗粒的非定常归一化启动速度随颗粒雷诺数变化规律

Figure  17.  Dimensionless unsteady start-up velocity of particles along the particle Reynolds number for different density ratios

下载: 全尺寸图片 幻灯片

由图17可知, 无量纲的颗粒启动速度随着颗粒雷诺数Re_p增加而增加. 通过拟合, 进一步得到颗粒启动速度的归一化表达式

变形后, 得

其中, 无量纲剪切率$ B{\text{ = }}{{{B^*}} \mathord{\left/ {\vphantom {{{B^*}} {{B_\infty }}}} \right. } {{B_\infty }}}{\text{ = }}1 $. 式(19)清晰地表明了颗粒和流体惯性作用对颗粒输运的影响. 右边第一项为流体微团或中性颗粒作为示踪粒子在流体中的启动速度^[37-40]; 右边第二项为惯性修正项, 体现了流体惯性和颗粒惯性的耦合修正效应, ${U_{{\text{p0}}}} = g\left( {R{e_{\text{p}}},St} \right)$, 其中颗粒惯性St对颗粒启动速度U_p0起到了阻滞作用.

4.   结论

本文通过单个球形惯性颗粒在Couette流动的PR-DNS模拟, 获得了不同颗粒/流体密度比、不同粒径条件下烧蚀颗粒在近壁流动中的输运规律.

(1)获得了烧蚀颗粒关键特征参量对颗粒输运的动力学规律: 输运速度与颗粒/流体密度比ρ_r和颗粒直径d_p有关. 烧蚀颗粒密度比ρ_r增大, 输运速度因颗粒惯性St增大而减小; 颗粒粒径d_p增大, 水平输运速度因颗粒惯性St增大而减小, 法向速度和位移均因受到的Saffman升力的增大而增大. 烧蚀颗粒法向位移远小于水平位移, 颗粒以水平输运为主.

(2) PR-DNS得到烧蚀颗粒启动速度的归一化表达式, 与理论分析一致, 归一化的颗粒启动速度是颗粒和流体惯性的函数${U_{{\text{p0}}}} = g\left( {R{e_{\text{p}}},St} \right)$, 这与前人的经典理论和数值模拟结果定性吻合^{[33, 35, 37]}. 具体地, 颗粒水平输运速度等于流体微团或中性颗粒的速度减去颗粒和流体的惯性修正项.

本文研究对烧蚀颗粒动力学问题作了单个颗粒、球形颗粒、Couette流动等假设, 未考虑颗粒间相互作用、复杂颗粒形貌、近壁真实流动条件下的颗粒动力学机理. 后续研究将进一步重点关注可压缩性、颗粒群作用、椭球形等复杂颗粒形状、真实流动等更符合实际物理过程的烧蚀颗粒动力学问题.

致谢

感谢浙江大学航空航天学院余钊圣教授提供Fortran并行 PR-DNS代码支持.