化学工业与工程杂志投稿格式参考范文：计算流体力学在结晶过程模拟仿真中的应用

　　结晶是一种重要的分离和纯化单元操作，已广泛地应用于化工、医药和精细化工等行业。特别是在制药行业中，大部分的产品以固体的形式进行生产和销售，对这类产品质量的一致性要求也越来越高。在药物晶体产品的流程开发上，生产者面临着来自监管和经济的双重挑战。从监管的角度看，晶体产品的开发应基于对系统的过程理解和控制上，也确保晶体产品的质量符合 “质量源于设计 (QbD)” 的理念，因此，必须获得与过程相关的充分信息，以确保质量和设计之间的相关性;从经济的角度看，结晶仍是一门 “半科学与半艺术” 的学科，结晶过程缺乏相应的成本效益控制和可靠的数学建模工具，可能需要消耗大量的时间和资源。因此，有必要建立结晶过程的模型，使其应用到结晶的设计与过程控制中，并据此制定一个更加高效而稳健的流程。

　　结晶过程主要包括 2 个基本步骤，即晶体成核过程和晶体生长过程。溶质从过饱和溶液析出到晶体的进一步生长，显著影响晶体的纯度、形态和粒度分布 (CSD) 等性质。而结晶系统的混合性能、质量和热量传递效果会直接影响结晶动力学过程 (成核与生长)。因此，流体力学研究对于结晶过程以及晶体产品的精准控制至关重要。近些年来，随着计算机计算能力的提高，计算流体力学 (CFD) 依靠其独特优势，被广泛应用于结晶领域。此外，基于描述颗粒数量密度函数的粒数衡算方程 (又称群体平衡方程 PBE)，结合结晶动力学模型能够预测晶体的尺寸分布。同时，为了解决非理想混合问题，可以将 PBE 与质量、动量和能量输运方程相结合，共同模拟晶体的尺寸与形态的演化过程，建立 CFD-PBE 耦合模型，为深入理解结晶过程提供新的思路和方法。

　　1 粒数衡算方程及数值解法

　　粒数衡算方程是一个用来描述颗粒数量密度函数的偏微分方程，在结晶过程中通常引入数密度函数 (Number density function，NDF) 的概念来表示颗粒群。NDF 表示颗粒在任何时间和空间上所关注的属性上的数量分布 (n (t,x,ξ))，这些属性可以是颗粒的尺寸与组成，也可以是温度 T、质量 m 等，这些属性称为颗粒的内坐标 ξ。值得一提的是，如果内坐标为颗粒的速度 v 时，则粒数衡算方程 (PBE) 转变为普适性粒数衡算方程 (GPBE)。

　　颗粒的聚集与破碎受流体力学环境的影响很大，如在搅拌釜中结晶时，在叶轮附近的高剪切力区域，晶体颗粒会经历更强的破碎过程，而在其他低剪切力的区域，晶体颗粒可能经历更多的聚集过程。所以聚集核和破碎核函数高度依赖于外坐标，此时破碎和聚集可能对方程 (1) 起到主要的控制作用。

　　一般而言，PBE 足以完成结晶过程的模拟，但结合方程 (1) 和方程 (2) 来看，PBE 是一个典型的双曲型积分 - 偏微分方程，只有在非常严格的假定下，PBE 具有解析解。所以通常情况下，需要采用数值解法对 PBE 进行求解。目前开发的方法包括矩量法 (MOM)、分组法 (CM)、加权残差法 (WRM) 和一些随机方法等。下面重点综述了矩量法和分组法在求解粒数衡算方程中的应用。

　　1.1 矩量法 MOM：矩量法最开始由 Randolph 及 Hulburt 等提出。该方法并不直接跟踪 NDF，而是通过在内坐标上积分的矩量来跟踪，将 PBE 转换为一系列的 NDF 矩量输运方程。利用这些矩量，可以计算得到颗粒分布的整体属性，如颗粒数量、表面积、体积、平均尺寸和变异系数等信息。

　　虽然 QMOM 方法多数情况下是为气溶胶的研究而开发，但可以扩展到结晶过程的研究中，并在耦合 CFD 与 PBE 方面显示出了显著优势。如 Marchisio 等使用 QMOM 对晶体破碎进行建模，得到了不同聚集破碎核函数、子分布和初始条件下的计算结果，并与 Vanni 等利用 CM 方法得到的结果进行了比较。结果表明，QMOM 能够以较小的误差预测多个矩量。在上述研究基础上，Jin 等对 QMOM 算法进行改进，该算法提高了 QMOM 求解的鲁棒性，并且在高阶矩量的求解上表现优异。

