本节包含以下内容：

12.2.1. 颗粒运动方程
12.2.2. 颗粒湍流扩散
12.2.3. 颗粒运动方程的积分

12.2.1 颗粒运动方程

12.2.1.1 颗粒受力平衡

Ansys Fluent 通过在拉格朗日参考系下对颗粒上的力平衡进行积分，来预测离散相颗粒（或液滴或气泡）的运动轨迹。这种力平衡将颗粒惯性与作用在颗粒上的力相等，可以写成

$m_{p} \frac{d u _{p}}{d t} = m_{p} \frac{u - u _{p}}{τ _{r}} + m_{p} \frac{g ( ρ _{p} - ρ )}{ρ _{p}} + F (12.1)$

其中， $m_{p}$ 表示颗粒质量， $u$ 是流体相速度， $u_{p}$ 是颗粒速度， $ρ$ 是流体密度， $ρ_{p}$ 是颗粒密度， $F$ 是附加力， $m_{p} \frac{u - u _{p}}{τ _{r}}$ 是阻力，而 $τ_{r}$ 是液滴或颗粒的松弛时间 [215]（第1069页），计算公式如下：

$τ_{r} = \frac{ρ _{p} d _{p}^{2}}{18 μ} \frac{24}{C _{d} Re} (12.2)$

此处， $μ$ 表示流体的分子粘度， $d_{p}$ 为颗粒直径，而 Re 则为相对雷诺数，其定义为

$R e \equiv \frac{ρ d _{p} ∣ u _{p} - u ∣}{μ} (12.3)$

12.2.1.2 粒子扭矩平衡

粒子旋转是粒子运动的自然组成部分，对粒子在流体中移动的轨迹有显著影响。对于具有高转动惯量的大而重的粒子，这种影响更为明显。在这种情况下，如果在模拟研究中忽略粒子旋转，得到的粒子轨迹可能与实际路径有显著差异。为了考虑粒子旋转，需要求解粒子角动量的一个额外的常微分方程（ODE）：

$I_{p} \frac{d ω _{p}}{d t} = \frac{ρ _{f}}{2} (\frac{d _{p}}{2})^{5} C_{ω} Ω \cdot Ω = T (12.4)$

其中， $I_{p}$ 表示转动惯量， $ω_{p}$ 表示粒子角速度， $ρ_{f}$ 表示流体密度， $d_{p}$ 表示粒子直径， $C_{ω}$ 表示转动阻力系数， $T$ 表示施加在流体域中粒子上的扭矩，而 $Ω$ 则表示通过以下公式计算的粒子与流体间的相对角速度：

$Ω = \frac{1}{2} \nabla \times u_{f} - ω_{p} (12.5)$

对于一个球形粒子，其转动惯量 $I_{p}$ 的计算公式如下：

$I_{p} = \frac{π}{60} ρ_{p} d_{p}^{5} (12.6)$

从公式12.4（第449页）可以看出，扭矩 $T$ 是粒子惯性和阻力之间达到平衡的结果。

有关如何使用粒子旋转功能的详细信息，请参阅《Fluent用户指南》中的粒子旋转部分。

12.2.1.3 包含重力项

虽然公式12.1（第448页）包含了作用在粒子上的重力，但需要注意的是，在Ansys Fluent中，默认的重力加速度为零。如果您希望包含重力作用，必须记得在操作条件对话框中定义重力矢量的大小和方向。

12.2.1.4 其他力

公式12.1（第448页）在粒子受力平衡中考虑了其他力 $(F)$ ，这些力在特定情况下可能很重要。其中第一个是“虚拟质量”力，即加速粒子周围流体所需的力。这个力可以表示为

$F = C_{v m} m_{p} \frac{ρ}{ρ _{p}} (u_{p} \nabla u - \frac{d u _{p}}{d t}) (12.7)$

其中， $C_{v m}$ 为虚拟质量因子，默认值为 0.5。

此外，由于流体中的压力梯度，还会产生一个额外的力：

$F = m_{p} \frac{ρ}{ρ _{p}} u \nabla u (12.8)$

