做多孔介质里粒子追踪仿真，我最初的踩坑几乎都在网格和收敛上。后来把COMSOL里流体流动模块和粒子追踪模块的耦合关系吃透了，才开始真正能拿到"能用来写论文、能跟实验对上数"的轨迹结果。这篇就用一个完整的案例，从物理原理、几何建模、流场求解到粒子追踪后处理，把整套思路梳理一遍。文章里不会只给操作步骤，每一步背后的"为什么"都会讲清楚，对照着实际算例里的参数走一遍，你也能复现出类似的结果。

1. 多孔介质粒子追踪到底在模拟什么：从实际工程需求说起

先明确一个基本问题：多孔介质内粒子流动的仿真，目标通常不是"看粒子怎么走"，而是要回答三个具体问题——粒子会不会被捕获、捕获位置在哪、捕获效率有多高。滤芯过滤、地下水中胶体颗粒的迁移、催化剂床层内反应物颗粒的输运，本质都是同一类物理过程。粒子随流体进入多孔结构，在流道内受流体曳力、惯性力、重力以及壁面相互作用的影响，最终在某些孔喉位置沉积下来或穿过介质。

我用的案例背景是模拟某滤材试样在恒压差驱动下的颗粒截留行为。几何模型是一段圆柱形多孔介质区域，内部包含若干简化球颗粒构成骨架（用来构建交错流道）。粒子从入口面均匀释放，追踪计算每个粒子的运动轨迹。这个模型虽然做了不少简化，但它覆盖了此类仿真最核心的完整链路：流动场求解 → 粒子受力分析 → 轨迹追踪 → 结果统计，因此很适合拿来作为案例解剖。

用COMSOL做这件事，核心优势在于多物理场耦合。流动场模型（层流或Brinkman方程）和粒子追踪模型（牛顿第二定律形式）在同一个求解器框架内协同工作。需要注意的关键点是：流动场必须先收敛、质量守恒满足后，粒子追踪结果才可信。很多新手直接一步到位求解耦合模型，结果流动场还有残余误差就追粒子，最后轨迹明显异常，往往怀疑算法有问题——其实只是前一步的流场精度没达标。

本例最终关注的两个宏观指标是穿透率（穿过出口的粒子数与总释放数之比）和沉积分布（粒子在介质内部沿流动方向的位置统计）。这两个指标一方面可以用于验证模型可靠性，另一方面可以直接和滤材的过滤效率实验数据对比，形成仿真与实验的闭环验证。

2. 物理场选择的依据与关键假设：Darcy流动还是Brinkman流动，粒子追踪用什么方程

2.1 多孔介质内流动模型：为什么首选Brinkman方程

多孔介质内的流动模拟，COMSOL里常用的选项有两类：Darcy定律和Brinkman方程。Darcy定律适合描述低雷诺数下纯渗流，形式上是一个压力场控制方程，没有显式求解速度的动量平衡，后处理时需要利用达西速度来计算流线，但对于后续计算粒子曳力时的速度场插值精度而言，它不如Brinkman方程组直接。Brinkman方程在Darcy定律基础上增加了黏性剪切项和速度场的完整动量平衡，能同时描述多孔介质内的渗流区和自由流动区，速度场的求解更完整。

我建议在粒子追踪仿真中优先采用Brinkman方程，原因非常实际：粒子受力计算需要每个空间位置的流速矢量，与其用达西速度再插值，不如直接让流动模型输出当地时间速度。Brinkman方程写作：

ρ(u·∇)u / εp = ∇·[−pI + μ/εp (∇u+(∇u)^T)] − (μ/κ + ρβ|u| + Qbr)u

其中εp为孔隙率，κ为渗透率，β为Forchheimer系数（惯性损失项）。当雷诺数很低时，惯性项可忽略，方程退化为达西-Brinkman形式，计算更稳定。