　　然而，QMOM 存在 2 个比较严重的问题，一是对于内坐标和颗粒速度存在强耦合的多相分散体系，仅仅通过求解矩量方程难以得到准确的数值解;二是对于多维 PBE，使用 QMOM 方法来求解十分复杂和低效。因此，在这个基础上，又发展出了许多新的方法，如直接正交矩量法 (DQMOM) 和扩展正交矩量法 (EQMOM)。Selma 等在 OpenFOAM 中植入了 CM 和 DQMOM，对鼓泡塔和搅拌釜中的多相流动进行建模。结果表明，在平均轴向速度、气含率和气泡直径方面，数值计算与实验结果均吻合较好，且 DQMOM 的计算成本明显低于 CM。Askari 等在 OpenFOAM 中采用 EQMOM 对径向和轴向 2 种搅拌釜的气液两相流动进行了数值模拟，数值结果在轴向和径向液体速度分布和溶解氧的浓度变化方面与实验结果相吻合。此外，条件正交矩量法 (CQMOM) 也是一种基于条件密度函数 (CDF) 的数值解法，适用于多维 PBE 的算法，并可以处理具有多个内部坐标的复杂体系。

　　综上，目前已经开发了大量的矩量法来解决一维或多维 PBE，证明其是一种有效地求解 PBE 的方法，特别是在与 CFD 相结合的方面，矩量法可以在计算精度和计算资源之间取得了很好的平衡。但由于矩量法不能直接预测 NDF，其自身的局限性与面临的数值挑战依然存在，未来仍需进一步地研究更加健壮而高效的算法。

　　1.2 分组法 CM：分组法也称为离散法，其基本思想是将颗粒的内坐标划分为若干个相邻的子区间，每个小区间的颗粒群用一个代表节点表示，在每个子区间内分别进行积分，从而将 PBE 转化为一系列离散的方程组。

　　为此，Kumar 等提出了固定节点法 (Fixed pivot method)，通过将聚集或破碎产生的新颗粒按一定的比例将其分别分配到邻近的节点处，从而保证颗粒的总个数和总质量守恒，并且可以用于任意的网格类型，是目前使用最广泛的方法。另外，Alopaeus 等提出，将新颗粒分配到 2 个以上的节点处时，可以进一步提高守恒特性。然而，固定节点法虽然对低阶矩的数值计算结果比较精确，但是对大颗粒数密度和高阶矩的预测经常偏高，特别是在存在聚集的情况下。因此，Kumar 等随后提出了移动节点法 (Moving pivot) 和单元平均法 (Cell average method) 来解决这一问题，该方法在很大程度上缓解了固定节点法预测偏高的问题。

　　在数值求解 PBE 时，需要考虑这些尖锐的界面和不连续性对数值解的影响，以确保数值解的准确性和稳定性。所以，必须谨慎处理方程 (10) 的对流生长项，以避免引入人工黏性。一方面可以使用限制器，如 MUSCL 限制器和 WENO 限制器，这些限制器可以在界面处自适应地调整数值解，以避免数值振荡和人工黏性的出现，从而保持解的高精度。Jiang 等提出了一种与五阶 Runge-Kutta 方案结合的 WENO 格式，它允许以非常小的时间步长进行高效率求解，并在结晶过程中得到了一些应用。另一方面可以使用高分辨率数值方法，如高分辨率有限体积法和高分辨率有限差分法，这些方法可以有效地处理尖锐的界面或不连续性，同时保持解的高精度。Kurganov 等提出了一种高分辨率有限体积格式 (K-T 格式)，由于该方法满足 minmod 重构的全标量递减特性，因此不会产生非物理振荡。在此基础上，Cheng 等将其扩展到通用化的形式，如通用 K-T 高精度有限体积中心差分格式，以便用于均匀或几何离散的内网格中。相对于迎风差分格式，特别是在颗粒比较大时，K-T 高精度有限体积中心差分法计算精度的提升较为显著。