当流体密度远低于颗粒密度时，虚拟质量和压力梯度力并不重要，正如气体流动中的液固颗粒情况一样 $(ρ / ρ_{p} ≪ 1)$ 。对于接近1的 $ρ / ρ_{p}$ 值，虚拟质量和压力梯度力变得显著，建议当密度比大于0.1时包括这些力。有关如何在模型中包含这些力的详细信息，请参阅《Fluent用户指南》中的“在颗粒上包括虚拟质量力和压力梯度效应”。

12.2.1.5 移动参考系中的力

方程12.1（第448页）中的附加力项 $F$ 还包括由于参考系旋转而产生的颗粒上的力。这些力产生于您在移动参考系中建模流动时（参见第21页的“移动参考系中的流动”）。例如，对于围绕 $z$ 轴的旋转，颗粒在笛卡尔 $x$ 和 $y$ 方向上的力可以写为

$m_{p} (1 - \frac{ρ}{ρ _{p}}) Ω^{2} x + 2 m_{p} Ω (u_{p, y} - \frac{ρ}{ρ _{p}} u_{y}) (12.9)$

其中， $u_{p, y}$ 和 $u_{y}$ 分别表示在笛卡尔坐标系 $y$ 方向上的颗粒和流体速度， $Ω$ 表示转速（RPM）。

$m_{p} (1 - \frac{ρ}{ρ _{p}}) Ω^{2} y - 2 m_{p} Ω (u_{p, x} - \frac{ρ}{ρ _{p}} u_{x}) (12.10)$

其中， $u_{p, x}$ 和 $u_{x}$ 分别是粒子和流体在笛卡尔坐标系 $x$ 方向上的速度。

12.2.1.6 热泳力

悬浮在具有温度梯度的气体中的小颗粒会受到一个与梯度方向相反的力。这种现象被称为热泳。Ansys Fluent 可以选择性地在方程 12.1（第 448 页）中的附加力 $F$ 中包含对颗粒的热泳效应：

$F = - D_{T, p} \frac{1}{T} \nabla T (12.11)$

其中， $D_{T, p}$ 表示热泳系数。您可以将该系数定义为常数、多项式或用户自定义函数，也可以采用 Talbot [644]（第1094页）提出的形式：

$D_{T, p} = \frac{6 π d _{p} μ ^{2} C _{s} ( K + C _{t} K n )}{ρ ( 1 + 3 C _{m} K n ) ( 1 + 2 K + 2 C _{t} K n )} (12.12)$

其中：

$K n =$ Knudsen number $= 2 λ / d_{p}$

$λ =$ 流体的平均自由程

$K = k / k_{p}$

$k =$ 仅基于平移能量的流体热导率 $= (15 /4) μ R$

$k_{p} =$ 颗粒热导率

$C_{S} = 1.17$

$C_{t} = 2.18$

$C_{m} = 1.14$

$T =$ 当地流体温度

$μ =$ 流体粘度

此表达式假设颗粒为球形且流体为理想气体。

12.2.1.7 布朗力

对于亚微米颗粒，可以选择在附加力项中包括布朗运动的影响。布朗力的分量被建模为一个高斯白噪声过程，其频谱强度 $S_{n, ij}$ 由 [363] (第 1078 页) 给出：

$S_{n, ij} = S_{0} δ_{ij} (12.13)$

其中， $δ_{ij}$ 是克罗内克δ函数，并且

$S_{0} = \frac{216 v k _{B} T}{π ^{2} ρ d _{p}^{5} ( \frac{ρ _{p}}{ρ} ) ^{2} C _{c}} (12.14)$

其中， $T$ 表示流体的绝对温度， $v$ 是运动粘度， $C_{c}$ 是坎宁安修正系数（定义见公式12.55，第459页），而 $k_{B}$ 则是玻尔兹曼常数。布朗力分量的振幅形式为

$F_{b_{i}} = m_{p} ζ_{i} \frac{π S _{o}}{Δ t} (12.15)$

其中， $ζ_{i}$ 是均值为零、单位方差独立的正态随机数。布朗力分量的幅值在每个时间步长进行评估。为了使布朗力生效，必须启用能量方程。布朗力仅适用于层流模拟。

12.2.1.8 Saffman 升力

