在多个工业流体流动应用中，存在具有尺寸分布的次级相。颗粒的尺寸分布，包括固体颗粒、气泡或液滴，可以在多相系统中的传输和化学反应过程中演变。这种演变过程可以是诸如成核、生长、分散、溶解、聚集和破裂等不同现象的组合，从而产生分散体。因此，在涉及尺寸分布的多相流中，除了动量、质量和能量平衡外，还需要一个平衡方程来描述颗粒群的变化。这种平衡通常被称为群体平衡。群体平衡适用的案例包括结晶、气相或液相中的沉淀反应、气泡柱、气体分散、喷雾、流化床聚合、造粒、液-液乳化和分离以及气溶胶流动。

本章节描述了Ansys Fluent中可用的群体平衡模型的理论基础。本章内容组织如下：

14.9.1 简介

为了利用这一建模概念，引入了一个数密度函数来描述颗粒群。借助颗粒属性（例如，颗粒大小、孔隙度、组成等），可以区分群体中的不同颗粒并描述它们的行为。

Ansys Fluent提供了三种求解群体平衡方程的方法：离散群体平衡法、标准矩方法和正交矩方法。

14.9.1.1 离散方法
14.9.1.2 非均匀离散方法
14.9.1.3 标准矩方法
14.9.1.4 正交矩方法

14.9.1.1 离散方法

在离散方法中，颗粒群被离散化为有限数量的粒径区间。这种方法的优点在于可以直接计算颗粒粒径分布（PSD）。当颗粒尺寸范围已知且不超过两到三个数量级时，这种方法尤为适用。在这种情况下，可以用相对较少的粒径区间对颗粒群进行离散化，并计算与流体动力学耦合的粒径分布。离散方法的缺点是，如果需要大量的区间，计算成本会很高。

14.9.1.2 非均匀离散方法

现有均匀离散方法的一个局限性在于，所有区间都被分配给相同的次相，因此以相同的相动量进行平流。这对于模型中大颗粒和小颗粒由于不同的动量场而可能发生分离的情况是不适用的。非均匀离散方法通过允许将区间组以不同的相速度进行平流，克服了这一局限。因此，当激活非均匀离散模型时，群体平衡模型可以应用于多个次相。

离散区间分数 $f_{i}$ 的一般输运方程可以写为：

$\frac{\partial}{\partial t} (ρ α f_{i}) + \nabla \cdot (u_{p} α f_{i}) = S_{bi} (14.689)$

由于在均匀离散方法中，所有箱子都属于同一相，因此在发生破裂和团聚的情况下，该相的净质量源为零，可以表示为：

$i = 1 \sum M S_{bi} = 0 (14.690)$

如图14.16所示（第785页），均匀离散方法中所有箱子都以相同的相速度 $u_{p}$ 进行平流。相比之下，图14.17所示（第785页）的非均匀离散方法允许将箱子分配给多个相。在这里，每个相有 $M$ 个箱子，分布在 $N$ 个相上，总共是 $M \times N$ 个箱子。箱子 $f_{1}$ 和 $f_{M}$ 分别以相速度 $u_{P_{1}}$ 等进行平流。由于给定相中的箱子可以通过破碎或团聚迁移到另一个相，因此任何给定相中箱子源的总和不一定等于零，从而为该相创造了一个净质量源。

给定相的净质量源可以表示为属于该相的箱子源的总和。

$S_{i} = i = 1 \sum M S_{bi} (14.691)$

对于断裂和聚结，所有相源的总和为零。

$i = 1 \sum N S_{i} = 0 (14.692)$

与同质离散模型也颇为相似。

$i = 1 \sum M f_{i} = 1 (14.693)$

重要提示：非均匀离散方法目前仅限于破碎和聚结过程。

图14.16：均匀离散方法

各箱子通过破碎和聚结相互作用

所有箱子以相同的相速度平流

图14.17：非均匀离散方法

14.9.1.3 矩方法标准方法

矩方法标准方法（SMM）是离散群体平衡方法的一种高效替代方案。在这种方法中，群体平衡方程被转化为一组关于分布矩的输运方程。第 $i^{th}$ 矩是通过在整个粒子空间内对粒子属性提升到第 $i^{th}$ 次幂的数密度进行积分来定义的。通常只需要求解少数几个矩方程，通常是三到六个。这相比于离散方法可能显著减少了需要求解的方程数量。除了计算优势外，当不需要整个分布而只需要某些平均和总量来表示粒子分布时，SMM方法非常有用。通常，零阶矩表示总数密度，二阶矩表示单位体积的总表面积，三阶矩表示总质量密度。

在SMM方法中，不对尺寸分布做任何假设，并且矩方程以封闭形式表述，仅涉及矩本身的函数。然而，这种精确的封闭要求构成了严重限制，因为聚结（常数聚结核除外）和破碎现象不能写成矩的函数。

14.9.1.4 矩的正交方法

矩量配方法（QMOM）在计算成本方面与SMM具有类似的优势，但用近似闭包取代了SMM所需的精确闭包。这使得QMOM能够广泛应用于各种场景，没有任何限制。

14.9.2 群体平衡模型理论

本章概述了用于Ansys Fluent中预测颗粒生长和成核的方法的理论和控制方程。

14.9.2.1 颗粒状态向量
14.9.2.2 群体平衡方程（PBE）
14.9.2.3 求解方法
14.9.2.4 群体平衡统计

14.9.2.1 颗粒状态向量

颗粒状态向量由一组“外部坐标” $(x)$ （表示颗粒的空间位置）和“内部坐标” $(ϕ)$ （可能包括颗粒大小、成分和温度）来表征。基于这些坐标，可以假设一个数密度函数 $n (x, ϕ, t)$ ，其中 $ϕ \in Ω_{ϕ}, x \in Ω_{x}$ 。因此，在无穷小体积 $d V_{x} d V_{ϕ}$ 中的颗粒平均数量为 $n (x, ϕ, t) d V_{x} d V_{ϕ}$ 。相比之下，连续相状态向量由 $Y \equiv [Y_{1} (x, t), Y_{2} (x, t), \dots, Y_{c} (x, t)]$ 给出。

整个系统中的颗粒总数为

$\int_{Ω_{ϕ}} \int_{Ω} n d V_{x} d V_{ϕ} (14.694)$

物理空间中的局域平均数密度（即单位体积内的粒子总数）由以下公式给出：

$N (x, t) = \int_{Ω_{ϕ}} n d V_{ϕ} (14.695)$

所有颗粒的总体积分数由以下公式给出：

$α (x, t) = \int_{Ω_{ϕ}} nV (ϕ) d V_{ϕ} (14.696)$

其中， $V (ϕ)$ 表示处于状态 $ϕ$ 的粒子的体积。

14.9.2.2 群体平衡方程（PBE）

假设 $ϕ$ 代表粒子体积，那么数量密度函数的传输方程可表示为

