模拟流体中的粒子时应该使用哪种公式？

当你第一次使用 COMSOL 软件对流体中直径为几十微米或更小的粒子进行粒子追踪仿真时，你可能会注意到，瞬态求解器采用的时步比一般的情况要短得多。这通常是由于粒子的运动方程表现出 数值刚度 而导致。这篇博客，我们将介绍与粒子仿真有关的刚度的概念，并对基于粒子大小选择正确的方程提供一些建议。

示例：小球形粒子的重力沉降

以一个小球形粒子为例，当粒子以恒定速度 u （SI 单位：m/s）穿过周围流体时，遵循牛顿第二运动定律，

(1)

式中，

• m_p（SI单位：kg）是粒子质量
• q（SI单位：m）是粒子位置矢量
• F_t（SI单位：N）是作用在粒子上的净力或总力

对于在流体中沉降的粒子，其受到的总力是重力 F_g 和曳力 F_D 之和，

(2)

式中，

• ρ_p（SI 单位：kg/m³）为粒子密度
• ρ（SI 单位：kg/m³）为周围流体的密度
• g（SI单位：m/s²）为重力引起的加速度（海平面以下约为 9.8m/s²）
• μ（SI 单位：Pa s）为周围流体的动力黏度
• d_p（SI 单位：m）为粒径
• u（SI 单位：m/s）为周围流体的速度
• v（SI 单位：m/s）为粒子的速度（v≡dq/dt）

重力表达式中的（ρ_p-ρ）/ρ_p 项代表浮力。当粒子（如空气中的固体粒子）比流体重得多时会取代流体，浮力值接近 1。当粒子和周围的流体密度相同时，浮力值接近于零，在这种情况下，粒子被称为悬浮粒子。

此处使用的曳力表达式来自斯托克斯曳力定律（Stokes drag law）。当粒子的相对雷诺数非常小时，适用于曳力定律：

\textrm{Re}_\textrm{r} \equiv \frac{\rho d_\textrm{p}\left|\mathbf{u}-\mathbf{v}\right|}{\mu} \ll 1

此定律对于较小的粒子更有效。

假设粒子大小不变（d_p 和 m_p 为常数），球形粒子的质量可以表示为

(3)

结合方程1–3，我们得到粒子运动方程的简化表达式：

(4)

这里，引入了常数 τ_p，

τ_p 具有时间单位，通常被称为 拉格朗日时间尺度 或 粒子速度响应时间，我们在下文对其进行将解释。

进一步简化方程，假设周围的流体是静止的（û=0），并且所述粒子初始为静止状态（q=0, v=0, 时间 t= 0）。假设对齐坐标系，以使重力矢量指向 –y 方向。然后，根据方程4，粒子位置的 y 分量方程变为

(5)

给定初始条件为 q_y=0 和 v_y=0 时，方程 5 的精确解或解析 解为

式中，v_t 是自由沉降速度，

转换为无量纲变量

为了更好的理解粒子在最初的一段时间 τ_p 内是如何加速的，我们可以用相应的无量纲量(t´, q_y´, v_y´)来代替时间、位置和速度（t, q_y, v_y）变量，将其定义为

将这些无量纲变量代入解析解中，得到

在下图中，无量纲位置和速度被绘制为无量纲时间 t 的函数。该绘图表明粒子速度逐渐接近自由沉降速度，粒子加速主要发生在拉格朗日时间尺度 τ_p 的最初一段时间内。在此初始加速期之后，粒子位置呈线性变化。

重力作用下，粒子沉降过程中的无量纲位置和速度绘图，从静止状态开始。

一些典型粒径的时间尺度

为了更好地了解粒子加速所涉及的时间尺度，假设粒子为密度 2200 kg/m³ 的石英玻璃珠。下表列出了不同粒径的粒子在空气和水中的一些拉格朗日时间尺度值。

τ_p 与直径平方呈线性关系意味着大粒子比小粒子具有更长的速度响应时间和更快的自由沉降速度。这会产生两个主要结果：