Saffman 升力，或因剪切产生的升力，也可以作为选项包含在附加力项中。所使用的升力公式来自 Li 和 Ahmadi [363]（第 1078 页），是 Saffman [562]（第 1089 页）提供的表达式的一般化：

$F = m_{p} \frac{2 K v ^{1/2} ρ d _{ij}}{ρ _{p} d _{p} ( d _{l k} d _{k l} ) ^{1/4}} (u - u_{p}) (12.16)$

其中， $K = 2.594$ ，而 $d_{ij}$ 表示变形张量。这种升力形式适用于小颗粒雷诺数，仅推荐用于亚微米级颗粒。

12.2.1.9 马格努斯升力

当颗粒在流体中旋转时，会产生马格努斯或旋转升力。这种升力是由颗粒表面上的压力差引起的（[124]（第1064页））。对于高雷诺数情况，马格努斯力 $F_{R L}$ 通过旋转升力系数 $C_{R L}$ 进行缩放：

$F_{R L} = \frac{1}{2} A_{p} C_{R L} ρ_{f} \frac{V}{Ω} (V \times Ω) (12.17)$

其中，

$A_{p} =$ projected particle surface area

$V =$ relative fluid - particle velocity

$Ω =$ relative fluid-particle angular velocity

对于旋转升力系数 $C_{R L}$ ，文献中存在多种不同的方法。Ansys Fluent 提供了以下几种计算公式：

Oesterle and Bui Dinh（[484] ）

旋转升力系数 $C_{R L}$ 取决于旋转雷诺数 $R e_{ω}$ 和颗粒雷诺数 $R e_{p}$ ：

$C_{R L} = 0.45 + (4 \frac{R e _{ω}}{R e _{p}} - 0.45) exp (- 0.09896 R e_{ω}^{0.4} R e_{p}^{0.3}) (12.18)$

作者提供的相关性通过其实验测量得到了证实，这些测量针对的是 $R e_{p}$ 值高达140的情况。与其他作者提供的实验数据进行比较，也表明该相关性在 $R e_{p}$ 值高达2000时仍然有效。

Tsuji et al. ([660])

旋转升力系数 $C_{R L}$ 被定义为自旋参数 $S$ 的函数，具体如下：

$C_{R L} C_{R L} = 0.4, = (0.4 \pm 0.1) S, f or S \geq 1 f or S < 1 (12.19)$

自旋参数 $S$ 定义如下：

$S = \frac{ω _{p} d _{p}}{2 u _{f} - u _{p}} (12.20)$

这种表达方式被广泛使用，其有效范围可达 $R e_{p} < 1600$ 。

Rubinow and Keller ([561] (第 1089 页))

在这个基本假设中，旋转升力系数 $C_{R L}$ 与自旋参数 $S$ 成线性正比关系：

$C_{R L} = 2 S (12.21)$

此模型可用于比较目的或学术应用。

对于旋转粒子模拟，您可以在粒子力平衡方程（公式12.1（第448页））的附加力项中包含马格努斯升力。详细信息请参阅《Fluent用户指南》中的“粒子旋转”部分。

12.2.2 粒子的湍流分散

由于流体相中的湍流引起的粒子分散可以使用随机跟踪模型进行预测。随机跟踪（随机游走）模型通过使用随机方法（参见第453页的“随机跟踪”），考虑了瞬时湍流速度波动对粒子轨迹的影响。对于随机跟踪，有一个模型可以考虑连续相中湍流的生成或耗散（参见第583页的“离散相与连续相之间的耦合”）。

重要提示：如果使用Spalart-Allmaras湍流模型，则不能包含粒子的湍流分散。

12.2.2.1 随机跟踪

当流动是湍流时，Ansys Fluent将使用轨迹方程（公式12.1（第448页））中的平均流体相速度 $\overset{u}{ˉ}$ 来预测粒子的轨迹。可选地，您可以包含气体流动速度的瞬时波动值，