$\frac{\partial}{\partial t} [n (V, t)] + \nabla \cdot [u n (V, t)] + \nabla \cdot [G_{ν} n (V, t)] = Birth due to aggregation \frac{1}{2} \int_{0}^{V} a (V - V^{^{'}}, V^{^{'}}) n (V - V^{^{'}}, t) n (V^{^{'}}, t) d V^{^{'}} Death due to aggregation - 0 \int \infty a (V, V^{'}) n (V, t) n (V^{'}, t) d V^{'} Birth due to breakage + \int_{Ω_{ν}} p g (V^{^{'}}) β (V ∣ V^{^{'}}) n (V^{^{'}}, t) d V^{^{'}} Death due to breakage - g (V) n (V, t) (14.697)$

边界条件和初始条件设定如下：

$n (V, t = 0) = n_{v}; n (V = 0, t) G_{v} = \overset{n}{˙}_{0} (14.698)$

其中， $\overset{n}{˙}_{0}$ 表示成核速率，单位是粒子数每立方米每秒。

14.9.2.2.1 粒子生长

基于粒子体积的生长速率，记作 $G_{v}$ ，其单位为立方米每秒，定义为

$G_{v} = \frac{\partial V}{\partial t} (14.699)$

基于颗粒直径（或长度）的增长率定义为

$G = \frac{\partial L}{\partial t} (14.700)$

单个粒子的体积 ( V ) 被定义为 ( K_v L^3 )，因此 ( G_v ) 与 ( G ) 之间的关系是

$G_{v} = 3 K_{v} L^{2} G (14.701)$

单个颗粒的表面积 $A$ 定义为 $K_{a} L^{2}$ 。因此，对于立方体或球体， $K_{a} = 6 K_{v}$ 。

重要提示：

在颗粒溶解的情况下，由于颗粒体积的减小，颗粒体积增长率是负的。
在由于密度（压力）变化导致的气泡膨胀情况下，方程 14.697（第 787 页）中的增长率 $G_{v}$ 是气泡体积膨胀的速率。

考虑 $m$ 为一个体积为 $V$ 且密度为 $ρ$ 的气泡的质量，

$m = ρ V \frac{d m}{d t} = ρ \frac{d V}{d t} + V \frac{d ρ}{d t} (14.702)$

对于常数 $m$ （无质量转移的气泡膨胀）， $d m / d t = 0$ ，因此

$G_{V} = - \frac{V}{ρ} \frac{d ρ}{d t} (14.703)$

在 Mills 的论文 [443]（第1082页）中可以找到更详细的解释。

14.9.2.2.2 由于破碎和聚集导致的颗粒生成与消亡

颗粒的生成和消亡是由于破碎和聚集过程引起的。破碎过程的例子包括结晶器中的晶体断裂和泡罩塔中由于液体湍流引起的气泡破碎。同样，聚集可以由结晶器中的颗粒团聚和泡罩塔反应器中的气泡合并引起。

请注意，在 FLUENT 群体平衡模型中，所有的破碎和聚集事件都被建模为二元的，除了广义破碎核。

14.9.2.2.2.1 破碎

破碎速率表达式，或称为核 [400]（第1080页），被表示为

$Ω_{b r} = F_{B} g (V^{'}) β (V ∣ V^{'}) (14.704)$

其中

$g (V^{'}) =$ 破碎频率；即，体积为 $V^{'}$ 的颗粒在单位时间内破碎的分数 $(s^{- 1})$

$β (V ∣ V^{'}) =$ 概率密度函数（PDF），表示从体积 $V^{'}$ 破碎成体积为 $V$ 的颗粒的概率 $(1/ m^{3})$

$F_{B} =$ 破碎因子，用于校准破碎核

体积为 $V$ 的颗粒因破碎而产生的出生率由以下公式给出：

$B_{b r} = \int_{Ω_{v}} p g (V^{'}) β (V ∣ V^{'}) n (V^{'}) d V^{'} (14.705)$

在单位时间内，体积为 $V^{'}$ 的粒子中，有 $g (V^{'}) n (V^{'}) d V^{'}$ 个粒子发生破碎，产生 $p g (V^{'}) n (V^{'}) d V^{'}$ 个新粒子，其中比例为 $β (V ∣ V^{'}) d V$ 的粒子体积为 $V$ 。这里， $p$ 表示每个母粒子破碎后产生的子粒子数量（例如，二元破碎产生两个粒子）。

由于破碎导致的体积为 $V$ 的粒子的死亡率由以下公式给出：

$D_{b r} = g (V) n (V) (14.706)$

PDF $β (V ∣ V^{'})$ 也被称为颗粒破碎分布函数，或子颗粒大小分布。虽然已经提出了几种破碎分布函数的形式，但必须满足以下物理约束：归一化的破碎颗粒数量总和必须为1，碎片的质量总和必须等于原始颗粒的质量，并且形成的碎片数量必须正确表示。

数学上，这些约束可以写成如下形式：

对于归一化条件：

$\int_{0}^{V^{'}} β (V ∣ V^{'}) d V = 1 (14.707)$

质量守恒

$p \int_{0}^{V^{'}} m (V) β (V ∣ V^{'}) d V = m (V^{'}) (14.708)$

对于二元断裂， $β$ 相对于 $V / V^{'} = 0.5$ 是对称的；也就是说，

$β (V^{'} - V ∣ V^{'}) = β (V ∣ V^{'}) (14.709)$

以下是Ansys Fluent中可用于计算破碎频率的模型列表：

常数值
Luo模型
Lehr模型
Ghadiri模型
Laakkonen模型
Liao模型
用户自定义模型

Ansys Fluent提供了以下模型用于计算破碎概率密度函数（PDF）：

抛物线PDF
Laakkonen PDF
多破碎碎片通用PDF
用户自定义模型

破碎频率模型以及抛物线和通用PDF将在后续章节中详细描述。

14.9.2.2.2.2 Luo和Lehr破碎核

Luo和Lehr模型是综合核模型，涵盖了破碎频率和破碎颗粒的PDF。

其一般形式是对撞击直径为 $d$ （及体积为 $V$ ）的颗粒的涡流大小 $λ$ 进行积分。积分是在无量纲涡流大小 $ξ$ = $λ$ / $d$ 上进行的。一般形式为

$Ω_{b r} (V, V^{'}) = K \int_{ξ_{m i n}}^{1} \frac{( 1 + ξ ) ^{2}}{ξ ^{n}} exp (- b ξ^{m}) d ξ (14.710)$

参数如表14.4所示：Luo模型参数（第790页）和表14.5：Lehr模型参数（第790页）。

表14.4：Luo模型参数

$K$	$n$	$b$	m
$0.9238 ε^{1/3} d^{- 2/3} α$	11/3	$12 [f^{2/3} + (1 - f)^{2/3} - 1] σ ρ^{- 1} ε^{- 2/3} d^{- 5/3} β^{- 1}$ 其中 $f = V / V^{'}$	-1/3

详细符号说明请参见 Luo [400]（第 1080 页）。

表 14.5：Lehr 模型参数

$K$	$n$	$b$	$m$
$1.19 ε^{- 1/3} d^{- 7/3} σ ρ^{- 1} f^{- 1/3}$	13/3	$2 W e_{crit} σ ρ^{- 1} ε^{- 2/3} d^{- 5/3} f^{- 1/3}$	-2/3

其中， $W e_{crit}$ 通过图形用户界面输入。参见 Lehr [353]（第 1077 页）。

14.9.2.2.2.3 Ghadiri 破碎核

与 Luo 和 Lehr 模型相反，Ghadiri 模型 [201]（第 1068 页），[460]（第 1083 页）仅用于模拟固体颗粒的破碎频率。您需要指定概率密度函数（PDF）模型来定义子代分布。

破碎频率 $g$ 定义为：

$g (V^{'}) = K_{b} v^{2} d^{5/3} (14.711)$

其中， $v$ 表示颗粒速度， $d$ 表示破碎前的颗粒直径。破碎常数 $K_{b} [s^{3} m^{- 11 /3}]$ 决定了每秒钟有多少母颗粒发生破碎的总体规模。这一常数取决于具体情况，并与颗粒的特性成正比：

$K_{b} \sim \frac{ρ E ^{2/3}}{Γ ^{5/3}} (14.712)$

其中， $ρ$ 表示颗粒密度， $E$ 为颗粒的弹性模量， $Γ$ 则是界面能。请注意，您可能需要调整 $K_{b}$ ，例如通过与基准实验校准来进行调整。

14.9.2.2.2.4 Laakkonen 破碎核

Laakkonen 破碎核表示为破碎频率 $g (V^{'})$ 与子代概率密度函数 $β (V, V^{'})$ 的乘积，其中

$g (V^{'}) = C_{2} ε^{1/3} erfc (C_{3} \frac{σ}{ρ _{L} ε ^{2/3} d ^{5/3}} + C_{4} \frac{μ _{L}}{ρ _{L} ρ _{G} ε ^{1/3} d ^{4/3}}) (14.713)$

其中， $ε$ 表示液相涡流耗散率， $σ$ 为表面张力， $ρ_{l}$ 是液体密度， $ρ_{g}$ 是气体密度， $d$ 为母粒子直径， $μ_{L}$ 是液体粘度。

常数 $C_{2} = 2.52$ , $C_{3} = 0.04$ 和 $C_{4} = 0.01$ 。

子粒径分布函数（daughter PDF）由以下公式给出：

$β (V, V^{'}) = \frac{30}{V ^{'}} (\frac{V}{V ^{'}})^{2} (1 - \frac{V}{V ^{'}})^{2} (14.714)$

其中， $V$ 和 $V^{'}$ 分别表示子颗粒和母颗粒的体积。该模型是广泛使用的 Luo 模型的有用替代方案，因为它具有简单的子颗粒概率密度函数表达式，因此需要显著较少的计算工作量。

14.9.2.2.2.5 Liao 破碎核函数

在流动液体中的气泡会受到应力 $(τ)$ ，该应力会使气泡变形，一旦应力超过表面张力的恢复效应 $(τ_{crit})$ ，气泡最终会破碎，如图 14.18 所示：气泡破碎（第 791 页）（改编自 [369]（第 1078 页））。

图 14.18：气泡破碎

基于 Martínez-Bazán 等人的工作 [416]（第 1081 页），Liao [369]（第 1078 页）提出将尺寸为 $d_{i}$ 的气泡的破碎频率定义为应力 $τ$ 和临界应力 $τ_{crit}$ 的关系如下：

$Ω_{b r} (V, V^{'}) = {\frac{1}{d _{i}} \frac{τ ( d _{i} ) - τ _{cr i t} ( d _{i} , d _{j} )}{A} 0 for (τ > τ_{cr i t}) for (τ \leq τ_{cr i t}) (14.715)$

其中， $d_{j}$ 表示子气泡的直径。

结合破碎方法用于根据能量准则和力准则来确定临界应力。对于二元破碎，即气泡 $i$ 分裂成两个气泡 $j$ 和 $k$ 时， $τ_{crit}$ 可以表示为：

$τ_{crit} = max (τ_{crit, 1}, τ_{crit, 2}) (14.716)$

能量约束 $τ_{cr i t, 1}$ 由断裂过程中表面能的增加所决定：

$τ_{crit, 1} = E_{σ, j} + E_{σ, k} - E_{σ, i} = \frac{6 σ}{d _{i}} ((\frac{d _{j}}{d _{i}})^{2} + (\frac{d _{k}}{d _{i}})^{2} - 1) (14.717)$

其中， $E_{σ}$ 表示气泡的表面能量，而气泡 $k$ 的大小则是根据气泡 $i$ 和 $j$ 的大小计算得出的，具体如下：

$d_{k} = (d_{i}^{3} + d_{j}^{3})^{1/3} (14.718)$

力约束 $τ_{cr i t, 2}$ 由最小子气泡的毛细力决定：

$τ_{crit, 2} = \frac{σ}{min ( d _{k} , d _{j} )} (14.719)$

其中， $σ$ 是表面张力系数。

破坏应力 $τ_{crit}$ 导致气泡破裂，这是由气泡表面上的速度梯度引起的。在湍流中，速度变化可能由多种机制引起，例如湍流脉动、平均流动中的速度梯度、湍流涡旋和界面阻力 [369]（第 1078 页），如图 14.19 所示：湍流中的气泡破碎机制（第 792 页）（改编自 [369]（第 1078 页））。

图 14.19：湍流中的气泡破碎机制

这些破碎机制的应力表示为：

由湍流速度脉动引起的惯性力：

$τ_{turb} = B_{turb} ρ_{l} (ε_{l} d_{i})^{2/3} for (d_{i} > η) (14.720)$

由于主流流速梯度引起的粘性剪切力：

其中， $ρ_{l}$ 是液体密度， $ε_{l}$ 是液体相中的湍流耗散率， $η$ 是Kolmogorov微尺度。

$τ_{shear} = B_{shear} μ_{l} \overset{γ}{˙}_{b} (14.721)$

其中， $μ_{l}$ 表示液体的动力粘度， $\overset{γ}{˙}_{b}$ 则是体相中的剪切应变率。

由涡流中速度梯度引起的粘性剪切力：

$τ_{eddy} = B_{eddy} μ_{l} \overset{γ}{˙}_{eddy} for (d_{i} \leq η) (14.722)$

其中， $\overset{γ}{˙}_{e dd y}$ 表示最小涡流中流动的特征应变率。

界面阻力（或摩擦）：

$τ_{fric} = B_{fric} 0.5 ρ_{l} u_{T, i}^{2} C_{D, i} (14.723)$

其中， $u_{T, i}$ 和 $C_{D, i}$ 分别表示气泡 $i$ 的终端速度和阻力系数。

常数 $B_{turb}$ 、 $B_{shear}$ 、 $B_{eddy}$ 和 $B_{fric}$ 用于考虑不同机制引起的破坏应力估计中的不确定性。系数默认（并由 Liao [369] (第 1078 页) 推荐）值为： $B_{turb} = B_{shear} = B_{eddy} = 1$ ，以及 $B_{fric} = 0.25$ 。

14.9.2.2.2.6 抛物线型概率密度函数 (PDF)

破碎概率密度函数 (PDF) 包含了由破碎事件形成的碎片的概率信息。它提供了颗粒的数量以及从破碎中可能得到的大小分布。Ansys Fluent 中实现的抛物线型 PDF 形式使您能够定义破碎 PDF，从而...

$β (V ∣ V^{'}) = 0.5 [\frac{C}{V ^{'}} + \frac{1 - C /2}{V ^{'}} {24 (\frac{V}{V ^{'}})^{2} - 24 (\frac{V}{V ^{'}}) + 6}] (14.724)$

其中， $V$ 和 $V^{'}$ 分别表示子颗粒和母颗粒的体积。形状因子 $C$ 的取值不同，颗粒破碎分布函数的形状也会表现出不同的行为。例如，若 $C = 2$ ，颗粒破碎呈现均匀分布。若 $0 < C < 2$ ，则得到一个凹形抛物线，意味着获得大小不均的碎片比均等大小的碎片更为可能。相反，对于 $2 < C < 3$ 的情况，则更倾向于获得均等大小的碎片。超出 0 到 3 范围的值是不允许的，因为概率密度函数（PDF）不能取负值。

需要注意的是，方程 14.724（第 793 页）中定义的 PDF 关于 $V / V^{'} = 0.5$ 对称。

14.9.2.2.2.7 广义 PDF

Ansys Fluent 中实现的广义 PDF 形式使您能够模拟多个破碎碎片（即大于两个），并指定子分布的形式（例如，均匀分布、等尺寸、磨损、幂律、抛物线、二元贝塔分布）。该模型既可应用于离散方法，也可用于 QMOM。

考虑到自相似公式 [539]（第 1088 页），其中相似性 $z$ 是子颗粒与母颗粒大小的比值（即 $z \equiv \frac{V}{V ^{'}}$ ），那么广义 PDF 可表示为

$pβ (V ∣ V^{'}) = \frac{θ ( z )}{V ^{'}} (14.725)$

其中， $θ (z)$ 是自相似的女儿分布 [141]（第1065页）。 $θ (z)$ 的第 $k$ 阶矩，记作 $b_{k}$ ，

$b_{k} = \int_{0}^{1} z^{k} θ (z) d z = \frac{B _{k} ( V ^{'} )}{V ^{' k}} (14.726)$

哪里

$B_{k} (V^{'}) = \int_{0}^{V^{'}} V^{k} pβ (V ∣ V^{'}) d V (14.727)$

随后，数量守恒和质量守恒的条件可以表述为

$b_{0} = \int_{0}^{1} θ (z) d z = p (14.728)$

$b_{1} = \int_{0}^{1} z θ (z) d z = 1 (14.729)$

$θ (z)$ [141]（第1065页）的广义形式可以表示为

$θ (z) = i \sum w_{i} p_{i} \frac{z ^{q_{i} - 1} ( 1 - z ) ^{r_{i} - 1}}{β ( q _{i} , r _{i} )} (14.730)$

其中， $i$ 可以是0或1，分别表示 $θ (z)$ 由1个或2个项组成。对于每个项， $w_{i}$ 是权重因子， $p_{i}$ 是子颗粒的平均数量， $q_{i}$ 和 $r_{i}$ 是指数，而 $β (q_{i}, r_{i})$ 是beta函数。在方程14.730（第794页）中，对参数施加了以下约束：

$i \sum w_{i} = 1 (14.731)$

$i \sum w_{i} p_{i} = p (14.732)$

$i \sum w_{i} (\frac{p _{i} q _{i}}{q _{i} + r _{i}}) = 1 (14.733)$

为了展示如何将广义PDF转换为表示适当的子分布，请考虑以下表格中所示的表达式：

表14.6：子分布

类型	$θ (z)$
等大小 [321] (p. 1075)	$p δ (z - \frac{1}{p})$
磨损 [321] (p. 1075)	$δ (z - (1 - ε)) + δ (z - ε)$
幂律 [674] (p. 1096)	$(ν + 1) z^{ν - 1}$
抛物线 -a [674] (p. 1096)	$(ν + 2) (ν + 1) z^{ν - 1} (1 - z)$
奥斯汀 [34] (p. 1059)	$w (ν_{_{1}} + 1) z^{ν_{_{1}} - 1} + (1 - w) (ν_{_{2}} + 1) z^{ν_{_{2}} - 1}$
二元贝塔 -a [259] (p. 1072)	$60 z^{2} (1 - z)^{2}$
二元贝塔 -b [420] (p. 1081)	$\frac{2}{β ( ν , ν )} z^{ν - 1} (1 - z)^{ν - 1}$
均匀 [674] (p. 1096)	$p (p - 1) (1 - z)^{p - 2}$

表14.7：女儿分布（续）

类型	$p$	约束条件
等大小 [321] (p. 1075)	p	$p \geq 2$
磨损 [321] (p. 1075)	2	$ε ≪ 1$
幂律 [674] (p. 1096)	$\frac{ν + 1}{ν}$	$0 < v \leq 1$
抛物线 -a [674] (p. 1096)	$\frac{v + 2}{v}$	$0 < v \leq 2$
奥斯汀 [34] (p. 1059)	$w (1 + \frac{1}{ν _{1}}) + (1 - w) (1 + \frac{1}{ν _{2}})$	$v_{1}, v_{2} > 0; 1 \geq w \geq v_{1} (\frac{v _{2} - 1}{v _{2} - v _{1}})$
二元贝塔 -a [259] (p. 1072)	2	N/A
二元贝塔 -b [420] (p. 1081)	2	$v 0$
均匀 [674] (p. 1096)	$p$	$p \geq 2$

在表14.6：子分布（第795页）中， $δ$ 是狄拉克δ函数， $w$ 是权重系数，而 $ε, v, v_{1}$ 和 $v_{2}$ 是用户定义的参数。

通过使用表14.8：广义形式中子分布的值（第796页）所示的数值，广义形式可以表示表14.6：子分布（第795页）中的子分布。

表14.8：广义形式中子分布的值

类型	w0	$p 0$	$q 0$	$r 0$	w1	$p 1$	$q 1$	$r 1$	约束条件
等尺寸	1	$p$	0	00*	N/A	N/A	N/A	N/A	$p \geq 2$
损耗	0.5	2	$ε$	1	0.5	2	1	$ε$	$ε ≪ 1$
幂律	1	$\frac{v + 1}{v}$	$v$	1	N/A	N/A	N/A	N/A	$0 < v \leq 1$
抛物线	1	$\frac{v + 2}{v}$	$v$	2	N/A	N/A	N/A	N/A	$0 < v \leq 2$
奥斯汀	w	$\frac{v _{1} + 1}{v _{1}}$	$V_{1}$	1	$1 - w$	$\frac{v _{2} + 1}{v _{2}}$	$V_{2}$	1	$v_{1}, v_{2} > 0$ ; $1 \geq w \geq ν_{_{1}} (\frac{ν _{_{2}} - 1}{ν _{_{2}} - ν _{_{1}}})$
二元贝塔 **	1	2	$v$	$v$	N/A	N/A	N/A	N/A	$v > 0$
均匀分布	1	$p$	1	$p - 1$	N/A	N/A	N/A	N/A	$p \geq 2$

(*) 你可以通过使用一个非常大的数字来近似无穷大，例如 1e20。

(**) 当 $v = 3$ 时，Binary Beta -a 是 Binary Beta -b 的一个特例。

重要提示：请注意，对于 Ansys Fluent 中广义形式的 PDF 实现，您只需输入 $w_{0}, p_{0}, q_{0}, r_{0}$ 和 $q_{1}$ 的值，其余的值 $(w_{1}, p_{1}$ 和 $r_{1}$ ) 将自动计算。

14.9.2.2.2.8. Martinez-Bazan 破碎核

Martinez-Bazan 破碎频率 [416]（第 1081 页）通常用于模拟搅拌槽反应器中的气泡破碎。该核提供了一个基于湍流应力和表面张力局部值的气泡破碎速率的现象学表达式。破碎速率表示为：

$g (V^{'}) = {K_{g} \frac{β ( ε d _{i} ) ^{2/3} - 12 σ / ( ρ _{l} d _{i} )}{d _{i}} 0 for d_{i} > d_{m a x} for d_{i} \leq d_{m a x} (14.734)$

最大稳定气泡直径 $d_{m a x}$ 由以下公式给出：

$d_{m a x} = 2^{1/5} \frac{σ ^{3/5}}{ρ _{l}^{3/5} ε ^{2/5}} (14.735)$

在上述方程中，

$K_{g} = 0.25$ 和 $β = 8.2 =$ 实验确定的常数

$d_{i} =$ 母泡直径

$ε =$ 湍流耗散率

$σ =$ 表面张力

$ρ_{l} =$ 液体密度

14.9.2.2.2.9 倒U形概率密度函数

大多数气泡破碎的子分布概率密度函数（PDFs）假设了一个对称的U形分布，其中等尺寸子气泡的破碎概率最小。尽管这种子分布函数在许多应用中效果良好，但[416]（第1081页）中还提出了一种替代的倒U形分布，与U形分布相反，等尺寸子气泡的形成概率被假设为最大。在Ansys Fluent中，这一分布实现如下：

$β (V, V^{'}) = \frac{( X ^{2/3} - Λ ^{5/3} ) [ ( 1 - X ^{3} ) ^{2/9} - Λ ^{5/3} ]}{\int _{0}^{V^{'}} ( X ^{2/3} - Λ ^{5/3} ) [ ( 1 - X ^{3} ) ^{2/9} - Λ ^{5/3} ] d V} (14.736)$

其中， $Λ = D_{c} / D_{0}$ ， $D_{c} = [12 σ / (βρ)]^{3/5} ε^{- 2/5}$ ， $X = (\frac{V}{V ^{'}})^{1/3}$ ，而 $β (V, V^{'})$ 的定义见破裂部分（第788页）。

14.9.2.2.2.10. 聚合

聚合核函数[400]（第1080页）表示为

$Ω_{a g} = F_{A} a (V, V^{'}) (14.737)$

其中， $F_{A}$ 是一个用于校准聚合核的聚合因子。

聚合核的单位是 $m^{3} / s$ ，有时被定义为两个量的乘积：

体积为 $V$ 的粒子与体积为 $V^{'}$ 的粒子之间的碰撞频率
“聚合效率”（即体积为 $V$ 的粒子与体积为 $V^{'}$ 的粒子聚合的概率）。

体积为 $V$ 的粒子因聚合而产生的生成速率由下式给出：

$B_{a g} = \frac{1}{2} \int_{0}^{V} a (V - V^{'}, V^{'}) n (V - V^{'}) n (V^{'}) d V^{'} (14.738)$

体积为 $V - V^{'}$ 的粒子与体积为 $V^{'}$ 的粒子聚集，形成体积为 $V$ 的粒子。因子 $1/2$ 的引入是为了避免对每个碰撞事件进行重复计数。

体积为 $V$ 的粒子因聚集导致的消亡速率由以下公式给出：

$D_{a g} = \int_{0}^{\infty} a (V, V^{'}) n (V) n (V^{'}) d V^{'} (14.739)$

重要提示：破碎和聚结核函数取决于物理应用的性质。例如，在气液分散体系中，这些核函数是局部液相湍流耗散率的函数。

以下是Ansys Fluent中可用的聚结函数列表：

常数
Luo模型
自由分子模型
湍流模型
Prince和Blanch模型
Liao模型
用户自定义模型

Luo模型、自由分子模型、湍流模型、Prince和Blanch模型以及Liao模型的聚结函数将在后续章节中详细描述。

14.9.2.2.2.11.Luo聚结核函数

对于Luo模型[398]（第1080页），一般的聚结核函数定义为由于体积为 $V_{i}$ 和 $V_{j}$ 的颗粒二元碰撞而形成的颗粒体积速率：

$a (V_{i}, V_{j}) = ω_{a g} (V_{i}, V_{j}) P_{a g} (V_{i}, V_{j}) [m^{3} / sec] (14.740)$

其中， $ω_{a g} (V_{i}, V_{j}) [m^{3} / sec]$ 表示碰撞频率，而 $P_{a g} (V_{i}, V_{j})$ 则是碰撞导致聚结的概率。碰撞频率的定义如下：

$ω_{a g} (V_{i}, V_{j}) = \frac{π}{4} (d_{i} + d_{j})^{2} \overset{u}{ˉ}_{ij} (14.741)$

其中， $\overset{u}{ˉ}_{ij}$ 表示直径为 $d_{i}$ 和 $d_{j}$ 的两个粒子之间的特征碰撞速度。

$\overset{u}{ˉ}_{ij} = (\overset{u}{ˉ}_{i}^{2} + \overset{u}{ˉ}_{j}^{2})^{1/2} (14.742)$

以下文本的中文翻译：

$\overset{u}{ˉ}_{i} = 1.43 (ε d_{i})^{1/3} (14.743)$

聚合概率的表达式为：

$P_{a g} = exp {- c_{1} \frac{[ 0.75 ( 1 + x _{ij}^{2} ) ( 1 + x _{ij}^{3} ) ] ^{1/2}}{( ρ _{2} / ρ _{1} + 0.5 ) ^{1/2} ( 1 + x _{ij} ) ^{3}} W e_{ij}^{1/2}} (14.744)$

其中， $c_{1}$ 是一个数量级为1的常数， $x_{ij} = d_{i} / d_{j}$ ， $ρ_{1}$ 和 $ρ_{2}$ 分别表示主相和次相的密度，而韦伯数定义为

$W e_{ij} = \frac{ρ _{l} d _{i} ( u ˉ _{ij} ) ^{2}}{σ} (14.745)$

14.9.2.2.2.12 自由分子聚集核

真实粒子以频率（或核）聚集和破碎，这些频率（或核）具有复杂的依赖关系，涉及粒子的内部坐标 [584]（第1090页）。特别是，非常小的粒子（例如，直径可达 $1 μ m$ ）由于布朗运动引起的碰撞而聚集。在这种情况下，碰撞频率与粒子大小有关，通常采用以下核函数：

$a (V_{i}, V_{j}) = \frac{2 k _{B} T}{3 μ} \frac{( L _{i} + L _{j} ) ^{2}}{L _{i} L _{j}} (14.746)$

其中， $k_{B}$ 是玻尔兹曼常数， $T$ 是绝对温度， $μ$ 是悬浮流体的粘度， $L_{i}$ 和 $L_{j}$ 分别是粒子 $i$ 和 $j$ 的直径（或长度）。这个核函数也被称为布朗核或邻位核。

14.9.2.2.2.13 湍流聚合核

在混合过程中，机械能被提供给流体。这种能量在流体中产生湍流。湍流产生涡旋，进而有助于能量耗散。能量从最大的涡旋传递到最小的涡旋，在其中通过粘性相互作用耗散。最小的涡旋的尺寸是科尔莫戈罗夫微尺度， $η$ ，它表示为运动粘度和湍流动能耗散率的函数：

$η = (\frac{V ^{3}}{ε})^{1/4} (14.736)$

在湍流流场中，聚集可以通过两种机制发生：

粘性子范围机制：当颗粒小于科尔莫戈罗夫微尺度 $η$ 时，适用此机制。
惯性子范围机制：当颗粒大于科尔莫戈罗夫微尺度时，适用此机制。在这种情况下，颗粒具有独立的速率。

对于粘性子范围，颗粒碰撞受涡流内局部剪切的影响。根据萨夫曼和特纳 [564]（第1089页）的工作，碰撞率表示为，

$a (V_{i}, V_{j}) = ς_{T} \frac{8 π}{15} \overset{γ}{˙} \frac{( L _{i} + L _{j} ) ^{3}}{8} (14.748)$

其中， $ς_{T}$ 是一个预因子，考虑了湍流碰撞的捕获效率系数，而 $\overset{γ}{˙}$ 是剪切速率：

$\overset{γ}{˙} = \frac{ε}{v}^{0.5} (14.749)$

在惯性子区，颗粒大于最小的涡流，因此它们受到流场中速度波动的影响而被拖拽。在这种情况下，聚并速率采用Abrahamson模型[6]（第1057页）来表达。

$a (V_{i}, V_{j}) = ς_{T} 2^{3/2} π \frac{( L _{i} + L _{j} ) ^{2}}{4} (U_{i}^{2} + U_{j}^{2}) (14.750)$

其中， $U_{i}^{2}$ 表示粒子 $i$ 的均方速度。

湍流碰撞的经验捕获效率系数描述了碰撞粒子之间的水动力和吸引相互作用。Higashitani等人[245]（第1071页）提出了以下关系：

$ς_{T} = 0.732 (\frac{5}{N _{T}})^{0.242}; N_{T} \geq 5 (14.751)$

其中， $N_{T}$ 表示黏性力与范德华力之间的比值。

$N_{T} = \frac{6 π μ ( L _{i} + L _{j} ) ^{3} λ ˙}{8 H} (14.752)$

其中， $H$ 为哈梅克常数，是颗粒材料的函数，而 $\dot{λ}$ 则为变形速率。

$\dot{λ} = (\frac{4 ε}{15 π v})^{0.5} (14.753)$

14.9.2.2.2.14 Prince和Blanch聚合核

Prince和Blanch的聚合模型[535]（第1088页）假设两个气泡的聚合过程分为三个步骤：

气泡碰撞，其间夹带少量液体。
分离气泡的液膜排液，直至达到临界厚度。
液膜破裂，气泡合并成一个大气泡。

此聚合核通过两个气泡的碰撞速率以及与聚合所需时间相关的碰撞效率来建模。

$a (V_{i}, V_{j}) = (θ_{ij}^{T} + θ_{ij}^{B} + θ_{ij}^{S}) η_{ij} (14.754)$

碰撞效率通过比较聚并所需时间 ( t_{ij} ) 与碰撞过程中实际接触时间 ( \tau_{ij} ) 来进行建模：

$η_{ij} = e^{- t_{ij} / τ_{ij}} (14.755)$

其中：

$t_{ij} = (\frac{ρ _{l} r _{ij}^{3}}{16 σ})^{1/2} ln (\frac{h _{o}}{h _{f}}) (14.756)$

及：

$τ_{ij} = \frac{r _{ij}^{2/3}}{ε _{l}^{1/3}} (14.757)$

在公式14.756（第800页）和公式14.757（第801页）中， $ρ_{1}$ 表示液体密度， $σ$ 代表表面张力， $h_{o}$ 是初始膜厚度， $h_{f}$ 是破裂发生时的临界膜厚度， $ε_{l}$ 表示液体的湍流涡流耗散，而 $r_{ij}$ 则是定义为等效半径的参数：

$r_{ij} = (\frac{1}{2} (\frac{1}{r _{i}} + \frac{1}{r _{j}}))^{- 1} (14.758)$

其中， $r_{i}$ 和 $r_{j}$ 分别表示气泡 $i$ 和 $j$ 的半径。

湍流对碰撞频率的贡献被建模为：

$θ_{ij}^{T} = F_{CT} S_{ij} (u_{t i}^{2} + u_{t j}^{2})^{1/2} (14.759)$

碰撞粒子横截面积的定义如下：

$S_{ij} = \frac{π}{4} (d_{i} + d_{j})^{2} (14.760)$

湍流速度由以下公式给出：

$u_{t i} = 2 ε_{l}^{1/3} d_{i}^{1/3} (14.761)$

$F_{CT}$ 是一个校准因子，用于调整湍流贡献。

浮力对碰撞频率的贡献计算方式如下：

$θ_{ij}^{B} = F_{CB} S_{ij} ∣ U_{r j} - U_{r i} ∣ (14.762)$

其中：

$U_{r i} = \frac{2.14 σ}{ρ _{l}} + 0.505 g d_{i} (14.763)$

而 $F_{CB}$ 是浮力的校准因子。

在公式 14.760（第 801 页）、公式 14.761（第 801 页）和公式 14.763（第 801 页）中， $d_{i}$ 和 $d_{j}$ 分别是气泡 $i$ 和 $j$ 的直径， $g$ 是重力加速度。

目前忽略了剪切对碰撞频率的贡献（公式 14.754（第 800 页）中的 $θ_{ij}^{S}$ ）。

14.9.2.2.2.15 廖氏聚合核

聚结发生在两个或多个流体颗粒碰撞时。这可能由于多种原因发生，如湍流、涡流捕获、平均流动中的速度梯度、体积力以及尾流卷吸等，如图 14.20 所示：湍流中的气泡碰撞（第 802 页）（改编自 [369]（第 1078 页））。

图 14.20：湍流中的气泡碰撞

假设所有聚结机制的累积贡献，廖 [369]（第 1078 页）提出了一种计算两个体积分别为 $V_{i}$ 和 $V_{j}$ 的气泡之间二元碰撞总聚结频率 $a$ 的公式。

$a (V_{i}, V_{j}) = \frac{α _{g, max}}{α _{g, max} + α _{g}} (Γ_{turb} + Γ_{eddy} + Γ_{shear} + Γ_{buoy} + Γ_{wake}) (14.764)$

其中， $α_{g}$ 表示气泡填充极限，而 $α_{g, m a x}$ 是其最大值。当气泡相的局部体积分数接近最大填充极限 $α_{g, m a x}$ 时，气泡合并速率无限增大。

由于气泡间的碰撞可能导致合并或反弹，因此使用合并效率（或概率）。合并频率 $a$ 表示为碰撞频率 $h$ 与每次碰撞事件效率 $λ$ 的乘积：

$a = hλ (14.765)$

类似于动力气体理论，碰撞频率 ( h ) 近似为两个粒子在单位时间内扫过的体积：

$h = A u_{re l} (14.766)$

其中， $u_{re l}$ 表示两气泡间的相对速度，而 $A$ 则为碰撞的有效横截面积。

对于由湍流波动引起的随机气泡碰撞，其有效横截面积为：

$A_{turb} = \frac{π}{4} (d_{i} + d_{j})^{2} (14.767)$

其中， $d_{i}$ 和 $d_{j}$ 分别代表气泡 $i$ 和 $j$ 的直径。

气泡 $i$ 和 $j$ 之间的相对速度估计如下：

$u_{re l, t u r b} = (u_{t u r b, i}^{2} + u_{t u r b, j}^{2})^{1/2} (14.768)$

其中， $u_{turb}$ 是在距离 $d_{i}$ 上的湍流速度，具体由以下公式给出：

$u_{turb} = 2 (ε_{l} d_{i})^{1/3} (14.769)$

其中， $ε_{l}$ 表示液相中的湍流动能耗散率。

对于由湍流引起的碰撞，气泡 $i$ 与气泡 $j$ 之间的聚并频率 $Γ_{turb}$ 可表示为：

$Γ_{turb} = \frac{π}{4} (d_{i} + d_{j})^{2} C_{turb} 2 (d_{i}^{2/3} + d_{j}^{2/3})^{1/2} ε_{l}^{1/3} λ_{inertial} (14.770)$

其中， $C_{turb}$ 是一个常数，而 $λ_{inertial}$ 是在惯性碰撞区域中的聚并效率。方程 14.770（第 803 页）适用于远大于 Kolmogorov 长度尺度的气泡。

当考虑浮力效应时，只有在向上流动中，较快的气泡从下方接近较慢的气泡时，才会发生碰撞。聚并频率 $Γ_{buoy}$ 被建模为：

$Γ_{buoy} = 0.5 \frac{π}{4} (d_{i} + d_{j})^{2} C_{buoy} ∣ u_{T, i} - u_{T, j} ∣ λ_{eff} (14.771)$

其中， $C_{buoy}$ 是一个常数，而 $u_{T, i}$ 和 $u_{T, j}$ 分别表示气泡 $i$ 和气泡 $j$ 的终端速度。有效聚并效率 $λ_{e ff}$ 的计算公式如下：

$λ_{e ff} = {λ_{inertial} λ_{viscous} for (d_{i} + d_{j} > η) for (d_{i} + d_{j} \leq η) (14.772)$

其中， $η$ 是 Kolmogorov 长度尺度。

惯性碰撞的聚并效率 $λ_{inertial}$ 表示为：

$λ_{inertial} = exp (C_{eff} W e_{ij, max}^{0.5}) (14.773)$

系数 $C_{e ff}$ 的默认值为 5.0，

$W e_{ij, max} = \frac{ρ _{l} d _{e q}}{σ} max (u_{re l, turb}, u_{re l, buoy}, u_{re l, shear})^{2} (14.774)$

其中， $ρ_{l}$ 表示液体密度， $σ$ 代表表面张力系数，而

$d_{e q} = \frac{2 d _{i} d _{j}}{d _{i} + d _{j}} (14.775)$

粘性碰撞的聚并效率 $λ_{viscous}$ 表示为：

$λ_{viscous} = exp - \frac{3 μ _{l} d _{e q} Y ˙ _{e dd y}}{4 σ} ln (\frac{πσ d _{e q}^{2}}{32 A _{H}})^{1/3} (14.776)$

其中， $μ_{l}$ 是液体的动力粘度， $A_{H}$ 是哈梅克常数，表示宏观物体间范德华相互作用的强度， $\overset{γ}{˙}_{e dd y}$ 是流体在最小涡流中的特征应变速率。对于空气-水界面，其值近似为 $3.7 \times 10^{- 20} J$ 。

在剪切诱导的碰撞中，气泡只能从高速度侧沿着主流速度方向与它跟随的另一个气泡发生碰撞。由速度梯度引起的气泡 $i$ 和气泡 $j$ 之间的合并频率 $Γ_{s h e a r}$ 表示为：

$Γ_{shear} = 0.5 \frac{π}{4} (d_{i} + d_{j})^{2} C_{shear} [\frac{0.5}{π} (d_{i} + d_{j}) \overset{γ}{˙}_{b}] λ_{eff} (14.777)$

其中， $C_{shear}$ 是一个常数， $\overset{γ}{˙}_{b}$ 是整体流动的剪切应变率，而 $λ_{e ff}$ 由方程 14.772（第 803 页）定义。

对于远小于 Kolmogorov 微尺度的气泡，其碰撞被认为是受粘性力控制的。聚并频率 $Γ_{e dd y}$ 被建模为：

$Γ_{eddy} = 0.5 \frac{π}{4} (d_{i} + d_{j})^{2} C_{eddy} [\frac{0.5}{π} (d_{i} + d_{j}) \overset{γ}{˙}_{eddy}] λ_{viscous} (14.778)$

其中， $C_{eddy}$ 是一个常数。它可以从 Kolmogorov 尺度推导得出：

$\overset{γ}{˙}_{e dd y} = \frac{ρ _{l}}{μ _{l}} (14.779)$

当一个直径大于临界尺寸的气泡在液体中移动时，它会在其后方的液体相中产生尾流。位于前导气泡尾流区域内的跟随气泡可以被加速并达到前导气泡。尾流裹挟的聚并效率 $Γ_{wake}$ 被建模为：

$Γ_{w ak e} = C_{w ak e} 0.5 \frac{π}{4} [d_{i}^{2} u_{T, i} C_{D, i}^{1/3} Θ_{i} + d_{j}^{2} u_{T, j} C_{D, j}^{1/3} Θ_{j}] (14.780)$

其中， $C_{wake}$ 是一个常数， $u_{T, i}$ 和 $u_{T, j}$ 分别是气泡 $i$ 和 $j$ 的终端速度，而 $C_{D, i}$ 和 $C_{D, j}$ 分别是气泡 $i$ 和 $j$ 的阻力系数。气泡 $i$ 的混合函数 $Θ_{i}$ 定义为：

$Θ_{i} = {\frac{( d _{i} - 0.5 d _{cr i t} ) ^{6}}{( d _{i} - 0.5 d _{cr i t} ) ^{6} + ( 0.5 d _{cr i t} ) ^{6}} 0 for d_{i} \geq 0.5 d_{cr i t} for d_{i} < 0.5 d_{cr i t} (14.781)$

临界直径 $d_{crit}$ 的计算公式如下：

$d_{cr i t} = 4 \frac{σ}{g ( ρ _{l} - ρ _{g} )} (14.782)$

其中 $ρ_{g}$ 是气体密度。

常数 $C_{turb}$ 、 $C_{buoy}$ 、 $C_{shear}$ 、 $C_{eddy}$ 和 $C_{wake}$ 用于考虑在不同条件下相对速度估计的不确定性。Liao [369]（第1078页）建议将这些系数设置为1.0，这也是Ansys Fluent中的默认设置。

14.9.2.2.3 粒子成核生成

根据应用的不同，粒子可能会由于主要相中分子的转移而自发成核。例如，在溶液结晶过程中，第一步是相分离或“新生”新晶体。在沸腾应用中，第一个蒸汽泡的产生是一个成核过程，称为核沸腾。

成核速率通过边界条件定义，如公式14.698（第787页）所示。

14.9.2.3 求解方法

如前所述，Ansys Fluent 中的人口平衡方程可以通过四种不同的方法求解：离散方法、非均匀离散方法、标准矩方法（SMM）和四阶矩方法（QMOM）。对于每种方法，Ansys Fluent 的实现仅限于对应于粒子大小的单一内部坐标。以下小节描述了每种方法的理论背景，并列出了它们的优缺点。

14.9.2.3.1 离散方法和非均匀离散方法

14.9.2.3.2 标准矩方法（SMM）

14.9.2.3.3 四阶矩方法（QMOM）

14.9.2.3.4 直接四阶矩方法（DQMOM）

14.9.2.3.1 离散方法和非均匀离散方法

离散方法（也称为分类或分段方法）由Hounslow [256]（第1071页）、Litster [384]（第1079页）和Ramkrishna [539]（第1088页）开发。该方法基于将连续的颗粒尺寸分布（PSD）表示为一组离散的尺寸类别或区间，如图14.21所示：离散方法表示的颗粒尺寸分布（第805页）。这种方法的优点在于其稳健的数值计算和直接给出PSD的能力。缺点是需要预先定义区间，并且可能需要大量的类别。

非均匀离散方法的求解方法基于离散方法，因此共享许多相同的基本原理。

图14.21：离散方法表示的颗粒尺寸分布

14.9.2.3.1.1 数值方法

在Ansys Fluent中，PBE以颗粒尺寸 $i$ 的体积分数形式编写：

$\frac{\partial}{\partial t} (ρ_{s} α_{i}) + \nabla \cdot (ρ_{s} u_{i} α_{i}) + \frac{\partial}{\partial V} (\frac{G _{v} ρ _{s} α _{i}}{V}) = ρ_{s} V_{i} (B_{a g, i} - D_{a g, i} + B_{b r, i} - D_{b r, i}) + 0^{i} ρ_{s} V_{o} \overset{n}{˙}_{0} (14.783)$

其中， $ρ_{s}$ 表示次相的密度， $α_{i}$ 表示颗粒尺寸 $i$ 的体积分数，定义为

$α_{i} = N_{i} V_{i} i = 0, 1, \dots, N - 1 (14.784)$

以下是翻译成中文的文本：

在哪里

$N_{i} (t) = \int_{V_{i}}^{V_{i + 1}} n (V, t) d V (14.785)$

而 $V_{i}$ 表示粒径为 $i$ 的颗粒体积。在 Ansys Fluent 中，引入了一个称为 $f_{i}$ 的分数作为解变量，这个分数被定义为

$f_{i} = \frac{α _{i}}{α} (14.786)$

其中， $α$ 表示次级相的总体积分数。

成核速率 $\overset{n_{0}}{˙}$ 出现在描述最小尺寸 $V_{0}$ 的体积分数的离散化方程中。符号 $0^{i}$ 表示这个特定项，即 $ρ_{s} V_{0} \overset{n}{˙}_{0}$ ，仅在最小颗粒尺寸的情况下出现在方程 14.783（第 805 页）中。

生长速率，如方程 14.783（第 805 页）[256]（第 1071 页）所示，其离散化如下：

$\frac{\partial}{\partial V} (\frac{G _{v} ρ _{s}^{’} α _{i}}{V}) = ρ_{s}^{’} V_{i} [(\frac{G _{v, i - 1}^{’} N _{i - 1}}{V _{i} - V _{i - 1}}) - (\frac{G _{v, i}^{’} N _{i}}{V _{i + 1} - V _{i}})] (14.787)$

体积坐标被离散化为 [256]（第 1071 页） $V_{i + 1} / V_{i} = 2^{q}$ ，其中 $q = 1, 2, \dots$ ，并被称为“比率因子”。

粒子的生成和消亡率定义如下：

$B_{a g, i} = k = 1 \sum N j = 1 \sum N a_{kj} N_{k} N_{j} x_{kj} ξ_{kj} (14.788)$

$D_{a g, i} = j = 1 \sum N a_{ij} N_{i} N_{j} (14.789)$

$B_{b r, i} = j = i + 1 \sum N g (V_{j}) N_{j} β (V_{i} ∣ V_{j}) (14.790)$

$D_{b r, i} = g (V_{i}) N_{i} (14.791)$

其中， $a_{ij} = a (V_{i}, V_{j})$ ，

$ξ_{kj} = {10 for V_{i} < V_{a g} < V_{i + 1}, where i \leq N - 1 otherwise (14.792)$

$V_{a g}$ 表示由粒子 $k$ 和 $j$ 聚集形成的粒子体积，其定义如下：

$V_{a g} = [x_{kj} V_{i} + (1 - x_{kj}) V_{i + 1}] (14.793)$

以下文本的中文翻译：

$x_{kj} = \frac{V _{a g} - V _{i + 1}}{V _{i} - V _{i + 1}} (14.794)$

如果 $V_{a g}$ 大于或等于最大颗粒尺寸 $V_{N}$ ，那么对类别 $N - 1$ 的贡献为

$x_{kj} = \frac{V _{a g}}{V _{N}} (14.795)$

重要提示：请注意，对于最小的颗粒类别，不存在破碎现象。

14.9.2.3.1.2 离散方法的破碎公式

Ansys Fluent 中离散方法的默认破碎公式基于 Hagesather 方法 [224]（第 1070 页）。在该方法中，破碎源被分配到相应的尺寸区间，同时保持质量和数量密度不变。当相邻区间尺寸比可以表示为 $2^{n}$ ，其中 $n = 1, 2, \dots$ 时，第 $i$ 个区间（ $i = 1, \dots, N$ ）的源项可以表示为

$B_{b} (i) = N \sum k = i + 1, i \neq = N Ω_{b} (ν_{k}, ν_{i}) + i \sum k = i, i \neq = N x_{i + 1, k} Ω_{b} (ν_{i + 1}, ν_{k}) + i - 1 \sum k = 1, i \neq = 1 (1 - x_{i, k}) Ω_{b} (v_{i}, v_{k}) (14.796)$

这里

$Ω_{b} (v_{k}, v_{i}) = N_{k} g (v_{k}) β (v_{k}, v_{i}) (14.797)$

Ramkrishna [326]（第1076页）提供了一种更为严谨的数学表达方式，其中破碎率被表示为...

$B_{b} (i) = i \sum N n_{i, k} g (v_{k}) N_{k} (14.798)$

哪里

$n_{i, k} = \int_{v_{i}} v_{i + 1} \frac{v _{i + 1} - v}{v _{i + 1} - v _{i}} β (v_{k}, v) d v + \int_{v_{i - 1}}^{v_{i}} \frac{v - v _{i - 1}}{v _{i} - v _{i - 1}} β (v_{k}, v) d v (14.799)$

Ramkrishna公式可能由于所需积分点数量庞大而计算较慢。然而，对于简单的 $β$ 形式，积分可以相对容易地进行。Hagesather公式所需的积分点较少，通过适当选择分箱大小，其与Ramkrishna公式的精度差异可以得到校正。

重要提示：为了保持计算时间合理，当Luo模型与Ramkrishna公式结合使用时，采用体积平均值来表示湍流涡流耗散。

注意：非均匀离散相应用Hagesather公式。

14.9.2.3.2 标准矩方法（SMM）

SMM由Randolph和Larson提出[540]（第1088页），是求解PBE的另一种方法。其优点在于降低了问题的维度，并且求解低阶矩的传输方程相对简单。缺点是只有在恒定聚集和尺寸无关增长的情况下才能实现右侧的精确闭合，并且无法进行破碎建模。然而，通过QMOM（参见第809页的矩量法（QMOM））克服了闭合约束。

14.9.2.3.2.1 数值方法

SMM方法基于对PBE关于内部坐标（在此情况下为颗粒尺寸 $L$ ）取矩。

定义 $k^{th}$ 矩为

$m_{k} (x, t) = \int_{0}^{\infty} n (L; x, t) L^{k} d L k = 0, 1, \dots, N - 1 (14.800)$

假设颗粒持续增长，其传输方程可以表示为：

$\frac{\partial}{\partial t} (ρ_{s} m_{k}) + \nabla \cdot (ρ_{s} u m_{k}) = ρ_{s} (\overset{ˉ}{B}_{ag, k} - \overset{ˉ}{D}_{ag, k} + \overset{ˉ}{B}_{br, k} - \overset{ˉ}{D}_{br, k}) + ρ_{s} (0^{k} \overset{n}{˙}_{0} + Growth) (14.801)$

此处

$\overset{ˉ}{B}_{a g, k} = \frac{1}{2} \int_{0}^{\infty} n (λ) \int_{0}^{\infty} a (u, λ) (u, λ) (u^{3} + λ^{3})^{k /3} n (u) d u d λ (14.802)$

$\overset{ˉ}{D}_{a g, k} = \int_{0}^{\infty} L^{k} n (L) \int_{0}^{\infty} a (L, λ) n (λ) d λ d L (14.803)$

$\overset{ˉ}{B}_{b r, k} = \int_{0}^{\infty} L^{k} \int_{0}^{\infty} g (λ) β (L ∣ λ) n (λ) d λ d L (14.804)$

$\overset{ˉ}{D}_{b r, k} = \int_{0}^{k} L^{k} g (L) n (L) d (14.805)$

$N$ 表示指定的时刻数量，而 $\overset{n_{0}}{˙}$ 则是成核速率。生长项被定义为

$Growth \equiv \int_{0}^{\infty} k L^{k - 1} G (L) n (L, t) d L (14.806)$

以下文本翻译成中文为：而持续增长则表现为

$k G m_{k - 1} (14.807)$

公式14.802可以通过以下方法推导出：

$u^{3} = L^{3} - λ^{3}; d L = \frac{u ^{2}}{L ^{2}} d u$

并反转积分顺序。基于这些矩量，可以推导出描述粒子群体总体特性的参数。

$N_{total} = m_{0} (14.808)$

$L_{total} = m_{1} (14.809)$

$A_{total} = K_{a} m_{2} (14.810)$

$V_{total} = K_{v} m_{3} (14.811)$

$d_{32} = \frac{m _{3}}{m _{2}} (14.812)$

这些特性与混合悬浮液每单位体积中固态颗粒的总数、长度、面积和体积有关。通常使用索特平均直径 $d_{32}$ 作为平均粒径。

为了封闭方程 14.801（第 808 页），需要用正在求解的矩来表示方程 14.802（第 808 页）至方程 14.805（第 808 页）中的量。为此，一种方法是假设与粒径无关的破碎和聚集核，以及其他简化，如对项 $(u^{3} + λ^{3})^{k /3}$ 进行泰勒级数展开。或者，可以假设粒度分布（PSD）的轮廓，以便可以对方程 14.802（第 808 页）至方程 14.805（第 808 页）进行积分，并用正在求解的矩来表示。

在 Ansys Fluent 中，通过限制 SMM 应用于粒径无关生长和恒定聚集核的情况，实现了精确封闭。

14.9.2.3.3 矩量法（QMOM）

矩量法（QMOM）由 McGraw [422]（第 1081 页）首次提出，用于模拟气溶胶演化和凝聚问题。Marchisio 等人 [413]（第 1081 页）的应用表明，该方法需要相对较少的标量方程来跟踪人口矩，且误差较小。

当需要聚集量而不是精确的 PSD 时，QMOM 提供了一种有吸引力的离散方法替代方案。其优点是变量较少（通常只有六个或八个矩），并且动态计算尺寸区间。缺点是可能没有足够的横坐标来描述 PSD，并且求解乘积差分算法可能耗时较长。

14.9.2.3.3.1 数值方法

正交多项式序列的确定基于与 $n (L)$ （即颗粒尺寸分布）正交的原理。若正交近似的横坐标为阶数为 $N$ 的多项式的节点，则该正交近似即为所求。

$\int_{0}^{\infty} f (L) n (L) d L \approx i = 1 \sum N f (L_{i}) w_{i} (14.813)$

如果 $f (L)$ 是一个次数为 $N$ 或更小的多项式 [140]（第1065页），则该方法是精确的。在其他情况下， $f (L)$ 越接近一个多项式，近似就越准确。

计算正交近似的一种直接方法是利用其通过矩的定义：

$m_{k} = i = 1 \sum N w_{i} L_{i}^{k} (14.814)$

阶数为 $N$ 的正交逼近由其 $N$ 个权重 $w_{i}$ 和 $N$ 个横坐标 $L_{i}$ 定义，并且可以通过其前 $2 N$ 个矩 $m_{0}, \dots, m_{2 N - 1}$ 来计算，方法是根据矩 $m_{k}$ 写出多项式的递归关系。一旦这种关系以矩阵形式写出，就很容易证明多项式的根对应于雅可比矩阵 [534]（第 1088 页）的特征值。此过程称为矩逆算法，它使用乘积差算法 [213]（第 1069 页）或惠勒算法 [705]（第 1098 页）。一旦知道了权重和横坐标，就可以计算由于聚结和断裂产生的源项，从而可以求解矩的输运方程。

应用公式 14.813（第 809 页）和公式 14.814（第 810 页），公式 14.801（第 808 页）中的产生和消亡项可以重写为

$\overset{ˉ}{B}_{a g, k} = \frac{1}{2} i = 1 \sum N w_{i} j = 1 \sum N w_{j} (L_{i}^{3} + L_{j}^{3})^{k /3} a (L_{i}, L_{j}) (14.815)$

$\overset{ˉ}{D}_{a g, k} = i = 1 \sum N L_{i}^{k} w_{i} j = 1 \sum N w_{j} a (L_{i}, L_{j}) (14.816)$

$\overset{ˉ}{B}_{b r, k} = i = 1 \sum N w_{i} \int_{0}^{\infty} L_{k} g (L_{i}) β (L ∣ L_{i}) d L (14.817)$

$\overset{ˉ}{D}_{b r, k} = i = 1 \sum N w_{i} L_{i}^{k} g (L_{i}) (14.818)$

从理论上讲，在使用QMOM时，对破碎和聚结核函数的表达没有任何限制。

成核速率的定义方式与SMM相同。QMOM的生长速率由方程14.806（第808页）定义并表示为

$i = 1 \sum N w_{i} L_{i}^{k - 1} G (L_{i}) (14.819)$

14.9.2.3.4 直接矩量法（DQMOM）

DQMOM方程通过矩传递方法从基本数密度函数方程推导而来，类似于QMOM。不同之处在于，DQMOM假设每个积分点占据一个独立的流场，而QMOM假设所有积分点在同一流场中移动。这一差异使得DQMOM能够预测由于颗粒相互作用引起的颗粒分选。

在此DQMOM实现中，必须指定四个相：一个主相和三个次相，即DQMOM相。与QMOM相比，对于一个具有三个积分点的系统，DQMOM方法仅需要三个额外的方程来求解颗粒的有效长度，但每个DQMOM相的体积分数方程有额外的源项。在Ansys Fluent中，考虑了三种颗粒相互作用：生长、聚集和破碎。未考虑成核。

14.9.2.3.4.1 数值方法

描述一个经历聚集、破碎和生长的多分散颗粒系统的DQMOM方程可以写成如下形式（DQMOM公式的详细信息可以在[169]（第1066页）中找到）：

$\frac{\partial ε _{i} ρ _{s}}{\partial t} + \nabla \cdot (u_{s i} ε_{i} ρ_{s}) = 3 k_{v} ρ_{s} L_{i}^{2} (b_{i} + w_{i} G_{i}) - 2 k_{v} ρ_{s} L_{i}^{3} a_{i} (14.820)$

$\frac{\partial ε _{i} L _{i} ρ _{s}}{\partial t} + \nabla \cdot (u_{s i} ε_{i} L_{i} ρ_{s}) = 4 k_{v} ρ_{s} L_{i}^{3} (b_{i} + w_{i} G_{i}) - 3 k_{v} ρ_{s} L_{i}^{4} a_{i} (14.821)$

其中， $ε_{i}$ 和 $ε_{i} L_{i}$ 分别代表体积分数和颗粒相的有效长度。 $w_{i}$ 表示单位体积内的颗粒数量，而 $G_{i}$ 是在第 $i$ 个积分点处的生长速率。至于 $a_{i}$ 和 $b_{i}$ ，则可以通过使用 $N$ 个积分点对颗粒数密度输运方程进行矩变换后得到的线性系统来计算。该线性系统可以用矩阵形式表示为

$A α = d (14.822)$

定义的 $2 N \times 2 N$ 系数矩阵 $A = [A_{1} A_{2}]$ 如下：

$A_{1} = 10 - L_{1}^{2} ⋮ 2 (1 - N) L_{1}^{2 N - 1} \dots \dots \dots ⋱ \dots 10 - L_{N}^{2} ⋮ 2 (1 - N) L_{N}^{2 N - 1} (14.823)$

$A_{2} = 01 2 L_{1} ⋮ (2 N - 1) L_{1}^{2 N - 2} \dots \dots \dots ⋱ \dots 01 L_{N} ⋮ (2 N - 1) L_{N}^{2 N - 2} (14.824)$

未知向量 $α$ 是一个由 $2 N$ 个元素组成的向量，其定义如下：

$α = [a_{1} \dots a_{N} b_{1} \dots b_{N}]^{T} = [a b] (14.825)$

方程14.822（第811页）右侧为仅涉及聚并和破碎现象的已知源项。生长项直接在方程14.820（第811页）和方程14.821（第811页）中考虑。目前，成核现象尚未纳入考虑。

$d = [S_{0}^{(N)} \dots S_{2 N - 1}^{(N)}]^{T} (14.826)$

第 $k^{t h}$ 阶矩的源项 $S_{k}^{(N)} (k = 0, \dots, 2 N - 1)$ 定义为

$S_{k}^{(N)} (x, t) = \int_{0}^{\infty} L^{k} S (x, t) d L (14.827)$

当正交点 $L_{i}$ 的横坐标不同时，矩阵 $A$ 是良定义的，并且可以求得方程 14.822（第 811 页）的唯一解。否则，矩阵 $A$ 不是满秩的，无法求逆以找到 $α$ 的唯一解。Ansys Fluent 采用的解决方法是使用扰动技术。例如，对于当前的三个正交点系统，扰动技术会向横坐标添加一个小的值，以确保矩阵 $A$ 是满秩的。需要注意的是，扰动技术仅用于矩阵 $A$ 的定义，并未对公式 14.827（第 811 页）的源项向量进行任何修改。因此，由聚并和破碎产生的权重和总体源项均不受扰动方法的影响。模拟测试发现，扰动方法可以稳定公式 14.822（第 811 页）的解，并减少那些横坐标值过于接近的两相的物理上不切实际的巨大源项。然而，对于横坐标 $L_{i}$ 与其他两相明显不同的相，该技术影响甚微。