一种带有动态域映射的结晶群体平衡方程自适应网格方法

《Computers & Chemical Engineering》：An adaptive grid method for crystallization population balance equations with dynamic domain mapping

　　

## 摘要

群体平衡方程(PBEs)被广泛用于预测颗粒尺寸分布(PSD)，但固定域离散化在颗粒数量可忽略的尺段上会浪费分辨率，并且除非使用精细网格，否则可能遭受数值扩散的影响。本研究提出了柔性域群体平衡方程(FlexPBE)，一种自适应结晶PBE框架，其中尺度域通过移动边界随PSD演化。FlexPBE采用既有的动态域映射方法，而非引入新的坐标变换。边界基于由矩导出的PSD统计量使用Z分数准则进行更新，同时采用指数移动平均来防止剧烈运动和数值刚度。该映射修改了尺度导数项，使FlexPBE能够作为传统离散化方案的包装器运行，而无需改变其结构。FlexPBE针对解析批次生长和成核-生长基准进行了验证。成核和晶体生长被视为基本结晶机制，而团聚和破碎被排除在外，以隔离动态域映射的数值效应；采用与颗粒尺寸成正比的线性生长定律。采用一阶上风向差分，FlexPBE相较于其固定域对应方案将误差指标降低了约95%，并且在相同设置下优于固定网格二阶上风向差分和HR-Van Leer方案。在计算成本相当的情况下，FlexPBE使用1000个网格所获得的误差低于固定域仿真中使用14,000至15,000个网格的结果。最后，FlexPBE被集成到基于OpenFOAM的MP-PIC-PBE求解器中。耦合仿真保持数值稳定，并且在相同的30网格条件下，FlexPBE 30产生的PSD形状和矩特征比固定网格PBE 30更接近于更高分辨率固定网格结果。

## 引言

颗粒在多种科学与工程领域中发挥着至关重要的作用，其重要性随时间推移而日益凸显(Rhodes, 2008)。颗粒尺寸分布(PSD)被视为许多过程中的关键状态变量，因为它决定了反应和输运特性、产品质量、压降和设备性能(Shekunov等人, 2007)。群体平衡方程(PBE)描述了含有分散相系统中关于内部坐标(如颗粒尺寸、体积)和空间坐标的数密度分布的时间演化，已成为预测各种分散相过程中PSD的成熟建模工具(结晶、团聚、研磨、气溶胶、气泡和液滴分散体系等)。PBE的概念最初由Hulburt和Katz(1964)提出，他们建立了一个统计框架来描述由于成核、生长、团聚和破碎而变化的颗粒尺寸分布(PSD)。其最初的公式主要考虑颗粒尺寸分布作为随时间演化的变量；然而，后来建立了一种考虑三维物理空间和单一颗粒尺寸坐标的一般多维PBE形式(Ramkrishna, 2000)。这项努力使得PBE能够通过计算流体力学(CFD)的耦合来考虑空间异质性，Marchisio和Fox(2005)提出了最广泛使用的CFD–PBE实现方法。

在PBE中，颗粒尺寸的变化主要归因于四种机制，各机制的解释和发展进程如下：

**• 成核**

成核是指新颗粒的形成过程，通常仅作为源项添加在某一特定尺寸以下区域。20世纪60至70年代早期的研究引入了简单的成核常数(Randolph和Larson, 1971)。后来，通过与经典成核理论(CNT)耦合，发展成为结晶或凝华过程中的成核速率模型(Kashchiev, 2000; Mullin, 2001)。最近，人们报道了"两步成核"机制，其中先形成亚稳态(高密度)团簇，然后在其内部生成晶体晶核(Vekilov, 2004)。该机制更真实地描述了成核过程，已被讨论为可能补充传统CNT无法充分解释的方面，例如非晶态液相前驱体的存在或成核发生得比预期更快(Gebauer和C?lfen, 2011)。

**• 生长**

生长是指从周围相(溶液、气体、连续相)中供应的物质(如溶质、可凝结蒸气或单体)被输送并附着到颗粒表面，导致颗粒长大的过程。在PBE中，它表现为内部坐标空间中的对流项。早期，生长用与尺寸无关的零阶生长速率公式来表示(McCabe, 1929)，但逐渐演变为尺寸依赖型生长速率模型(Abbegg等人, 1968)。后来，Randolph和Larson(1971)在数学上建立了结晶理论并系统地定义了生长项。

**• 团聚**