　　综上，相比于矩量法，分组法的优势在于可以直接预测连续的 NDF，缺点在于为了捕捉完整的 NDF 分布，可能需要将内坐标离散成大量的组，计算成本高。但随着计算机技术的飞速发展，在可预见的未来内算力将得到快速提升，分组法也将获得更好的应用前景。

　　2 CFD-PBE 在结晶过程数值模拟中应用

　　对结晶过程的 CFD 模拟主要集中在溶析结晶和反应结晶过程中，CFD 模型包括不同尺度的混合，即微观混合、介观混合和宏观混合。以搅拌釜为例，展示了湍流状态下不同尺度的混合效应 。由于微观混合能够在分子水平上影响着溶质与溶剂分子的扩散程度和体系的过饱和度，因此其直接影响晶体的成核和生长过程，并决定产品的纯度、粒度晶体形貌等质量指标。

　　CFD 中的宏观混合模型主要由湍流输运方程组成，对微观混合有重要的间接影响，即微观混合和宏观混合是相互依赖。此外，在 CFD 框架中介观混合通常被归入宏观混合，而微观混合模型成为宏观混合的子网格建模组件。因此，当采用 CFD 方法来描述宏观输送时，必须包括相互兼容的微观混合建模。在反应和溶析结晶过程研究中，常用的微混合模型包括涡流卷吸模型 (E 模型) 和基于联合组成概率密度函数模型 (PDF) 模型。

　　2.1 溶析结晶

　　溶析结晶由于不产生较大的温度变化而有利于热敏性药物的结晶，尤其广泛应用于制药行业中以生产粒径小而分布窄的晶体。由于此方法需要将溶析剂与溶液快速而充分地混合，因此基于 CFD 可以获得对结晶过程中混合效应的深入理解。目前，科研人员已经对混合在溶析结晶过程中的影响进行了细致的研究工作，主要集中在结晶器基础设计和模拟方法优化改进 2 个方面。

　　结晶器设计方面，近年来先后开发了 Y 型结晶器、同轴射流结晶器以及多种横向、径向撞击射流结晶器等。Choi 等采用 CFD-PBE 的方法，利用矩量法求解了 PBE 并模拟了环四亚甲基四硝胺在小型的 Y 型结晶器中的稳态溶析结晶。在不考虑微混合的条件下研究了速度比、过饱和度、平均晶体尺寸和变异系数的空间分布，发现即使在小型的结晶器中，混合效应依旧显著。Woo 等采用 CFD-PBE-PDF 方法，研究了洛伐他汀在横向撞击射流结晶器的溶析结晶。结果表明，微混合对 CSD 的影响取决于微混合特征时间和结晶诱导时间的相对大小。Pirkle 等在 Woo 的基础上，同样采用 CFD-PBE-PDF 的方法考虑了混合热与结晶热的影响，模拟了洛伐他汀在同轴射流结晶器中的溶析结晶过程。结果表明，较高的流量带来较短的停留时间和高雷诺数，从而得到较窄的 CSD 分布和较小的颗粒平均尺寸。Rosa 等提出了一种多入口径向结晶器的新设计，与只有 1 个入口的径向结晶器相比，多个入口的设计提供了更好的混合。Chen 等在此基础上进行改进，将多入口径向结晶器扩展到多位置入口的径向结晶器的设计，结果表明，多位置进料点的设计有助于在保持良好混合的同时，降低局部过饱和度从而获得更窄的 CSD 和更大的平均晶体尺寸。上述研究表明，采用 CFD-PBE 耦合方法是溶析结晶器开发的有力工具，为获得粒度均匀的晶体，径向撞击射流结晶器是撞击流结晶器发展的新方向。

　　上述工作虽然采用了 CFD-PBE 耦合方法对溶析结晶过程进行了数值模拟，然而模拟过程均是在拟均相假设下进行的。为了进一步提升模拟准确性，Cheng 等在 CFD-PBE-PDF 方法基础上，进一步考虑了颗粒和溶液间的相互作用，采用 Mixture - 欧拉多相流模型模拟了洛伐他汀在横向射流结晶器中的溶析结晶过程。结果表明，随颗粒尺寸增大，颗粒与溶液间的相互作用不可忽视，拟均相假设的偏离逐渐增大。在此基础上，Farias 等提出了耦合 PBE 的欧拉 - 欧拉模型来模拟同轴结晶器中的冷却 - 溶析结晶过程。结果表明，溶剂与溶析剂的质量分数、成核生长速率的分布可能是不对称的，且入口处的进料速度对最终粒度分布产生显著影响。总之，在工业规模的计算中，拟均相假设因其计算量小而被大规模使用，但这是以牺牲计算精度为前提的。随着计算机能力的提升，在计算中充分考虑颗粒间相互作用，无论是采用基于欧拉框架的 CFD-PBE-PDF 方法还是基于拉格朗日框架下的 MP-PIC-PE 方法，都将是未来溶析结晶建模的主流方法。