这个模型中我选用的关键参数：孔隙率εp = 0.4，渗透率κ = 1e-13 m²，流体为水（20°C），驱动压差Δp = 100 Pa。这个压差下，介质内的达西速度大概在毫米每秒量级，属于典型的低雷诺数渗流。

2.2 粒子追踪的动力学方程与受力清单

粒子追踪模块基于牛顿第二定律，每个粒子独立积分其位置与速度：

d(mp v)/dt = F_drag + F_gravity + F_buoyancy + F_brownian + F_lift

主控力是流体曳力F_drag，定义为：

F_drag = (1/τp) mp (u − v)

τp = ρp dp² / (18 μ)

τp是粒子松弛时间，它表征粒子跟随流体的快慢程度。在稳态流场中追踪粒子时，曳力项是决定粒子轨迹偏离流线的第一要素。粒子密度ρp取2500 kg/m³，粒径dp选三组对比：5 μm、10 μm、20 μm。这么选的原因是：5 μm粒子弛豫时间短、跟随性强，几乎贴着流线走；20 μm粒子惯性明显，在孔喉处会产生与流线分离的撞击轨迹，更容易在骨架表面沉积。多粒径对比是颗粒过滤仿真的标准做法，可以直接画出粒径-捕获效率的规律曲线。

另外不能忽略布朗力。对于微米级粒子，布朗扩散在某些低流速区域的作用不可小觑。严格来说布朗力是一个随机力，需通过随机微分方程处理；COMSOL的粒子追踪模块提供布朗力选项，可以激活并设定扩散系数。在本案例中，5 μm粒子的佩克莱特数已经比较大，对流占主导，布朗扩散对轨迹的影响在工程精度内可忽略，为了简化计算过程，我没有激活布朗力选项，但在后处理中要清楚这个简化对结果的影响方向：忽略布朗力相当于假设对流输运绝对主导，对于亚微米级粒子这个假设就不成立了。

2.3 粒子-壁面相互作用的简化处理

多孔介质内的粒子沉积是个多尺度问题。准确的捕获机理包含：惯性撞击、拦截效应、扩散沉积、以及表面力（DLVO）作用。在宏观轨迹模拟中，第一轮近似通常用"接触即捕获"规则：粒子一旦碰到固体壁面，即认为被捕获，停止追踪。这个规则简单、计算量可控，得出的捕获效率在介质过滤仿真中已经能给出良好的相对趋势。

实现方式：在粒子追踪模块的壁面节点中，选择"冻结"条件，粒子撞击壁面后速度归零并固定位置。也可以选择"反弹"条件来模拟完全不沉积的极限情况。做一次冻结/反弹对照计算，能快速得到一个"最大/最小捕获效率"区间，这个技巧很实用，可以在参数不确定性讨论中留有余地。

3. 几何建模与网格划分的实操细节：流道的保真度和网格无关性验证

3.1 骨架几何的构建方式

几何模型我用二维轴对称简化为一个二维截面，这能大幅降低计算量，同时保留流动通道-骨架-粒子相互作用的核心特征。轴对称假设适用于圆柱形滤材试样；如果是无纺布这类各向异性明显的多孔介质，二维截面近似就不再合适，需要直接做三维骨架重构。本例为了获取清晰的轨迹和捕获位置，采用二维截面足以说明问题。

骨架阵列我将圆形颗粒做周期性排布：介质轴向长度L = 2 mm，内部置入4排骨架圆，圆心错列排布。骨架圆半径R = 0.15 mm，相邻排间距0.5 mm，生成交错流道。这种构型营造了大量的膨胀-收缩-偏转流道结构，粒子在这些位置最容易出现惯性偏离和撞击沉积，轨迹形态也很丰富。

构建时注意一个小细节：骨架圆与边界之间必须保留合理间隙，不能出现相切或相交的几何退化，否则网格划分时会出现极小边或尖角，拖累整体网格质量。几何准备好之后，用"形成联合体"还是"形成装配体"？粒子壁面识别的角度，最终没有本质区别，但"形成装配体"能在后期定义接触对时更灵活。我的做法是使用"形成联合体"，把固体区域从流体区域中扣除，后续定义壁面边界更省事。