团聚涉及两个或更多颗粒结合形成较大颗粒，同时包含产生和消亡两种效应。后来，基于团聚理论引入了碰撞核(Smoluchowski, 1918)，并随后被纳入PBE框架。根据核的形式，出现了各种模型，如布朗团聚、湍流团聚和电性团聚(Saffman和Turner, 1956; Thomas等人, 1999)。此外，一种确定性方法基于Smoluchowski方程以群体平衡方程的形式解释碰撞核，有利于预测平均晶体尺寸分布(CSD)并与CFD耦合(Marchisio等人, 2003)。相比之下，基于Monte Carlo的随机方法能够在没有数值扩散的情况下自然地捕捉复杂的颗粒特性和变异性，但方差控制和计算成本被认为是关键挑战(Xu等人, 2014)。

**• 破碎**

破碎描述了颗粒在外部力的作用下碎裂成较小颗粒的过程，与团聚机制一样，同时包含产生和消亡两种效应。它在研磨和冲击现象的背景下得到了发展，Epstein(1947)在数学上建立了破碎机制。后来，Tavares和King(1998)以及Broadbent和Callcott(1956)等工作中进一步发展了详细函数，并引入了破碎函数和母粒分布函数等概念。由于控制机制多样且强烈依赖于颗粒的材料特性和内部缺陷，并且有许多因素使得从第一性原理中完整地闭合核变得困难，因此破碎不像生长或团聚那样，往往依赖于基于实验的经验关联(Austin等人, 1984)。

近年来，许多研究致力于将CFD与PBE耦合。这些工作源于预测系统内空间异质性的需要，这会影响最终的PSD。由于CFD–PBE耦合因分布函数的高维性增加了计算成本，因此需要高效的数值方法来用低阶矩概括分布。从计算方法学角度，矩方法(MOM)已被用于解决这一耦合问题。然而，由于MOM仅求解有限数量的低阶矩，当非线性积分项主导分布时(如团聚/破碎中需要高阶矩的情况)，可能会出现闭合问题。为了解决这一问题，提出了QMOM，其中分布由少量求积节点和权重近似，并通过求积闭合矩积分(McGraw, 1997)。此外，DQMOM通过将(wi(权重), ξi(求积节点))本身作为未知数并直接以输运方程形式求解，从而改善了在CFD–PBE耦合中的适用性和计算效率(Fox, 2003; Marchisio和Fox, 2005; Fan等人, 2004)。然而，由于基于MOM的方法将分布压缩为少量自由度，与直接离散化尺寸轴并计算n(x,t)的方法相比(Dürr和Bück, 2020)，其再现多峰分布、长尾或尖锋等形状的能力可能有限。此外，由于满足同一组矩的分布不唯一，难以唯一重建PSD形状(Yuan等人, 2012)。另外，如果数值误差破坏了矩序列的实性，可能出现负权重或非物理节点，需要额外的稳定和校正(Vikas等人, 2011)。因此，为了可靠地预测不仅矩而且颗粒分布的形状信息，需要一种基于直接尺寸轴离散化的CFD–PBE数值框架来补充或替代基于矩的约简技术。

在对混合流敏感的结晶过程中(如反溶剂结晶)，局部微观混合决定了瞬时过饱和度分布，这通过成核和生长速率的空间不均匀性直接影响最终的CSD。为了定量预测这种关联，提出了CFD–PDF–PBE框架：使用基于RANS的CFD分析宏观混合，通过多环境假设PDF模型引入亚网格微观混合，并耦合空间变化的PBE来计算CSD(Woo等人, 2006; Pirkle等人, 2015)。然而，虽然该方法主要在Euler框架下描述了混合-反应-结晶的空间异质性，但最近有尝试从Lagrangian视角追踪分散相以加强流场与颗粒分布之间的多尺度耦合。Kim提出了MP-PIC-PBE框架，在其中对每个Lagrangian包迹求解PBE，实现了流场与颗粒分布之间的多尺度耦合(Kim等人, 2020, 2021, 2024)。

