颗粒运动方程

颗粒的动量守恒方程在拉格朗日架构中编写。动量的变化由作用于颗粒的表面力和体积力来平衡。

质量为 $m_{p}$ 的材料或 DEM 颗粒的（线性）动量守恒方程由以下公式给出：

m_{p} \frac{d v_{p}}{d t} = F_{s} + F_{b}

(2986)

其中， $v_{p}$ 表示瞬时颗粒速度， $F_{s}$ 是作用于颗粒表面的力的合力， $F_{b}$ 是体积力的合力。这些力进而分解为：

F_{s} = F_{d} + F_{p} + F_{v m}

(2987)

且

F_{b} = F_{g} + F_{M R F} + F_{u} + F_{c} + F_{Co}

(2988)

其中：

$F_{d}$ 为曳力。
$F_{p}$ 为压力梯度力。
$F_{v m}$ 为虚拟质量力
$F_{g}$ 为重力。
$F_{u}$ 为用户自定义的体积力。
$F_{c}$ 为接触力（仅适用于 DEM 建模）。
$F_{Co}$ 为库仑力。

表面力 $F_{s}$ 的合力表示从连续相到颗粒的动量传递。当使用双向耦合建模方法（请参见与连续相的相互作用）时， $F_{s}$ 在所有粒子束间累积并应用于连续相动量方程。

在 DEM 建模中，接触力 $F_{c}$ 表示颗粒间以及颗粒与边界的相互作用。此力在颗粒具有的接触间累积并作为体积力应用：

F_{c} = \sum_{c o n t a c t s} F_{c m}

(2989)

其中， $F_{c m}$ 是接触模型力。Simcenter STAR-CCM+ 提供了多个接触力模型来描述接触力 $F_{c m}$ 。有关详细信息，请参见接触力。

此外，颗粒具有方向，因此其角动量也必须守恒。

对于 DPM 颗粒，设定为：

I_{p} \frac{{d ω}_{p}}{d t} = M_{b} + M_{u}

(2990)

其中：

$I_{p}$ 为颗粒惯性矩，由二阶张量描述。
$ω_{p}$ 为颗粒角速度。
$M_{b}$ 为曳力扭矩，即因旋转曳力而作用于颗粒的力矩（曳力扭矩）。
$M_{u}$ 为用户自定义扭矩。

对于 DEM 颗粒，公式为：

I_{p} \frac{d ω_{p}}{d t} = M_{b} + M_{u} + M_{c}

(2991)

且：

M_{c} = \sum_{c o n t a c t s} (r_{c} \times F_{c m} + M_{c m})

(2992)

其中：

$M_{c}$ is the total moment from contact forces, both from rolling resistance models ( $M_{c m}$ ) and any other contact forces ( $r_{c} \times F_{c m}$ )。
$r_{c}$ 为从颗粒重心到接触点的位置矢量。
$M_{c m}$ 为因滚动阻力而作用于颗粒的力矩（滚动阻力）。

曳力

Eqn. (2987) 中的曳力定义如下：

F_{d} = \frac{1}{2} C_{d} ρ A_{p} | v_{s} | v_{s}

(2993)

其中：

$C_{d}$ 为颗粒的曳力系数。
$ρ$ 为连续相的密度。
$v_{s} = v - v_{p}$ 为颗粒滑移速度， $v$ 为连续相的瞬时速度。
$A_{p}$ 为颗粒的投影面积。

Eqn. (2993) 也可以编写为：

F_{d} = \frac{m_{p} v_{s}}{τ_{v}}

(2994)

其中， $τ_{v}$ 是动量松弛时间尺度：

τ_{v} = \frac{2 m_{p}}{C_{d} ρ A_{p} | v_{s} |}

(2995)

曳力系数

Eqn. (2993) 中的曳力系数

C_{d}

是围绕单个颗粒的小尺度流体特征的函数。

这些特征在空间上求解不切实际，所以一般的做法是从相关性（通常从实验或理论研究衍生）获取曳力系数。这些相关性取决于离散相的性质。液滴、气泡和固体颗粒具有不同的相关性。颗粒的形状、相间质量的存在和能量传递等影响这些相关性。

Simcenter STAR-CCM+ 提供了以下方法来定义曳力系数：