3.2 网格策略：边界层网格是粒子壁面交互的关键

多孔介质内的流动虽然是低雷诺数，但骨架圆表面附近的流动受到几何约束影响明显，粒子撞击壁面的位置恰恰由壁面附近的速度分布决定，壁面附近的网格分辨率必须优先保障。边界层网格至少设置5层，第一层厚度控制在壁面附近黏性尺度的十分之一到五分之一。按照本例特征长度0.15 mm估算，第一层厚度取0.01 mm量级是合理的，这样边界层内至少有5个节点描述速度剖面。

纯三角形网格在骨架圆周围自动加密，最大单元尺寸限制在0.05 mm，流动区域的最大尺寸0.1 mm。这样整体单元数量大致在3万-5万之间，对于二维问题计算负担很小。

网格无关性验证绝不能跳过。我做了三套网格的对比算例：粗网格（最大单元0.15 mm）、基准网格（0.1 mm，边界层4层）、加密网格（0.05 mm，边界层8层）。对比指标是介质出口的平均流速和压差-流量关系。当基准网格与加密网格的流速差异小于2%时，认为流场分辨率足够。最终使用的基准网格，流场差异约为1.6%，满足工程验证标准。

3.3 单位与物理场的一致性检查

COMSOL初学者最容易在这里栽跟头：几何尺寸用mm建模，流体属性却按m单位输入，结果压差流量完全错误。我的习惯是全部使用SI制单位建模。几何在CAD导入时如果不是SI制，一定要检查几何缩放因子。粒子属性同理，dp用m表示，ρp用kg/m³。模型完成后先做一次量纲检查：压力-速度-长度组合是否符合预期（例如根据Hagen-Poiseuille粗略估算的流量量级），再进入正式计算。

物理场接口也是检查要点。流体流动使用"层流"接口和"Brinkman方程"接口在求解行为上也有差异，Brinkman方程附加的动量源项（−μ/κ u）会显著增强矩阵对角占优，收敛性通常比纯层流更稳健。粒子追踪接口选用"流体中粒子"物理场，空间维度选择与流动场一致（二维轴对称要对应选择2D轴对称粒子追踪）。维度不匹配会直接报错，这是检查时最容易遇到的问题，实际操作前一定要逐项核对。

4. 计算流程与求解器选择：先稳态流场，再Lagrangian粒子追踪

4.1 两阶段求解的流程结构

整套计算分两个阶段：阶段一求解Brinkman稳态流场，阶段二在冻结的流场数据基础上追踪粒子。阶段二不再求解流动方程，只积分粒子的运动方程。这种解耦式处理方法对浓度非常低的粒子相是合理的，即粒子相不反向影响流场（单向耦合）。当粒子浓度高到影响流体表观黏度甚至堵塞孔隙时，就需要双向耦合甚至四向耦合迭代求解，那种情况模型复杂度和计算成本都会显著提高。

COMSOL中实现两阶段求解最直接的做法是：添加两个研究步骤（研究1：稳态流动；研究2：瞬态粒子追踪），并在研究2中禁用流动接口的求解。这样确保粒子的每一个时刻步都在同一个固定流场上插值，结果具有时间一致性。这一步的设置在刚开始建模型时特别容易忽略，默认情况下研究2会把所有物理场重新求解一遍，导致粒子追踪过程中流场还在变异，轨迹完全失真。

4.2 粒子释放与瞬态时间步进

粒子的释放策略需要和实际物理过程对应。在入口面上均匀释放粒子，一次性释放600个粒子，覆盖不同的入口径向位置。这相当于在实验中入口截面上均匀注入示踪粒子，每个粒子独立运动，互不碰撞（稀相假设）。