1. 大粒子落地的速度比小粒子快得多。
2. 当大粒子以一定的初始速度射入流体时，会沿着入射轨迹，并能够在曳力使其减速之前行进相当长的距离。相反，较小的粒子将更快地适应流体速度。粒子散开很可能是由于周围流体的湍流扩散所致。

数值粒子追踪仿真

上文中，方程4 具有精确的解析解。之所以能够获得精确解是因为我们做了许多简化，最主要的是流体速度 u 处处为零。在大多数真实情况中，周围流体的速度不仅不为零，而且在空间上也不均匀，因此仅靠公式不可能得到精确解。

对于更一般的问题，我们可以通过数值仿真获得近似解，其主要思想是：在初始时间 t=0 时，给定初始粒子位置 q₀ 和速度 v₀，可以使用数值时步算法来计算一组离散时步 t₁，t₂，t₃，……的解。为此，我们设计了各种不同的时步算法，其中有许多是在 COMSOL Multiphysics^® 软件中可用。

用数值方法求解一组微分方程会引入一定量的误差,即实际粒子运动与计算得到的数值解之间的差异。虽然我们通常无法从数值仿真中获得完美的解，但当t₁，t₂–t₁，t₃–t₂ 等时间间隔减小时，可以将模拟的粒子运动变得更加精确作为首要目标。

需要权衡的是，如果时步较小，则需要花更多的时步才能达到相同的输出时间。最终，这可能会导致 实际运行时间 显著增加，这意味着必须等待较长的时间才能完成仿真。数值仿真工程师必须始终在求解精度和时间之间寻求合理的平衡。

COMSOL Multiphysics^® 中的粒子追踪模块的 流体流动颗粒跟踪 接口通过数值求解牛顿第二定律来模拟周围流体中单个粒子的运动。基本上，此接口可求解方程1，同时允许在方程右侧添加各种不同的力。该接口还包括用于设置初始粒子位置和速度，以及检测和处理粒子与周围几何的表面碰撞的各种选项。

处理小粒子和长时间尺度

在许多实际应用中，粒子追踪模型所需求解时间的范围远大于拉格朗日时间尺度 τ_p。假设我们要在 1s 的总仿真时间内追踪直径约 20μm 的石英玻璃颗粒在水中的运动。由上述表格可知，这样的小粒子的拉格朗日响应时间约为 5×10^-5s，所以总仿真时间约为 2000τ_p。如果我们想在几分钟或几小时时间跨度内追踪更小的粒子，那么总仿真时间可能比 τ_p 大几百万倍。

下图显示了瞬态求解器在跟踪这些 20μm 的粒子时所采取的时步日志。在 步骤1 中输出时间的范围：瞬态 节点已被设置为 range(0,0.1,1)，这意味着它将仅以 0.1s 的倍数存储输出。但是，这并不妨碍求解器在必要时采用更短的时步来获得精确的解。如如所示，求解器先从 1ms 或更小的时间步开始，然后在粒子接近其最终速度时逐渐采用更大的时步。

如下图中的步骤 24 所示，在 COMSOL Multiphysics 中，粒子追踪物理场接口通常使用严格的时间步算法，该算法至少要求求解器所采取的某些步长与输出时间一致。但这并不是所有物理场的要求；对于某些物理场接口，可以通过在求解器所采用的最近步长之间进行插值获得输出时间。

在研究接近尾声时，粒子几乎不再加速，所以时步可能会很大。最终，求解器需要 24 个时步才能在 0.1s 获得第一个输出时间，但是只需要再增加 12 个时步就能在 1s 最终求解。

重力沉降中的粒子运动方程是 刚性常微分方程（刚性 ODE广义 α，这是一种二阶隐式时步法，非常适合用于处理刚性问题。如果需要额外的稳定性，可以在瞬态 求解器 设置中调整 放大高频 的数值阻尼项。因此，随着粒子速度接近自由沉降速度，时步变得更大。（显式 Runge–Kutta 方法 RK34 采用了 7425 个时步来求解相同的问题！）

