結晶是一種重要的分離和純化單元操作，已廣泛地應用于化工、醫藥和精細化工等行業。特別是在制藥行業中，大部分的產品以固體的形式進行生產和銷售，對這類產品質量的一致性要求也越來越高。在藥物晶體產品的流程開發上，生產者面臨著來自監管和經濟的雙重挑戰。從監管的角度看，晶體產品的開發應基于對系統的過程理解和控制上，也確保晶體產品的質量符合 “質量源于設計 (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