要数值求解完整的PBE，需要对颗粒尺寸分布进行离散化。在此设置中，分布域长度和离散节点数作为数值参数，在计算时间和精度之间产生权衡(Kumar和Ramkrishna, 1996a)。为解决此问题，可应用高阶离散化方案、通量限制器或自适应网格加密(AMR)，后者局部调整网格分辨率(Berger和Colella, 1989)。然而，这些技术可能难以实现且可能计算昂贵(Gunawan等人, 2004)，或者可能面临通量耗散过大等限制(LeVeque, 2002)。同时，在以团聚和破碎为主导的PBE系统中，移动枢轴技术等自适应离散化方法被广泛采用。移动枢轴技术允许每个离散区间内的代表性颗粒尺寸(枢轴)根据该区间内数密度不均匀性的演化动态移动，从而精确保持任意选择的两个矩(如颗粒数和颗粒质量)，同时减轻团聚和破碎导致的大颗粒群体过预测(Kumar和Ramkrishna, 1996b)。然而，由于枢轴运动规则与团聚和破碎等特定机制的控制方程紧密关联推导，难以根据给定情况下存在哪些机制自由添加或移除枢轴再分配程序。相比之下，本研究提出的FlexPBE框架被设计为一个包装器，独立于离散化方案的内部结构和单个机制的控制方程运行，允许相同的边界更新程序在不同机制组合(如生长和成核)中一致地应用。与移动枢轴和其他紧密耦合于特定机制(如团聚和破碎)的自适应技术不同，以生长为主导的系统需要一种方法来解决上述问题(实现复杂性、计算成本、数值扩散)而无需这种机制特定的耦合。

为此，本研究提出了一个建立在动态域映射(Dynamic Domain Mapping)既有概念之上的自适应计算框架。所提出的方法通过压缩计算域并将网格仅放置在粒子存在的区域，在有限的计算资源下提高效率。计算域的边界位置使用颗粒尺寸分布的Z分数计算。虽然计算域被归一化至[0,1]，但离散化方案考虑了边界运动的变化，将PBE仿真所需的微分项转换为反映动态建模变化的微分形式。我们称转换后的方程为FlexPBE(柔性域PBE)。

对于本评估，成核和晶体生长被视为控制本研究所考虑结晶过程中PSD演化的基本机制，FlexPBE使用采用单一线性正比于颗粒尺寸的定律的批次生长和批次成核-生长基准进行验证。为了清晰评估边界更新和域映射程序的效果，引入额外产生-消亡和颗粒再分配过程的团聚和破碎被排除在本验证之外；其纳入留待未来工作。最后，FlexPBE被集成到用于反溶剂结晶的MP-PIC-PBE求解器中，以评估其在有限内部坐标分辨率下的数值稳定性和PSD预测的精度-成本比。

## 各节片段

**群体平衡方程**

群体平衡方程是描述颗粒尺寸分布(PSD)如何随空间和时间变化的数学模型。

公式(1)给出了PBE的一般形式（Ramkrishna, 2000）：

$$\frac{\partial n(x,r,t)}{\partial t} + \nabla \cdot \left(u(r,t) n(x,r,t)\right) + \frac{\partial}{\partial v}\left(G(x,r,t) n(x,r,t)\right) = B_{nuc}(t) \delta(x - x_0) + B(x) - D(x)$$

其中 $n(x,r,t)$ 为具有尺寸坐标 $x$、空间位置 $r$ 和时间 $t$ 的颗粒数密度函数；$u(r,t)$ 为流体速度矢量场；$G(x,r,t)$ 为颗粒生长速率；$B_{nuc}(t)$ 为成核速率。

**案例研究方案**

FlexPBE通过与其他现有数值方法的比较进行了验证。所有算法最初均在MATLAB R2024b中实现并进行数值评估。本研究通过将其应用于两个排除聚集和破碎的案例来评估所提模型的有效性：案例1为间歇生长基准测试，案例2为间歇成核-生长基准测试。案例1考虑一个间歇反应器，其中生长速率与颗粒尺寸成线性正比关系，以便分析……

**结论**

本研究提出了一种新开发的计算框架FlexPBE，通过动态适应计算域以匹配不断演化的分布，从而改进颗粒尺寸分布的数值预测。所提方法将分辨率重新分配至颗粒实际存在的区域，从而克服了固定域离散化中常见的分辨率损失和数值扩散问题。其性能通过具有解析解的基准测试进行了验证。

**CRediT作者贡献声明**

Junho Chung：写作—审阅与编辑，写作—初稿撰写，可视化，验证，软件，方法论，调研，资金获取，数据管理，概念构思。Hyeonseung Jo：写作—审阅与编辑，验证，软件，形式分析，数据管理。Shin Hyuk Kim：写作—审阅与编辑，验证，指导，项目管理，方法论，调研，资金获取，概念构思。

Junho Chung | Hyeonseung Jo | Shin Hyuk Kim