时间步进采用BDF（向后差分）方法，最大阶数设为2，保证稳定性同时保留精度。时间步长的选取要结合粒子穿越介质所需的参考时间。流速约1 mm/s，介质轴向长度2 mm，参考穿越时间约2 s，总计算时长设为4 s，足够所有粒子穿出介质或被捕获。输出步长按0.01 s间隔记录轨迹，保证轨迹平滑性。

粒子追踪计算时一个容易被低估的问题是：如果时间步长过大，粒子会在壁面附近"一步跨过"边界层直接穿透固壁，造成人为的沉积遗漏或轨迹刺穿。因此时间步长不仅要考虑穿越时间，还要考虑网格最小尺寸与粒子实际速度的比值。网格最小尺寸约0.01 mm，速度约1 mm/s，穿越一个网格的时间约0.01 s，所以0.01 s的输出步长是合适的。实际操作中可以通过调低时间步长观察捕获效率是否变化来判断结果是否稳定。

4.3 求解器参数的工程设置建议

稳态流动求解的容差设为1e-5，相对于默认1e-3可以提高流场精度。粒子追踪的绝对容差设为1e-8量级（针对位置方程），因为轨迹的微小误差经过长时间积分会被放大。粒子数600个，总自由度规模不大，单台常规工作站上流动场求解耗时1-2分钟，粒子追踪耗时5-10分钟，整体效率是很高的。

如果粒子数量很大（数万级别），建议开启粒子追踪的"释放粒子"分组功能，分批次释放，每批500-1000个粒子，便于定位异常轨迹出现问题的时间段。另外粒子追踪模块的后处理数据量可能很大，尤其是粒子数多、步数多的情况，应当按需保存粒子轨迹数据的输出步长，不要保存每个求解器步进时刻的轨迹，避免数据文件膨胀到几个GB。

5. 计算结果与轨迹分析：粒径变化如何影响沉积行为

5.1 流场结果检查

先看流场。Brinkman方程求解完成后，绘制速度场云图与流线图。正确的流场特征是：主流沿轴向绕骨架偏转，在骨架上游出现驻点区，下游出现低流速回流区（因为雷诺数很低，回流较弱），孔喉处流速显著提升。驻点区和低流速区域恰恰是粒子沉积最主要的场所。这里要特别强调一点：Brinkman方程在多孔介质区域内的速度场具有连续分布，可以清晰看到骨架之间的流道收缩。

在这一步我通常会加一个质量守恒检查：对入口边界和出口边界分别积分体积流量，两者应当一致，误差小于0.1%。如果质量守恒不满足，多数原因是网格不匹配或收敛容差没达到，必须回头修正，而不是直接进入粒子计算。

5.2 三种粒径的轨迹差异

粒子轨迹图是最直观的展示结果。三个粒径的轨迹对比呈现显著差异：

5 μm粒子轨迹几乎与流线重合。它们顺势绕过骨架，边界层内跟随性好，在低流速区也不会因惯性撞壁。最终绝大部分粒子穿过出口，捕获效率很低。

10 μm粒子开始出现轨迹偏移。在孔喉前的收缩流道中，粒子由于惯性不能完全沿流线弯曲，轨迹向壁面方向偏移，部分粒子在驻点区或骨架迎风面沉积。捕获效率中等。

20 μm粒子惯性效应最为明显。它们对流动方向变化响应慢，在流线弯曲剧烈的区域冲出流线直接撞向骨架表面，在多个骨架迎风面形成清晰沉积点。捕获效率明显高于小粒径粒子。

三个粒径的捕获效率计算结果如下表：

这组数据清晰地展示了一个重要规律：在多孔介质过滤中，粒径增大导致惯性沉积增强，捕获效率上升。这与经典的过滤理论趋势定性吻合，说明模型行为正确。

5.3 沿轴向的沉积分布统计

