当你第一次使用 COMSOL 软件对流体中直径为几十微米或更小的微小粒子进行粒子追踪仿真时，你可能会注意到，瞬态求解器采用的时步较一般的情况要短得多。这通常是由于粒子的运动方程表现出数值刚度而导致的。这篇文章，我们将介绍与粒子仿真有关的刚度的概念，并对基于粒子大小选择正确的方程提供一些建议。
示例：小球形粒子的重力沉降
以一个小球形粒子为例，当粒子以恒定速度 u（SI 单位：m/s）穿过周围流体时，遵循牛顿第二运动定律，• F_t（SI单位：N）是作用在粒子上的净力或总力对于在流体中沉降的粒子，其受到的总力是重力 F_g 和曳力 F_D 之和，• ρ（SI 单位：kg/m^3）为周围流体的密度• g（SI单位：m/s^2）为重力引起的加速度（海平面以下约为 9.8m/s^2）• μ（SI 单位：Pa s）为周围流体的动力黏度• v（SI 单位：m/s）为粒子的速度（v≡dq/dt）重力表达式中的（ρp-ρ）/ρp 项代表浮力。当粒子（如空气中的固体粒子）比流体重得多时会取代流体，浮力值接近 1。当粒子和周围的流体密度相同时，浮力值接近于零，在这种情况下，粒子被称为悬浮粒子。
此处使用的曳力表达式来自斯托克斯曳力定律（Stokes drag law）。当粒子的相对雷诺数非常小时，适用于曳力定律：假设粒子大小不变（d_p 和 m_p 为常数），球形粒子的质量可以表示为结合方程1–3，我们得到粒子运动方程的简化表达式：τp 具有时间单位，通常被称为拉格朗日时间尺度或粒子速度响应时间，我们将在下文对其进行解释。进一步简化方程，假设周围的流体是静止的（û=0），并且所述粒子初始为静止状态（q=0, v=0 时间 t=0）。假设对齐坐标系，使重力矢量指向–y 方向。然后，根据方程4，粒子位置的 y 分量方程变为给定初始条件 qy=0 和 vy=0，方程5的精确解或解析解为转换为无量纲变量
为了更好的理解粒子在最初的一段时间 τp 内是如何加速的，我们可以用相应的无量纲量(t´, qy´, vy´)来代替时间、位置和速度(t, qy, vy)变量，将其定义为在下图中，将无量纲位置和速度变量绘制为无量纲时间 t 的函数。该绘图表明粒子速度逐渐接近自由沉降速度，粒子加速主要发生在拉格朗日时间尺度 τp 的最初一段时间内。在此初始加速期之后，粒子位置呈线性变化。重力的作用下粒子沉降过程中的无量纲位置和速度绘图，从静止状态开始。
一些典型粒径的时间尺度
为了更好地了解粒子加速所涉及的时间尺度，假设粒子为密度 2200kg/m^3 的石英玻璃珠。下表列出了不同粒径的粒子在空气和水中的一些拉格朗日时间尺度值。τp 与直径平方呈线性关系意味着，大粒子比小粒子具有更长的速度响应时间和更快的自由沉降速度。这会产生两个主要结果：- 当大粒子以一定的初始速度射入流体时，会沿着入射轨迹，并能够在曳力使其减速之前行进相当长的距离。相反，较小的粒子将更快地适应流体速度。粒子散开很可能是由于周围流体的湍流扩散所致。
数值粒子追踪仿真
上文中，方程4具有精确的解析解。之所以能够获得精确解是因为我们做了许多简化，最主要的是流体速度 u 处处为零。在大多数真实情况中，周围流体的速度不仅不为零，而且在空间上也不均匀，因此仅靠公式不可能得到精确解。对于一般的问题，我们可以通过数值仿真获得近似解，其主要思想是：在初始时间 t=0，给定初始粒子位置 q0 和速度 v0，可以使用数值时步算法计算一组离散时步 t1，t2，t3…… 的解。为此，我们设计了各种不同的时步算法，其中有许多在 COMSOL Multiphysics® 软件中可用。用数值方法求解一组微分方程会引入一定量的误差,即实际粒子运动与计算得到的数值解之间的差异。虽然我们通常无法从数值仿真中获得完美的解，但当 t1，t2–t1，t3–t2 等时间间隔减小时，我们可以将模拟的粒子运动变得更加精确作为首要目标。需要权衡的是，如果时步较小，则需要花更多的时步才能达到相同的输出时间。最终，这可能会导致实际运行时间显著增加，这意味着必须等待较长的时间才能完成仿真。数值仿真工程师必须始终在求解精度和时间之间寻求合理的平衡。
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® 软件新用户感到困惑。假设这些粒子被分配了初始速度并指向仿真域。从上文中的截图可以看出，初始时步大小（总仿真时间为1s）为 1ms。如果初始时步仍然远大于 τp，则曳力可能会过度补偿，导致粒子速度短暂改变方向并指向入口边界。如果发生这种情况，粒子可能会错误地检测到与入口边界的碰撞，使它们卡在此处。处理粒子追踪模型中的数值刚度
有两种方法可以求解流体中粒子运动的数值刚度模型，即输出时间间隔比 τp 大几个数量级的模型。第一种是我们所说的“强力”方法：只需告诉求解器采用更小的时步即可。如果不想生成大量的输出（这可能会创建大量文件大小），可以不考虑输出时间，而是在求解器序列的瞬态求解器设置中指定一个较小的时间步长或最大时间步长。另一种方法从 COMSOL Multiphysics® 5.6 版本开始引入，即从方程4中删除惯性项。首先，我们把方程4改写成一对耦合的一阶方程，![]()
然后，仅假设曳力始终与其他作用力处于动态平衡，而不是在 τp 最初一段时间完全解析粒子运动，换句话说，即仅假设粒子立即达到其自由沉降速度。如果达到自由沉降速度所需的时间比总仿真时间小很多数量级，那么这是一个合理的近似值。方程6可用于求解 v，其中，Fother 是除曳力以外的其他所有作用力的总和。然后，我们要做的就是把 v 的表达式对时间进行积分以获得粒子位置 q。请在粒子释放和传播部分选择流体流动颗粒跟踪接口要求解的方程组。从公式列表中，可以选择以下选项之一：
需要注意的是，牛顿和牛顿，一阶公式可用的内置力的数量略多于牛顿，忽略惯性项公式。这些力明显取决于粒子速度或已被排除的其他粒子的相对位置。
![]()
牛顿公式中可用的力。
牛顿，忽略惯性项公式中可用的力。
以下是 COMSOL 官网案例库中使用牛顿，忽略惯性项公式追踪需要较长求解时间的微小粒子的示例：由于粒子足够大导致惯性对粒子运动具有显著影响，以下示例使用牛顿公式追踪：结语
使用流体流动接口的粒子追踪模拟流体中的小颗粒运动时，通常应从估计与粒子相关的拉格朗日时间尺度 τp 开始，如果具有不同粒径的分布，应基于最小粒径进行估算，因为模型中最小惯性粒子决定了运动方程的数值刚度。如果要在比速度响应时间大得多的时间范围内预测粒子运动（比如几千倍甚至更多），则应该考虑惯性是否在粒子运动中实际起着重要作用。如果不是，可以从列表中选择牛顿，忽略惯性项（从 5.6 版本开始可用）。如果仍要考虑惯性，可以使用牛顿或牛顿一阶公式。但是，请注意，要求解的方程组是数值刚性的，可能需要手动减小求解器采取的时步大小，以防止粒子位置和速度发生非物理振荡。本文内容来自 COMSOL 博客，点击底部“阅读原文”，可查看更多延伸文章。如果您有相关问题，或者文中介绍的内容没有涉及您所关注的问题，欢迎留言讨论。