　　2.2 反应结晶

　　早在 1958 年，Danckwerts 便讨论了混合对反应结晶的影响。当特征反应时间小于特征混合时间或两者处在相同数量级时，混合效应对反应结晶中悬浮液的性质、颗粒产品的产率、纯度等许多方面产生显著影响。如 Wei 等采用 CFD-PBE 模型模拟了半间歇搅拌釜中 BaSO4 的反应结晶过程，得到的晶体平均尺寸随转速变化与实验结果相差接近 4 倍。Jaworski 等采用 CFD-PBE 方法对双进料连续搅拌釜中 BaSO4 反应结晶过程进行了建模，但与实验结果存在较大偏差。以上模拟结果的偏差主要是因为未考虑微混合效应且采用 SMOM 方法求解 PBE，忽略了结晶过程中的聚集和破碎过程。

　　针对这一问题，Cheng 等对同样双进料搅拌釜中的 BaSO4 反应结晶进行 CFD-PBE 建模，利用 QMOM 方法对 PBE 求解，考虑了聚集和破碎的过程。结果表明，在较高进料浓度下，模拟值与实验值偏差较大。随后 Cheng 等利用离散法求解 PBE，得到了完整的 CSD，但与实验结果的 CSD 相比峰值略大，且分布略宽。以上研究工作中模拟结果与实验结果的偏差，主要是由于双进料口距离较远和预混合进料的原因，并未考虑微混合效应对实验的影响。在此基础上，Marchisio 等采用 CFD-PBE-PDF 方法对 BaSO4 反应结晶进行建模，结果表明，在不考虑颗粒聚集时，高浓度反应物下的模拟结果较差，而考虑颗粒间的聚集效应后，得到了更好的结果。这也证明 BaSO4 反应结晶中，聚集破碎现象明显，在高浓度体系下甚至主导了产品质量。为了更好地描述聚集破碎问题对结晶过程的影响，Bal 等开发了一个包含混合、反应成核、生长、聚集和 Ostwald 熟化的耦合 CFD-PBE 模型，能够在广泛的操作条件下预测二氧化硅纳米颗粒的大小。

　　从上述研究工作结果中可以看出，相比于溶析结晶，反应结晶过程中的聚集和破碎现象更加常见，因此除了要关注微混合效应，对晶体聚集和破碎的建模也更加重要。因此，对 CFD-PBE 耦合模型的完善和计算方法的改进仍需加以深入探讨和研究。

　　3 CFD-PBE 耦合模拟在结晶器放大设计中的应用

　　在实验室规模的结晶器中，通常能够达到良好的混合进而有利于质量和热量的传递，且整个溶液中的过饱和度几乎保持一致。然而，随着结晶器从实验室规模增大到工业规模，流动、混合、传质和传热的一致性放大成为一个巨大的挑战。为了实现这一目标，这就需要对结晶器的工业放大实现模型化和数字化。

　　结晶过程的准确放大需要充分了解混合及其对结晶过程中关键性能的影响。对于冷却结晶过程，其主要受传热影响，微观和宏观混合影响通常较小。但与之相比，溶析结晶和反应结晶受微观混合的影响最大，结晶器中的流体力学决定了流动、传质、传热以及晶体的分散特性。同时，在溶析和反应结晶过程中涉及众多影响过饱和度分布和产品质量的关键工艺参数，如溶析剂添加方式、反应物混合方式、混合速度和结晶器型式等。因此，在溶析结晶和反应结晶放大过程需重点对流体力学性质进行系统研究。此外，结晶器的设计参数和操作参数决定了流体力学特性，因此，在相应结晶设备开发过程中同样需要重点关注设备特点对流体力学的影响，包括 (1) 结晶器几何形状 (尺寸、纵横比等);(2) 结晶器内部配置 (叶轮和挡板类型、尺寸、位置、数量等);(3) 溶液和晶体的物理化学属性;(4) 操作参数 (转速，填充率，进料速率、位置，流速等) 。前文主要综述了 CFD-PBE 耦合技术应用于辅助实验研究和实验过程设计，而 CFD 在辅助结晶器放大过程中同样发挥了十分重要的作用，有助于加深对结晶过程理解、缩短开发周期、降低开发成本、优化结晶器设计等。