Schiller-Naumann 相关性
刘氏动态曳力系数
Haider-Levenspiel

以下曳力系数可用于 DEM 曳力模型：

Clift
Di Felice
Gidaspow

Schiller-Naumann 相关性

Schiller-Naumann 相关性 [703] 适用于球形固体颗粒、液滴和小直径气泡。其设定为：

C_{d} = {\begin{matrix} \frac{24}{{Re}_{p}} (1 + 0.15 {Re}_{p}^{0.687}) & {Re}_{p} \leq 10^{3} \\ 0.44 & {Re}_{p} > 10^{3} \end{matrix}

(2996)

其中， ${Re}_{p}$ 是颗粒雷诺数，定义如下：

{Re}_{p} \equiv \frac{ρ | v_{s} | D_{p}}{μ}

(2997)

其中 $D_{p}$ 是颗粒直径， $μ$ 是动力粘度。仅当连续相为粘性时，此相关性才可用。

刘氏动态曳力系数

刘氏动态曳力系数 [679] 考虑液滴的曳力对其在空气动力作用下变形的依赖性。作为基础，它使用以下表达式计算未变形球体的曳力系数：

C_{d, s p h e r e} = {\begin{matrix} \frac{24}{{Re}_{p}} (1 + \frac{1}{6} {Re}_{p}^{2 / 3}) & {Re}_{p} \leq 1000 \\ 0.424 & {Re}_{p} > 1000 \end{matrix}

(2998)

随着液滴的变形不断增大，假设其形状变为一个盘体，该盘体的轴与滑移速度对齐。这种形状会增大对液滴的曳力。刘氏曳力系数通过指出盘体的曳力系数的雷诺数上限为 1.54 来对此效应建模。随后，假设盘体曳力比所有雷诺数下的球体曳力高 1.54/0.424，并且可以在这两个极限之间插入中间形状的曳力。所以：

C_{d} = C_{d, s p h e r e} (1 + 2.632 y)

(2999)

对于球体，插值因子 $y$ 为 0；对于盘体，则为 1。将它标识为 TAB 变形 Eqn. (3146)，因此刘氏曳力系数需要 TAB 变形模型来计算此物理量。

Clift 法

此方法仅适用于圆柱颗粒。曳力系数

C_{d}

取决于圆柱体的颗粒雷诺数和长宽比。Clift 曳力系数如下：

对于长圆柱颗粒（长宽比 > 1）：

$\begin{matrix} \begin{matrix} \begin{matrix} {Re}_{p} \leq 5 : C_{d} = 9.689 {Re}_{p}^{- 0.78} (1 + 0.147 {Re}_{p}^{0.82}) \end{matrix} \\ 5 < {Re}_{p} \leq 40 : C_{d} = 9.689 {Re}_{p}^{- 0.78} (1 + 0.227 {Re}_{p}^{0.55}) \end{matrix} \\ \begin{matrix} 40 < {Re}_{p} \leq 400 : C_{d} = 9.689 {Re}_{p}^{- 0.78} (1 + 0.0838 {Re}_{p}^{0.82}) \\ 400 < {Re}_{p} : C_{d} = 0.99 γ^{- 0.12} E^{- 0.08} \end{matrix} \end{matrix}$

(3000)

其中 $γ = ρ_{p} / ρ$ 是颗粒密度与流体密度之比，而 $E = h / d$ 是长宽比。
对于扁平圆柱颗粒（长宽比 < 1）：

$\begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} {Re}_{p} \leq 0.01 : C_{d} = \frac{64}{π {Re}_{p}} (1 + \frac{{Re}_{p}}{2 π}) \\ 0.01 < {Re}_{p} \leq 1.5 : C_{d} = \frac{64}{π {Re}_{p}} (1 + 10^{x}) \end{matrix} \\ (x = - 0.883 + 0.906 log ({Re}_{p}) - 0.025 {(log {Re}_{p})}^{2}) \end{matrix} \\ 1.5 < {Re}_{p} \leq 133 : C_{d} = \frac{64}{π {Re}_{p}} (1 + 0.138 {Re}_{p}^{0.792}) \end{matrix} \\ 133 < {Re}_{p} : C_{d} = 1.17 \end{matrix}$
(3001)

Di Felice 曳力系数法

此曳力法在流体曳力表达式 (Eqn. (2993)) 中引入了一个额外项来考虑颗粒周围存在其他颗粒。Di Felice 曳力系数如下：

C_{d} = {(0.63 + \frac{4.80}{\sqrt{ϵ_{i} {Re}_{p}}})}^{2} ϵ^{2 - ζ}

(3002)

其中， ${Re}_{p}$ 定义颗粒雷诺数， $ϵ_{i}$ 是颗粒周围的空隙率，并且：

ζ = 3.7 - 0.65 exp [- 0.5 {(1.5 - log [ϵ {Re}_{p}])}^{2}]

(3003)

项 $ϵ^{2 - ζ}$ 考虑因颗粒周围存在其他颗粒而使自身承受的曳力增强产生的效应。

Gidaspow 曳力系数法

Gidaspow 模型在包含颗粒密度较高的区域的系统中（例如，在流化床中）很有用。该模型通过将 Wen-Yu 与 Ergun 方程 [730] 相结合，涵盖各种不同的空隙率。

Simcenter STAR-CCM+ 包含曳力系数形式的方程，如下所示：

C_{d} = \frac{4}{3} (150 \frac{1 - ϵ}{ϵ {Re}_{p}} + 1.75) if ϵ < ϵ_{m i n}

(3004)

否则：

C_{d} = \begin{matrix}  \end{matrix} {\begin{matrix} \frac{24 (1 + 0.15 \cdot {Re}_{p}^{0.687})}{ϵ {Re}_{p}} ϵ^{1 - w} & if ϵ {Re}_{p} \leq 1000 \\ 0.44 & if ϵ {Re}_{p} > 1000 \end{matrix}

(3005)

其中：

$ϵ$ 为空隙率。
$ϵ_{m i n}$ 为截止空隙率。
$w$ 为 Wen-Yu 指数。

EMMS 曳力系数法

EMMS 模型在包含颗粒密度较高的区域的系统中（例如，在流化床中）很有用。该模型适用于颗粒集群的异质性影响变得显著的模拟。EMMS 曳力法通过考虑中尺度效应和各种相互作用尺度下颗粒交换的能量最小化方法，对流化床内的相互作用进行建模。

对于非均质流化床，EMMS 曳力系数 ([573]) 为：

C_{d} = \frac{3}{4} \frac{α_{p} ϵ ρ_{g}}{d_{p}} C_{D 0} | v_{p} - v | ϵ^{- 2.65} H_{D}

(3006)

其中：

$d_{p}$ 为颗粒直径。
$v_{p}$ 和 $v_{f}$ 分别为颗粒速度和流体速度。
$α_{p}$ 和 $ϵ$ 分别为颗粒的体积分数和空隙率。
$C_{D 0}$ 为单个颗粒的曳力系数，根据 Eqn. (2996) 中 Schiller 和 Naumann 相关性计算得出。
$H_{D} = f (α_{g})$ 为异质性指数，表示异质系统中颗粒所受的平均曳力与均质系统中的曳力之比。在 Simcenter STAR-CCM+ 中，Lu、Wang 和 Li 相关性 ([511]) 用于定义 Geldart A 和 B 颗粒组以及气泡、快速和气动流化区域的异质性指数。流区域根据系统中的气体速度和生成的流化特性来确定 ([467])。Geldart 分类系统根据颗粒的尺寸、密度和流化特性对颗粒进行分组。

Haider-Levenspiel 曳力系数法

曳力系数 $C_{d}$ 取决于颗粒雷诺数和颗粒球度：

C_{d} = \frac{24}{{Re}_{p}} (1 + A {Re}_{p}^{B}) + \frac{C}{(1 + \frac{D}{{Re}_{p}})}

(3007)

且：

\begin{matrix} A & = & 8.1716 e^{- 4.0665 ϕ} \\ B & = & 0.0964 + 0.5565 ϕ \\ C & = & 73.690 e^{- 5.0746 ϕ} \\ D & = & 5.3780 e^{6.2122 ϕ} \end{matrix}

(3008)

其中， $ϕ$ 是颗粒球度。

压力梯度力

压力梯度力定义如下：

F_{p} = - V_{p} \nabla p_{s t a t i c}

(3009)

其中， $V_{p}$ 为颗粒体积， $\nabla p_{s t a t i c}$ 为连续相中的静压力梯度。

虚拟质量力

虚拟质量力 [652] 如下：

F_{v m} = C_{v m} ρ V_{p} (\frac{D v}{D t} - \frac{d v_{p}}{d t})

(3010)

其中， $C_{v m}$ 是虚拟质量系数。此系数的默认值 0.5 针对均匀、无粘性、不可压缩流 [681] 中的球体。运算符 D/Dt 是材料（实质性）导数。

升力

DEM 模拟中的升力可能源自颗粒自旋和/或颗粒剪切。在这种情况下，升力指的是垂直于颗粒速度的力；它们不一定是向上的力。统计拉格朗日模拟中的升力只能源自颗粒剪切。

颗粒旋转升力

此力将应用于相对于流体移动的旋转颗粒。力由以下公式给出：

F_{L R} = \frac{ρ π}{8} D_{p}^{2} C_{L R} | v_{s} | \frac{Ω \times v_{s}}{| Ω |}

(3011)

Sommerfeld [100] 给定旋转升力系数 $C_{L R}$ 为：

C_{L R} = 0.45 + (\frac{{Re}_{R}}{{Re}_{p}} - 0.45) e^{- 0.05684 {Re}_{R}^{0.4} {Re}_{p}^{0.3}}

(3012)

Sommerfeld 指出， $C_{L R}$ 的这一公式对雷诺数低于 140 的颗粒有效。还可以使用场函数定义 $C_{L R}$ 。

颗粒剪切升力

此力将应用于相对于流体移动的颗粒，前提是流体中的速度梯度与相对运动正交。

Saffman [740] 给定剪切升力为：

F_{L S} = C_{L S} \frac{ρ π}{8} D_{p}^{3} (v_{S} \times ω)

(3013)

其中：

$C_{L S}$ 为剪切升力系数。
$ρ$ 为颗粒密度。
$v_{S} = v - v_{p}$ 为滑移速度，流体速度与颗粒速度之差。
$ω = \nabla \times v$ ，流体速度的旋度。

可以将 $C_{L S}$ 定义为场函数，或者使用发布的两个定义中的任意一个：

Saffman [740] 提供的定义可重新获得低雷诺数的原始 Saffman 渐进求解：

$C_{L S} = \frac{4.1126}{{Re}_{S}^{0.5}}$

(3014)
Sommerfeld [100] 提供的定义适用于更广泛的雷诺数：

$C_{L S} = \frac{4.1126}{{Re}_{S}^{0.5}} f {Re}_{p}, {Re}_{S}$

(3015)

其中：

$f {Re}_{p}, {Re}_{S} = {\begin{matrix} (1 - 0.3314 β^{0.5}) e^{- 0.1 {Re}_{p}} + 0.3314 β^{0.5}, & ({Re}_{p} \leq 40) \\ 0.0524 {(β {Re}_{p})}^{0.5}, & ({Re}_{p} > 40) \end{matrix}$

(3016)

和

$β = 0.5 \frac{{Re}_{S}}{{Re}_{p}}$

(3017)

其中：

${Re}_{S}$ 为剪切流的雷诺数：

${Re}_{S} = \frac{ρ D_{p}^{2} | ω |}{μ}$

重力

重力定义如下：

F_{g} = m_{p} g

(3018)

其中， $g$ 是重力加速度矢量。

移动参考坐标系中的力

旋转参考坐标系引入了一个体积力：

F_{M R F} = m_{p} [ω \times (ω \times r) + 2 (ω \times v_{p})]

(3019)

其中， $ω$ 是旋转参考坐标系的角速度矢量， $r$ 是到旋转轴的距离矢量。

库仑力

颗粒可以进入承受电荷的域，或在某些情况下，可以在进入之后获得电荷，从而导致额外的力作用于颗粒。带电颗粒出现在很多应用中，例如，空气污染控制领域的静电喷涂和静电除尘。

在电场 $E$ 中带电荷 $q$ 的带电颗粒将承受一个力。此力称为库仑力，由以下公式给出：

F_{Co} = q E

(3020)

场充电

场充电通过电场对来自载体相中的单极电荷分布的颗粒充电建模。对于大于 1 微米的颗粒，它是颗粒充电的主导机制。

当单极离子（正在沿电场运动）冲击颗粒并将其电荷传给颗粒时会发生场充电，直到颗粒上的电荷达到饱和水平。

暴露于电场的球形颗粒的饱和电荷为：

Q_{s} = \frac{3 π ϵ_{c} ϵ_{r}}{ϵ_{r} + 2} | E | D_{p}^{2}

(3021)

其中：

$ϵ_{c}$ 为介质的介电常数。
$ϵ_{r}$ 为颗粒的相对介电常数。
$E$ 为电场强度，根据 Maxwell 方程计算。

颗粒上的总表面电荷 $Q_{p}$ 演算如下：

\frac{ⅆ Q_{p}}{ⅆ t} = \frac{Q_{s}}{τ_{q}} {(1 - \frac{Q_{p}}{Q_{s}})}^{2}

(3022)

其中，时间常数 $τ_{q}$ 为：

τ_{q} = \frac{4 ϵ_{c}}{ρ K}

(3023)

$ρ$ 为电荷密度。
$K$ 为离子迁移率。

用户自定义的体积力

Eqn. (2988) 中用户自定义的体积力定义如下：

F_{u} = V_{p} f_{u}

(3024)

其中， $f_{u}$ 是用户指定的单位体积的体积力。

曳力扭矩

对于拉格朗日颗粒，曳力扭矩 $M_{b}$ （请参见 Eqn. (2991)）可减小颗粒与其浸入的流体之间旋转速度上的差异。

应用于颗粒的曳力扭矩 $M_{b}$ 定义为 [100]：

M_{b} = \frac{ρ}{2} {(\frac{D_{p}}{2})}^{5} C_{R} | Ω | Ω

(3025)

其中：

$C_{R}$ 为旋转曳力系数。
$Ω$ 为相对于流体的颗粒角速度。

相对于流体的颗粒角速度（或滑移旋转）由以下公式给出：

Ω = \frac{1}{2} \nabla \times v - ω_{p}

(3026)

其中：

$v$ 为流体速度。
$ω_{p}$ 为颗粒的角速度。

Eqn. (3025) 中的旋转曳力系数 $C_{R}$ 定义为 [100]：

C_{R} = {\begin{matrix} \frac{12.9}{{Re}_{R}^{0.5}} + \frac{128.4}{{Re}_{R}}, & 32 \leq {Re}_{R} < 1000 \\ \frac{64 π}{{Re}_{R}}, & {Re}_{R} < 32 \end{matrix}

(3027)

${Re}_{R}$ 为旋转雷诺数，定义如下：

{Re}_{R} = \frac{ρ D_{p}^{2} | Ω |}{μ}

(3028)

尽管 Sommerfeld 将旋转雷诺数限制为低于 1000 的值，但 Simcenter STAR-CCM+ 并未设置其上限。

用户自定义扭矩

Eqn. (2990) 中的用户自定义扭矩定义为：

M_{u} = V_{p} m_{u}

(3029)

其中

m_{u}

为用户指定的单位体积扭矩。

颗粒位置和速度

由于上面章节中列举的力，描述颗粒状态的变量以时间和空间的函数形式发生变化。在拉格朗日架构中，状态基于每个颗粒演变，而欧拉场处于冻结状态。可通过随时间对颗粒动量方程 Eqn. (2986) 和位置方程进行数值积分得出颗粒位置和颗粒速度：

\frac{d r_{p}}{d t} = v_{p}

(3030)

其中， $r_{p} (t)$ 是瞬时位置矢量。颗粒轨线是通过对位置方程进行积分得出的点的轨迹。

在统计拉格朗日方法中，会追踪单个颗粒的粒子束。“粒子束”一词用于表示全部具有相同的状态且以相同的方式演变的一组颗粒。对粒子束中的所有颗粒应用状态的单积分，与欧拉相进行任何交换都乘以粒子束所表示的颗粒数。

在后续章节中，将对于粒子束及其颗粒完全相同的物理量（如上一段中的速度）描述为“颗粒”物理量并为其指定下标 p。“粒子束”一词和下标 $π$ 为专门针对粒子束定义的物理量保留。

粒子束报告的连续相物理量是以网格单元为中心的值，但连续相速度将插值到粒子束位置。单元速度是插补的唯一连续相值。