$u = \overset{u}{ˉ} + u^{'} (12.22)$

预测颗粒由于湍流引起的分散。

在随机追踪方法中，Ansys Fluent通过积分单个颗粒的轨迹方程，使用瞬时流体速度 $\overset{u}{ˉ} + u^{'} (t)$ ，在积分过程中沿颗粒路径，来预测颗粒的湍流分散。通过以这种方式计算足够数量的代表性颗粒（称为“尝试次数”）的轨迹，可以包括湍流对颗粒分散的随机效应。

Ansys Fluent使用随机方法（随机游走模型）来确定瞬时气体速度。在离散随机游走（DRW）模型中，波动速度分量是时间的离散分段常数函数。它们的随机值在由涡流特征寿命给定的时间间隔内保持不变。

在强烈非均匀扩散主导的流动中，DRW模型可能会给出非物理结果，其中小颗粒应均匀分布。相反，DRW模型会显示这些颗粒倾向于集中在流动中的低湍流区域。此外，由于湍流，该模型已知对直径小于几微米的颗粒的壁面撞击率预测不佳。

12.2.2.1.1 积分时间

颗粒分散的预测利用了积分时间尺度 $T$ 的概念：

$T = \int_{0}^{\infty} \frac{u _{p} ^{'} ( t ) u _{p} ^{'} ( t - τ )}{u _{p} ^{'} ^{2}} d τ (12.23)$

积分时间与颗粒分散速率成正比，较大的值表明流体中存在更剧烈的湍动。可以证明，颗粒的扩散系数由 $\overline{u_{i} u_{j}} T$ 给出。

对于随流体一起运动（零漂移速度）的微小“示踪”颗粒，积分时间变为流体的拉格朗日积分时间， $T_{L}$ 。这一时间尺度可以近似为

$T_{L} = C_{L} \frac{k}{ε} (12.24)$

其中， $C_{L}$ 的值尚未明确，需要进一步确定。通过将示踪粒子的扩散率 $\overline{u_{i}^{'} u_{j}^{'} T_{L}}$ 与湍流模型预测的标量扩散速率 $v_{t} / σ$ 相匹配，我们可以求得 $C_{L}$ 的值。

$T_{L} \approx 0.15 \frac{k}{ε} (12.25)$

对于 $k - ε$ 模型及其变体，以及

$T_{L} \approx 0.30 \frac{k}{ε} (12.26)$

当使用雷诺应力模型（RSM）时 [132]（第1064页）。对于 $k$ - $ω$ 模型，将 $ε$ =0.09kω 代入方程 12.24（第453页）。LES模型使用等效的LES时间尺度。

12.2.2.1.2 离散随机游走模型

在离散随机游走（DRW）模型，或称为“涡流寿命”模型中，粒子与一系列离散的理想化流体相湍流涡的相互作用被模拟 [215]（第1069页）。每个涡流具有以下特征：

高斯分布的随机速度波动， $u^{'}, v^{'}$ 和 $w^{'}$
时间尺度， $τ_{e}$

在湍流涡的寿命期间占主导地位的 $u^{'}, v^{'}$ 和 $w^{'}$ 的值是通过假设它们服从高斯概率分布来采样的，因此

$u^{'} = ζ \overline{u^{'2}} (12.27)$

其中， $ζ$ 是一个正态分布的随机数，而等式右边的其余部分则是速度脉动的局部均方根（RMS）值。由于在流场中每一点处的湍流动能是已知的，因此这些RMS脉动分量值可以被定义（假设各向同性）为：

$\overline{u^{^{'} 2}} = \overline{v^{^{'} 2}} = \overline{w^{^{'} 2}} = 2 k /3 (12.28)$

对于 $k$ - $ε$ 模型、 $k$ - $ω$ 模型及其变体，当使用RSM时，应力非各向同性被纳入速度波动的推导中：

$u^{'} = ζ \overline{u^{'2}} (12.29)$

$v^{'} = ζ \overline{v^{'2}} (12.30)$

$w^{'} = ζ \overline{w^{'2}} (12.31)$

当在湍流的二阶矩为对角线的参考系中观察时 [731]（第1099页）。对于大涡模拟（LES）模型，速度波动在所有方向上都是等效的。详情请参阅尺度解析模拟的入口边界条件（第126页）。

涡流的特征寿命被定义为常数：

$τ_{e} = 2 T_{L} (12.32)$

其中， $T_{L}$ 通常由方程 12.24（第 453 页）给出（默认情况下为方程 12.25，第 454 页），或者作为关于 $T_{L}$ 的随机变化：

$τ_{e} = - T_{L} ln (r) (12.33)$

其中， $r$ 是一个大于零且小于1的均匀随机数，而 $T_{L}$ 由公式12.25（第454页）给出。选择随机计算 $τ_{e}$ 可以更真实地描述相关函数。

粒子涡流穿越时间定义为

$t_{cross} = - τ ln [1 - (\frac{L _{e}}{τ ∣ u - u _{p} ∣})] (12.34)$

其中， $τ$ 是颗粒松弛时间， $L_{e}$ 是涡流长度尺度， $∣ u - u_{p} ∣$ 是相对速度的大小。

假设颗粒在与涡流生命周期和涡流穿越时间中较小的一个时间内与流体相涡流相互作用。当达到这个时间时，通过在方程12.27（第454页）中应用新的 $ζ$ 值，得到瞬时速度的新值。

12.2.2.1.3 使用DRW模型

DRW模型所需的唯一输入是积分时间尺度常数 $C_{L}$ 的值（参见方程12.24（第453页）和方程12.32（第454页））以及用于预测涡流生命周期的方法选择。您可以选择使用常数值或随机值，通过在用户指南中描述的随机跟踪部分，在设置注入属性对话框中为每个注入选择适当的选项。

重要提示：如果使用Spalart-Allmaras湍流模型，则不能包括颗粒的湍流扩散。

12.2.3 颗粒运动方程的积分

轨迹方程以及任何描述颗粒与周围环境之间热量或质量传递的辅助方程，都是通过在离散时间步长上的逐步积分来求解的。时间积分在方程12.1（第448页）中进行，得到沿轨迹各点的颗粒速度，而轨迹本身则通过...

$\frac{d x}{d t} = u_{p} (12.35)$

请注意，方程12.1（第448页）与方程12.35（第455页）构成了一组耦合的常微分方程，而方程12.1（第448页）可以改写成以下通用形式：

$\frac{d u _{p}}{d t} = \frac{1}{τ _{p}} (u - u_{p}) + a (12.36)$

其中，项 $a$ 包含了除阻力外所有其他力引起的加速度。

通过解析积分，我们可以针对常数 $u$ 、 $a$ 和 $τ_{p}$ 求解这一组方程。对于粒子在新位置的速度 $u_{p}^{n + 1}$ ，我们得到

$u_{p}^{n + 1} = u^{n} + e^{- \frac{Δ t}{τ _{p}}} (u_{p}^{n} - u^{n}) - a τ_{p} (e^{- \frac{Δ t}{τ _{p}}} - 1) (12.37)$

新位置 $x_{p}^{n + 1}$ 可以通过类似的关联关系计算得出。

$x_{p}^{n + 1} = x_{p}^{n} + Δ t (u^{n} + a τ_{p}) + τ_{p} (1 - e^{- \frac{Δ t}{τ _{p}}}) (u_{p}^{n} - u^{n} - a τ_{p}) (12.38)$

在这些方程中， $u_{p}^{n}$ 和 $u^{n}$ 分别表示旧位置处的颗粒速度和流体速度。在使用解析离散化方案时，会应用方程 12.37（第 455 页）和方程 12.38（第 456 页）。

方程组 12.1（第 448 页）和方程 12.35（第 455 页）也可以通过数值离散化方案来求解。当对方程 12.36（第 455 页）应用欧拉隐式离散化时，我们得到...

$u_{p}^{n + 1} = \frac{u _{p}^{n} + Δ t ( a + \frac{u ^{n}}{τ _{p}} )}{1 + \frac{Δ t}{τ _{p}}} (12.39)$

在将梯形离散化方法应用于方程12.36（第455页）时，右侧的变量 $u_{p}$ 和 $u_{n}$ 取其平均值，而由于其他力引起的加速度 $a$ 则保持恒定。我们得到

$\frac{u _{p}^{n + 1} - u _{p}^{n}}{Δ t} = \frac{1}{τ _{p}} (u^{*} - u_{p}^{*}) + a^{n} (12.40)$

平均值 $u_{p}^{*}$ 和 $u^{*}$ 是通过计算得出的。

$u_{p}^{*} = \frac{1}{2} (u_{p}^{n} + u_{p}^{n + 1}) (12.41)$

$u^{*} = \frac{1}{2} (u^{n} + u^{n + 1}) (12.42)$

$u^{n + 1} = u^{n} + Δ t u_{p}^{n} \cdot \nabla u^{n} (12.43)$

在新的位置 $n + 1$ 处，粒子速度的计算方法如下：

$u_{p}^{n + 1} = \frac{u _{p}^{n} ( 1 - \frac{1}{2} \frac{Δ t}{τ _{p}} ) + \frac{Δ t}{τ _{p}} ( u ^{n} + \frac{1}{2} Δ t u _{p}^{n} \cdot \nabla u ^{n} ) + Δ t a ^{n}}{1 + \frac{1}{2} \frac{Δ t}{τ _{p}}} (12.44)$

对于隐式和梯形方案，新粒子位置总是通过方程12.35（第455页）的梯形离散化来计算。

$x_{p}^{n + 1} = x_{p}^{n} + \frac{1}{2} Δ t (u_{p}^{n} + u_{p}^{n + 1}) (12.45)$

方程12.36（第455页）和方程12.35（第455页）同样可以通过Cash和Karp[95]（第1062页）发布的Runge-Kutta方案来计算。这些常微分方程可以视为向量形式，其中左侧为导数 $y^{'}$ ，右侧则是任意函数 $f (t, y)$ 。

$y^{'} = f (t, y) (12.46)$

得到

$y^{n + 1} = y^{n} + c_{1} k_{1} + c_{2} k_{2} + c_{3} k_{3} + c_{4} k_{4} + c_{5} k_{5} + c_{6} k_{6} (12.47)$

其中

$k_{1} = Δ t f (t, y^{n}) (12.48)$

$k_{2} = Δ t f (t + a_{2} Δ t, y^{n} + b_{21} k_{1})$

$k_{3} = Δ t f (t + a_{3} Δ t, y^{n} + b_{31} k_{1} + b_{32} k_{2})$

$k_{4} = Δ t f (t + a_{4} Δ t, y^{n} + b_{41} k_{1} + b_{42} k_{2} + b_{43} k_{3})$

$k_{5} = Δ t f (t + a_{5} Δ t, y^{n} + b_{51} k_{1} + b_{52} k_{2} + b_{53} k_{3} + b_{54} k_{4})$

$k_{6} = Δ t f (t + a_{6} Δ t, y^{n} + b_{61} k_{1} + b_{62} k_{2} + b_{63} k_{3} + b_{64} k_{4} + b_{65} k_{5})$

系数 $a_{2} \dots a_{6}$ 、 $b_{21} \dots b_{65}$ 以及 $c_{1} \dots c_{6}$ 取自 Cash and Karp [95]）。

该方案提供了嵌入式误差控制，当未启用Accuracy Control时，该控制将被关闭。

对于移动参考系，积分在移动框架中进行，额外项如公式 12.9（第 450 页）和公式 12.10（第 450 页）所述，因此考虑了系统旋转。利用可用的精度控制机制，轨迹积分将随时间精确进行。

解析方案非常高效。但在大步长或粒子未与连续流体达到流体动力学平衡的情况下，可能会变得不准确。数值方案隐式和梯形，结合自动跟踪方案选择，考虑了作用在粒子上的大部分力的变化，并被选为默认方案。如果非阻力力在粒子积分步长内发生变化，建议使用龙格-库塔方案。

高阶方案梯形和龙格-库塔的积分步长受限于基于粒子动量响应时间的稳定范围。因此，建议将它们与自动跟踪方案选择结合使用。

对于无质量粒子类型，粒子速度等于连续相的速度，因此只需解轨迹方程 12.35（第 455 页），其中粒子速度 $u_{p} = u$ 。沿轨迹的新粒子位置总是通过方程 12.41（第 456 页）和方程 12.45（第 456 页）计算，其中 $u_{p} = u$ 。

旋转粒子运动方程（方程 12.195（第 493 页））采用欧拉隐式离散化方案求解。