　　Kramer 等对工业结晶器进行了系统研究，提出了工业结晶器的设计和放大策略：对于溶析结晶，溶析剂的添加速率起主要控制作用，并且需要注意由于溶析剂与溶液的混合产生的热效应。对于反应结晶过程，反应物的添加速率和反应热是影响放大的关键参数。以结晶中最常见的 DTB 结晶器为例，Zauner 等开发了一种半间歇和连续操作下，双进料反应结晶 DTB 结晶器的 CFD-PBE 模型，在实验室规模下测定了结晶动力学，并将模拟结果与不同规模的结晶器中获得的实验数据进行对比。结果表明，当结晶由微观或介观尺度的混合控制时，恒定搅拌速率、恒定叶尖速度和恒定的比功率等常规的放大策略不可行。针对这一问题，Oh 等通过实验和 CFD 模拟研究了 500 和 300L 的 DTB 溶析结晶器的放大，证明叶尖速度、Froude 数和比功率是主要放大标准。并指出双引流管的 DTB 比单引流管 DTB 轴向流动可获得充分发展，因而可以获得更均匀的混合效果，表现出更好的性能。

　　但从上述 2 个案例中看出，恒定叶尖速度、恒定比功率输入，在结晶器放大过程中也并不是都适用的，不同的操作条件，不同的结晶方式都对结晶器放大产生着影响。这是因为目前对于结晶器放大设计过程中，大多是基于实验 - 整体放大 - CFD 建模的策略，其中一个重要的因素是保持结晶器的几何相似性。理想情况下，混合、悬浮、流体剪切、传热等指标在放大时都应该保持恒定。但实际情况是，保持一项指标恒定的同时总会使其他指标发生变化。因此，这些指标往往是相互制衡的，很难实现所有指标的恒定。针对几何相似放大的问题，需要明确并非要对所有指标进行复现，许多结晶过程并不是对这些指标都很敏感，这可以大大简化放大工作。因此，未来的研究工作中，可以在放大之前通过中等规模的实验与 CFD 模拟，明确敏感变量，进而为放大工作提供理论基础。此外，还可以考虑放弃几何相似的制约，便于使不同结构、不同尺寸的结晶器中有更多的参数保持一致。

　　4 总结与展望

　　综述了用于结晶过程的数值模拟的 PBE 与 CFD 耦合求解的几种数值方法，重点介绍了矩量法与分组法。在此基础上，综述了 CFD-PBE 方法在溶析结晶与反应结晶过程模拟以及结晶器设备放大领域的研究进展。整体来看，虽然直接数值模拟可以获得不同结晶方式条件下全尺度上的混合信息，但在目前的计算能力下，仍然需要附加额外的微混合模型来建模。此外，尽管我们可以通过几何相似性来复制一些参数，但在不同规模下，流体力学特性、传质特性、混合过程等因素给设备放大开发带来诸多挑战。同时，晶体的形成过程也受到初始溶液浓度、温度、pH 值、添加剂等因素的影响，这使得从实验室到工业规模的放大更加困难。

　　因此，虽然 CFD 在辅助结晶过程开发和结晶器设计上是一个强有力的工具，然而目前研究仍存在一定不足，未来的研究中可重点关注以下几方面：首先，在利用矩量法求解 PBE 时，仍然需要建立更加稳健高效的算法来解决从低阶矩量重构 CSD 的问题;而在分组法上，虽然计算能力的提升促进了分组法的发展，但求解聚集破碎和生长的数值方法的提升依然是一个挑战。其次，将一维 PBE 扩展到多维 PBE 也是一项重要的工作，利用计算流体力学技术模拟晶体的尺寸和晶习是未来结晶过程模拟的一个重要研究方向。最后，在结晶器放大设计中，进行多尺度的 CFD - 微混合 - PBE 建模可以更深入地揭示结晶过程中流体动力学和结晶动力学之间的相互作用机理，帮助我们建立更准确的放大策略。

郜 源;花玉香;张 静,天津渤海职业技术学院;天津海关动植物与食品检测中心,202406