离散元方法

离散元方法 (DEM) 设计用于对沙粒、食品颗粒、粉末、胶囊和浆料等材料的颗粒流体进行建模。此类流体的特点是颗粒密度高，颗粒之间的相互作用非常重要。

DEM 模型由 Cundall 和 Strack [720] 构建，将拉格朗日公式延伸到在颗粒运动方程中考虑颗粒之间的相互作用。对于高负载流体（即，大量相互作用的颗粒），这些力不容忽视。

DEM 颗粒可以具有不同的形状和体积。Simcenter STAR-CCM+ 提供多种选项用于模拟非球形颗粒（请参见 DEM 颗粒类型），包括使用 CAD 几何定义的任意凸面或凹面形状。非球形颗粒和其他颗粒或壁面之间的碰撞可以有多个力传输接触点。

Simcenter STAR-CCM+ 将根据允许颗粒发生重叠的软颗粒公式对 DEM 颗粒进行建模。计算的接触力与重叠以及颗粒材料和几何属性成比例。此接触力将作为 $F_{c}$ 代入拉格朗日动量方程 Eqn. (2958)。

离散元方法颗粒时间尺度

要获取 DEM 颗粒轨迹，需要使用特征 DEM 颗粒时间尺度对 Eqn. (2958) 进行时间积分。

颗粒时间尺度（以及因此涉及的时间步）取决于与颗粒关联的以下三个时间中的某一个：

瑞利波传播时间 $τ_{1}$

假设在单个时间步期间，作用于颗粒的力仅受颗粒的近邻影响，这将限制 DEM 颗粒允许的最大时间步。因此，时间步受限于瑞利波穿过颗粒表面传播到另一侧所用的时间 [726]（如果是球体，则为从极点到极点）：

图 1. EQUATION_DISPLAY

τ_{1} = π \frac{R}{V_{R a y l e i g h}}

(3244)

其中， $R$ 是用于解释颗粒惯性矩的颗粒最小特征尺寸的一半（如果是球体，则为半径）。

瑞利波速取决于材料属性，通过特征方程 [731] 即可求解出精确值。Simcenter STAR-CCM+ 使用求解的近似值提供足够的精度，而不引起高昂的计算成本 [728]、[729]。

使用线性弹簧接触模型时，

图 2. EQUATION_DISPLAY

τ_{1} = 2 \sqrt{m / k}

(3245)

其中， $m$ 是颗粒质量， $k$ 是弹簧刚度。

最大时间步不得超过 $t_{scale} τ_{1}$ ，其中 $t_{scale}$ 是用户自定义的时间尺度。

冲击持续时间 $τ_{2}$

应用于移动颗粒的第二个时间步限制条件是两个完全弹性颗粒的冲击持续时间。对于两个半径为

R

的相同球体，假设符合 Hertz 接触理论，此持续时间由 Timoshenko [735] 推导得出：

图 3. EQUATION_DISPLAY

τ_{2} = 2.94 {(\frac{5 \sqrt{2} π ρ}{4} \frac{1 - ν^{2}}{E})}^{\frac{2}{5}} \frac{R}{\sqrt[5]{ν_{i m p a c t}}}

(3246)

最大时间步不得超过 $0.1 τ_{2}$ 。

颗粒过渡时间 $τ_{3}$

对 DEM 颗粒时间步的决定性限制是几何。它假设颗粒不可在时间步内移动太远。此条件可防止 DEM 颗粒之间以及颗粒与壁面之间缺少接触。因此，每个颗粒均受到限制，至少需要 10 个时间步才能移动半径全长。

图 4. EQUATION_DISPLAY

τ_{3} = (\frac{R}{ν_{p a r t i c l e}})

(3247)

其中， $R$ 是颗粒最小特征尺寸的一半（如果是球体，则为半径）。

最大时间步不得超过 $0.1 τ_{3}$ 。

最终颗粒时间步确定为

t_{scale} τ_{1}, 0.1 τ_{2}, 0.1 τ_{3}

中的最小值。实际上，

τ_{1}

通常为限制因子，而

τ_{2}

和

τ_{3}

仅约束快速移动的颗粒，或者为了加速模拟，将材料的杨氏模量设为较低的值。

人工粘度

当高速冲击导致数值模型过度预测碰撞对象之间的重叠时，可使用人工粘度模型来校正结果。人工粘度耗散的能量表示未求解的损失、变形过程和熵。

“人工粘度”一词仅指人工粘度模型公式与粘度公式之间的相似度。它与任何类似于真实粘度的物理效应无关。

此模型以 von Neumann-Richtmyer 人工粘度 [738] 为基础，在压缩期间（颗粒重叠且其中心正在接近）引入了力：

图 5. EQUATION_DISPLAY

F_{v} = - f (Q_{2} | v_{n} | - Q_{1} c) \frac{m_{e q}}{R_{e q}} v_{n}

(3248)

其中：

$F_{v}$ 为人工粘性力。
$v_{n}$ 为速度在接触法向上的分量。
$m_{e q}$ 为颗粒的等效（调和平均）质量。
$R_{e q}$ 为颗粒的等效（调和平均）半径。
$c$ 为材料中的声速。它根据等效质量和半径以及正在使用的特定接触模型的接触刚度估计得出。
$Q_{1}$ 为线性系数。
$Q_{2}$ 为二次项系数。

因子 $f$ 为按比例缩放到颗粒尺寸的线性跃升函数：

图 6. EQUATION_DISPLAY

f = {\begin{matrix} 0, & r / R_{e q} < p_{lower} \\ \frac{r / R_{e q} - p_{lower}}{p_{upper} - p_{lower}}, & p_{upper} \geq r / R_{e q} \geq p_{lower} \\ 1, & r / R_{e q} > p_{upper} \end{matrix}

(3249)

其中：

$r$ 为颗粒之间或颗粒与壁面之间的穿透（重叠）距离。
$p_{upper}$ 和 $p_{lower}$ 为跃升上限和下限。它们根据穿透 $r$ 与等效颗粒半径 $R_{e q}$ 之比来定义。

利用此跃升函数 $f$ 可以使用人工粘度，同时不会影响典型的低穿透接触。