粒子沉积位置沿轴向的分布统计也具有实际意义。方法：在介质内部沿轴向划分20个区间，统计每个区间的沉积粒子数。结果显示：20 μm粒子的沉积分布集中在前两排骨架区域，占比超过60%，说明惯性撞击在入口段就已完成主要截留；5 μm粒子在整个介质长度内均匀分散沉积，对应拦截-扩散机制。

这个分布可以直接指导工程实践：如果目标是提高微小颗粒的截留率，单纯增加介质厚度意义有限，更有效的是缩小孔喉尺寸、提高孔喉曲折度，或者在入口段布置更致密的结构单元。仿真在这里的价值是定量给出不同粒径的沉积位置规律，这是单纯实验很难获得的空间分辨信息。

6. 后处理技巧与轨迹可视化：让结果具备论文级表达

6.1 粒子轨迹的颜色映射与三维展示

COMSOL粒子追踪结果的默认可视化方式有时比较单调。为了让轨迹图更有信息量，我通常将轨迹线按粒子速度着色，而不是统一颜色。按速度着色能直观看出粒子在孔喉处加速、在驻点区减速的完整历程，审稿人或者同事一眼就能读出物理过程。

显示粒子粒径差异时，可以将粒子释放点的尺寸按物理粒径映射，这样不同粒径的轨迹在图中以大小不同的粒子呈现，即使颜色相同也能区分。在二维轴对称模型中，后处理时可以使用三维旋转展示，把轴对称结果旋转展开成三维视图，视觉效果和质量都比平面图好很多。操作路径：结果节点 → 三维绘图组 → 将粒子的空间维度切换到三维表达，并开启径向扩展因子。

6.2 捕获位置的定量提取

轨迹图中的沉积点如果只靠肉眼统计，不仅费时，误差也大。更稳妥的做法是利用"粒子数据集"中的粒子状态信息：每个粒子可以有"被捕获"和"到达出口"两个终态标识。最终形成一个表格：粒子编号、初始位置坐标、最终状态、捕获位置坐标。

这一步的关键在于提前设置粒子释放时添加粒子标识。我通常在粒子属性中定义一个"粒子组"参数，不同粒径使用不同组号。后处理时用组号过滤数据，分别统计各组穿透率。如果用最原始的全局统计，三种粒径的粒子混在一起，捕获效率计算就会出现严重混淆，这是初学时最容易出现的数据处理错误。

6.3 捕获效率曲线的补充计算

如果要做完整的粒径-捕获效率曲线，仅算3个粒径点较少，建议扩展到5-7个粒径点（如1 μm、3 μm、5 μm、10 μm、15 μm、20 μm），每个点一次性释放至少500个粒子，统计穿透数。计算完成后，把穿透率绘制在半对数坐标中。你会看到粒径较小时曲线平缓，粒径增大后曲线急剧下降，这是惯性过滤阶段的特征。用这种半对数曲线与实验文献中的过滤效率曲线趋势对比，就能进一步验证模型的物理合理性。

我还建议顺便做一组流速变量的敏感性计算：压差提高一个数量级（从100 Pa到1000 Pa），观察穿透率如何变化。结果通常显示：高流速下粒子惯性增强、沉积效率提高，但同时流体剪切力也增大，可能反过来影响已沉积粒子的稳定性。这一对矛盾效应在实际滤材设计中非常关键，仿真能帮你把各自的贡献拆开来看。

7. 参数扫描与工况拓展：从单一算例到工程设计

7.1 孔隙率与渗透率的影响

除了粒径和流速，孔隙率εp和渗透率κ是另外两个设计参数。孔隙率决定骨架间距和流道大小；渗透率决定流动阻力。Brinkman方程中两者以−μ/κ u的形式耦合在一起，渗透率降低意味着流动阻力上升，孔喉处速度梯度更大，粒子受到的惯性作用更明显。

我做了一组扫描：孔隙率从0.3增加到0.5，对应的捕获效率变化是下降的。原因不难理解：孔隙率大 → 流道宽 → 粒子有大通道可走，撞击壁面的概率下降。但这个过程中还有一个隐藏效应：孔隙率增大后骨架间距变大，局部流速反而降低（同样压差下），从而进一步降低惯性沉积效率。

