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

当你第一次尝试使用 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 具有时间单位，并且通常被称为拉格朗日时间尺度（Lagrangian time scale） 或粒子速度响应时间（particle velocity response time），下面我们将解释其原因。

进一步简化方程，假定周围的流体是静止的（û=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秒）为 1 毫秒。如果初始时步仍然远远大于 τ_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 版本开始可用）。

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