但是，如果粒子在几个不同的释放时间进入仿真域，或者背景流体速度在空间上不均匀（这样粒子在以后的研究中仍会加速），求解器可能直到最终时间都会一直采用如此小的时步。如果我们尝试在很长的仿真时间内追踪非常小的粒子，最终将需要大量时间才能完成这些研究，因为求解器可能需要成千上万甚至数百万时步。

使用 入口 边界条件将粒子释放到仿真域可能会使 COMSOL^® 软件新用户感到困惑。假设这些粒子被分配了初始速度并指向仿真域。从上文的截图可以看出，初始时步大小（总仿真时间为 1 s）为 1 ms。如果初始时步仍然远大于 τ_p，则曳力可能会过度补偿，导致粒子速度短暂改变方向并指向入口 边界。如果发生这种情况，粒子可能会错误地检测到与入口 边界的碰撞，使它们卡在此处。

处理粒子追踪模型中的数值刚度

有两种方法可以求解流体中粒子运动的数值刚度模型，即输出时间间隔比 τ_p 大几个数量级的模型。

第一种是我们所说的“强力”方法：只需告诉求解器采用更小的时步即可。如果不想生成大量的输出（可能会创建大量文件大小），可以不考虑输出时间，而是在求解器序列的瞬态求解器设置中指定一个较小的时间步或最大时间步。

另一种方法从 COMSOL Multiphysics^® 5.6 版本开始引入，即从 方程4 中删除惯性项。首先，我们把方程4 改写成一对耦合的一阶方程，

然后，仅假设曳力始终与其他作用力处于动态平衡，而不是在 τ_p 最初一段时间完全解析粒子运动，

(6)

换句话说，即 仅假设粒子立即达到其自由沉降速度。如果达到自由沉降速度所需的时间比总仿真时间小很多数量级，那么这是一个合理的近似值。方程6 可用于求解 v，

或者，使用更通用的形式，

(7)

其中，F_other 是除曳力以外的其他所有作用力的总和。

然后，我们要做的就是把 v 的表达式对时间进行积分以获得粒子位置 q。

请在 粒子释放和传播 部分，选择 流体流动颗粒跟踪 接口要求解的方程组。从公式列表中，可以选择以下选项之一：

• 牛顿：求解方程1
• 牛顿，一阶：将方程1 分离为 q 和 v 的一对耦合一阶方程，然后求解它们
• 牛顿，忽略惯性项（自版本 5.6 起可用）：使用方程7 定义速度的简化公式，然后求解 q
• 无质量：一种更简化的公式，其中直接指定 v 来求解 q

需要注意的是，牛顿 和 牛顿，一阶公式 可用的内置力的数量略多于 牛顿，忽略惯性项 公式。这些力明显取决于粒子速度或已被排除的其他粒子的相对位置。

牛顿 公式中可用的力。

牛顿，忽略惯性项 公式中可用的力。

下面是 COMSOL 案例库中使用 牛顿，忽略惯性项 公式来追踪长求解时间内的很小的粒子的示例：

• 层流静态混合器中的粒子轨迹
• 使用介电泳从红细胞中分离血小板

由于粒子足够大导致惯性对粒子运动具有显著影响，以下示例使用了 牛顿 公式追踪：

• 微混合器中的颗粒跟踪
• 污染物颗粒造成的管道冲蚀

结语

当使用流体流动接口的粒子追踪模拟流体中的小颗粒运动时，通常应从估计与粒子相关的拉格朗日时间尺度 τ_p 开始，

并将此时间尺度与要模拟的求解时间范围进行比较。

如果具有不同粒径的分布，应基于最小粒径进行此估算，因为模型中最小惯性粒子决定了运动方程的数值刚度。

如果要在比速度响应时间大得多的时间范围内预测粒子运动（比如几千倍甚至更多倍），则应该考虑惯性是否在粒子运动中实际起着重要作用。如果不是，则可以从列表中选择 牛顿，忽略惯性项（从 5.6 版本开始可用）。

如果仍要考虑惯性，可以使用 牛顿 或 牛顿，一阶 公式。但是，请注意，要求解的方程组是数值刚性的，可能需要手动减小求解器采取的时步大小，以防止粒子位置和速度发生非物理振荡。