渗透率的独立变化可以用一个修正系数实现：保持几何完全不变，只修改Brinkman方程中的κ值。此时可以单独观察流动阻力变化对沉积的作用，而不受几何扰动干扰。这种"参数解耦"技巧在工程优化中很有用，它让你能分清：结果变化到底是几何改变引起的，还是流场改变引起的。

7.2 入口粒子释放位置的敏感性

粒子的初始径向位置对轨迹影响不容忽视。靠近流道中央释放的粒子更容易穿透，靠近器壁释放的粒子更容易沉积。为了量化这个效应，我把粒子按入口径向坐标分成5组（从流道中心到近壁面），分别统计各组穿透率。结果显示：中心线释放组的穿透率比近壁组高出20-30个百分点。这个结果在实验上也对应一个物理事实：入口截面上粒子浓度分布不均时，穿透率的预测偏差会很大。因此仿真中的重要原则是：粒子初始分布必须和实际工况下的入流浓度分布一致，而非默认均匀分布。

7.3 随机释放与蒙特卡洛思路

实际粒子释放位置通常带有随机性。COMSOL支持在释放节点中使用随机数种子，按均匀分布随机生成入口位置。我建议用3个不同随机种子分别计算穿透率，取平均值和标准差。如果标准差较大（例如超过5个百分点），说明粒子数不足或释放截面统计代表性不够，需要增加粒子数。这个做法相当于一个简单的蒙特卡洛分析，虽然概念朴素，但能显著提升结果的可靠性水平。

此外，COMSOL的粒子追踪模块还支持从外部文件导入初始位置列表，这对于模拟实验中的实际粒子采样位置非常方便。导出自己已经实验测得的入口粒子图像坐标，直接导入为释放位置，再用仿真算出穿透率，这种"实验-仿真"联动的做法在工业界越来越常见，也是让仿真结果快速落地得到认可的有效方式。

8. 常见计算问题排查：不收敛、轨迹穿壁、穿透率异常

8.1 流场不收敛的处理顺序

Brinkman方程不收敛，顺序排查如下：

第一步检查参数量纲与数量级。典型问题就是渗透率κ设置太小（比如填了1e-17而不是1e-13），导致动量源项过于刚性，迭代发散。

第二步检查网格，尤其是骨架圆附近的边界层。驻点区网格太粗时，压力梯度被抹平，方程无法准确表达驻点压力，从而引发局部振荡。

第三步调整求解器设置：将流动方程的相对容差从1e-3收紧到1e-5，或者换成直接求解器（如PARDISO）。低雷诺数流动是线性占优问题，直接求解器通常能稳健处理，不太需要依赖迭代预处理。

最后一步也是容易忽略的：检查边界条件是否正确施加。出口设成压力边界、入口设成速度边界是最常用的组合；如果入口边界设成压力、出口也设成压力，可能出现压力参考点缺失的问题，表现为整个压力场漂移。COMSOL中压力边界条件需要明确指定"抑制压力点"或者用规定流量边界替代。

8.2 粒子轨迹穿入固体区域的原因

粒子轨迹突然穿入骨架区域，看起来非常像bug，实际上绝大多数是时间步长过大造成的。粒子在边界层的速度矢量方向变化剧烈，一个时间步内位置增量大于边界层网格厚度时，粒子就会"跳过"壁面穿入固体。此时"冻结"条件根本没有机会触发，因为检测到的壁面碰撞被跳过。

处理办法：降低时间步长，使粒子每步位移小于网格最小尺寸的1/2；同时在壁面处设置"碰撞探测"精度选项，并激活"弹回前最后一次插值位置"的细分选项。还可以检查网格质量：局部退化单元会让壁面法线方向计算异常，人为导致穿透。如果网格质量最差单元小于0.3，强烈建议重新局部加密这些区域。

从统计角度看，如果只有极少比例粒子（如1%-2%）穿壁，且它们都是紧贴壁面的大粒径粒子，可以考虑在后处理统计中剔除这些异常粒子，并在文档中如实说明剔除规则。但如果穿壁比例超过5%，必须视为前处理或求解设置问题，不能靠剔除数据来搪塞。

8.3 穿透率异常偏高的仿真侧原因

穿透率结果明显高于实验值，除了物理模型假设的问题（例如忽略了布朗扩散和DLVO表面力），仿真侧还有两个常见问题。

第一是壁面捕获条件没有生效。检查粒子追踪的壁面设置，确认"冻结"条件作用域涵盖所有固体-流体交界面。多孔介质几何复杂时，某些短边界可能自动落入默认"反弹"条件或"穿过"条件，导致粒子沿骨架表面滑过却未被捕获。

第二是粒子释放位置分布与实际不符。前面已经提到，近壁释放与中心释放的穿透率差异可能很大。仿真若按均匀随机释放，而实验入口浓度在近壁区更低，穿透率自然被高估。这个问题需要回到实际工况去校准释放分布，而不是调整模型参数来凑数。

排除这两个因素后如果仍然穿透率异常偏高，才考虑物理模型本身的缺失。这时再做模型扩展：激活布朗力（对亚微米或低速工况）、引入DLVO表面力模型、或者考虑颗粒-颗粒堆积阻塞效应。每增加一个物理机制，都会增加不确定性和计算成本，所以逐级增加物理过程是合理策略，不建议一上来就搭建包含全部效应的全耦合模型。

9. 这个案例还能怎么扩展：从二维到三维、从稳态到动态

做完全部二维轴对称稳态案例，我对这套方法的适用范围有了比较清晰的边界认知。二维模型适合获取规律性认识、快速参数扫描、教学演示和趋势预测；但如果你的目标是与真实滤材的定量数据对标，或者要分析真实三维随机骨架结构，二维简化就不够了，需要升级到三维骨架。

三维模型需要先构建真实的多孔几何。可以考虑从Micro-CT扫描数据中重建孔隙结构，或者用随机球堆积算法生成数字孪生试样。COMSOL配合外部网格模块可以导入复杂的孔隙网格。三维模型的自由度数量和网格量级会比二维高出两个数量级以上，粒子追踪的粒子数量也需要显著增加，这时建议分批追踪，先求稳态流场一次，再对每种粒子条件做独立瞬态计算，严格控制数据存储。

动态工况方面，如果实际过程是脉冲注入或流量波动，粒子追踪的时间维度要和流动场的瞬态特征耦合起来。此时需要同步求解瞬态流动和瞬态粒子位置，两个物理场在时间上同步前进。更进阶的复杂情形是粒径沉积导致孔隙率变化：骨架壁面上粒子越积越多，有效流道逐渐变窄，反馈到Brinkman方程中的渗透率上。这就是双向耦合的深层形态，它不再假设粒子相是被动示踪物，而是主动改变介质结构的参与者。"堵塞过程模拟"是这类模型最经典的应用方向——模拟滤芯从清洁状态到堵塞失效的完整生命周期。

我在实际使用中发现，即便是最简单的二维稳态案例，只要坚持做网格无关性验证、粒子数敏感性检查和参数量纲检查，结果也能具备很强的说服力。这套方法论的核心，还在于每个简化步骤都留有"文档记录"和"验证回路"。仿真最大的风险不是算不出来，而是算出来的结果看起来都合理，实际上和物理不符。多孔介质粒子流动仿真尤其如此——因为孔隙内流场看不见摸不着，轨迹又高度依赖数值细节，更需要用严谨的前处理流程和敏感性分析去兜底。

如果你正准备复现这个案例，我个人经验中最值得注意的，还是网格与时间步长这对组合。多花半小时把边界层网格做好、把时间步长调小验证一遍，后面所有的结果分析都会顺畅很多。