结晶是分离和提纯的重要操作，广泛应用于制药、精细化工和食品等行业。在结晶过程中，晶体尺寸分布(Crystal size distribution, CSD)或者弦长分布(Chord length distribution, CLD)通常被用作评价晶体产品质量的重要指标^[1]。

与传统间歇结晶操作相比，连续结晶过程可以节省9%~40%以上的生产成本^{[2, 3]}，并且稳态操作中较低的物料滞留率和对最终药品状态的优越控制使得连续生产在结晶过程中有很大应用价值^[4]。因此，如何使结晶过程由批次方式向连续方式转变，以更低的成本、更高的效率和更好的控制效果来提高晶体质量，已成为工业界和学术界的重要研究课题^[5]。目前最常使用的连续结晶器主要有多级连续冷却结晶器和管状结晶器。本研究使用一种管状结晶器，即连续振荡挡板式结晶器(Continuous oscillatory baffled crystallizer, COBC)，这是一种振荡流反应器(Oscillatory flow reactor, OFR)。COBC可在广泛的运行条件下以近平推流运行，使其内部环境中剪切速率、温度梯度和浓度分布更为均匀^[6]。

L-谷氨酸(L-glutamic acid, LGA)存在于许多食品中，并在许多神经退行性疾病治疗中发挥重要作用^[7]。LGA具有2种已知的具有不同特征形态的多晶型。亚稳态的α晶型呈棱柱状，而热力学稳定的β晶型呈板状或针状。LGA在连续振荡流中进行结晶的研究很少见^[8]。

在冷却结晶的生产过程中，冷却速率、过饱和度及混合程度等过程条件会对最终产品产生很大的影响^[9]。随着先进的在线传感器在结晶过程中的应用，基于过程分析技术的结晶过程控制与优化成为可能^[10]。为了对结晶过程进行模拟并优化过程操作条件，通常使用结晶热力学模型、溶质质量平衡方程、粒数衡算方程以及成核与生长方程^{[11, 12]}。然而，现有文献中的大部分机理模型都是基于理想化假设得出的简化模型^[13]，因此优化效果受到影响。相比之下，基于实际结晶过程检测数据的建模方法近年来受到了越来越多的关注^{[14, 15]}。

本研究提出1种基于实验设计的晶体产品尺寸分布预测模型和过程操作优化方法。将连续管式结晶器的振荡频率、2个区段的降温速率作为操作变量，给出1个响应面实验设计方案。通过检测不同实验条件下的晶体产品弦长分布，构建1个基于双层基函数的晶体产品尺寸分布预测模型。然后提出1个关于目标晶体尺寸分布和产品收率的目标函数，基于上述预测模型建立优化过程操作条件的方法。通过计算机仿真和实验，验证该方法的有效性和优点。

1 连续振荡挡板式结晶器

COBC由一系列挡板玻璃直管与弯管组成，如图 1所示。

在COBC起始管段连接1个振荡源。主进料被泵入系统以提供净流量。在各管段结晶区域，外接温度调节装置进行加热和冷却操作。需要说明，活塞推进器必须保证一定的振荡频率和振幅以防止晶体沉降或堵塞。连续振荡流可以用3个无量纲数来描述：振荡雷诺数Re₀、斯特劳哈尔数St和净流雷诺数Re_n^[16]。这些参数由式(1)~式(4)计算。

$ R e_0=\frac{2 \pi f \chi_0 \rho d}{\mu} $

(1)

$ S t=\frac{d}{4 \pi \chi_0} $

(2)

$ R e_{\mathrm{n}}=\frac{\rho u d}{\mu} $

(3)

$ \psi=\frac{R e_0}{R e_{\mathrm{n}}} $

(4)

式(1)~式(4)中：d是管径，m；χ₀是中心到峰值的振幅，m；f是振荡频率，Hz；ρ是密度，kg·m^-3；μ是流体黏度，kg·m^-1·s^-1；u是平均流速，m·s^-1。为达到混合效果，应满足Re₀≥Re_n。速度比ψ可由式(4)计算。关于振荡流的设计^[16-19]，需要满足Re₀≥100，St≤0.5，Re_n≥50，2≤ψ≤10。因此为实现溶液良好混合，并使流体形式接近平推流^[12]，可行的操作条件范围如表 1所示。

目前尚未有文献报道如何对COBC过程操作条件进行优化。为此，本研究提出一种基于实验设计的产品尺寸预测模型，由此给出过程优化方法。

2 实验设计与数据驱动建模 2.1 实验设计

响应面设计(Box-Behnken design, BBD)是一种基于不完全因子数的实验设计(Design of experiment, DOE)方法。一个3因素的响应面实验设计方案是一个由中心点和边界的中点组成的立方体，如图 2和表 2所示。

说明一下，BBD的实施效率高于中心组合实验设计方法以及3水平全因子设计方法^[20]。

2.2 双层基函数预测建模

基于近期提出的一种数据驱动建模方法^[21]，这里给出1种基于过程操作条件预测晶体产品尺寸分布的建模方法。假设所有的输入条件记为u=[u₁, …, u_i, …, u_I]^T，这里I代表实验批次的数目。γ(u_i, y_k) 代表输入条件u_i对应不同颗粒尺寸y的密度函数。假设各批次操作条件u_i下采集晶体颗粒尺寸区间都是相同的，记为y=[y₁, …, y_k, …, y_K]^T，这里K是各批次采集晶体尺寸区间的个数。CSD的密度函数记为：

$ \begin{array}{c} \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = \left[ {{\mathit{\boldsymbol{\gamma }}_1}, \cdots , {\mathit{\boldsymbol{\gamma }}_k}, \cdots , {\mathit{\boldsymbol{\gamma }}_K}} \right] = \\ \left[ {\begin{array}{*{20}{c}} {{\gamma _{1, 1}}}& \cdots &{{\gamma _{1, k}}}& \cdots &{{\gamma _{1, K}}}\\ \vdots &{}& \vdots &{}& \vdots \\ {{\gamma _{i, 1}}}& \cdots &{{\gamma _{i, k}}}& \cdots &{{\gamma _{i, K}}}\\ \vdots &{}& \vdots &{}& \vdots \\ {{\gamma _{B, 1}}}& \cdots &{{\gamma _{B, k}}}& \cdots &{{\gamma _{B, K}}} \end{array}} \right] \end{array} $

(5)

这里γ_{i, k}=γ(u_i, y_k)。定义过程操作条件(即模型输入变量) u的多项式基函数矩阵：

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \left[ {\begin{array}{*{20}{c}} {{\varphi _1}\left( {{\mathit{\boldsymbol{u}}_1}} \right)}& \cdots &{{\varphi _N}\left( {{\mathit{\boldsymbol{u}}_1}} \right)}\\ \vdots & \ddots & \vdots \\ {{\varphi _1}\left( {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{I}}}} \right)}& \cdots &{{\varphi _N}\left( {{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{I}}}} \right)} \end{array}} \right] $

(6)

这里N是多项式基函数的个数。再引入对密度函数小波分解的小波基函数矩阵：

$ \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \left[ {\begin{array}{*{20}{c}} {{\psi _1}\left( {{y_1}} \right)}& \cdots &{{\psi _1}\left( {{y_K}} \right)}\\ \vdots & \ddots & \vdots \\ {{\psi _M}\left( {{y_1}} \right)}& \cdots &{{\psi _M}\left( {{y_K}} \right)} \end{array}} \right] $

(7)

这里M是小波基函数个数。由此可以建立过程操作条件和晶体产品CSD之间的函数关系为：

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} + \mathit{\boldsymbol{E}} $

(8)

式(8)中：E是噪音矩阵，V是系数矩阵。注意Φ和Ψ是已知确定的基函数，只有系数矩阵V未知待求解。为了估计系数矩阵V，可以采用如下最小二乘法。

$ \min {\left\| {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right\|^2} $

(9)

从而得到系数矩阵V如式(10)。

$ \mathit{\boldsymbol{V}} = {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right)^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right)^{ - 1}} $

(10)

因此，将连续管式结晶器的振荡频率、2个区段的降温速率作为操作变量，通过式(8)可以预测晶体产品CSD。

3 晶体产品CSD评价指标

结晶产品CSD对其堆密度和流动性有关联作用，进而会影响下游工艺(过滤、干燥、造粒和压片)的效率^[22]。现有文献常使用产品CSD的平均尺寸与变异系数作为评价指标^{[23, 24]}。但考虑到实际生产中期望产品晶体集中在一定的目标尺寸范围内，为此定义CSD曲线以目标尺寸为中心的面积评估指标与总面积计算公式分别为：

$ A_0=\int_{d-r / 2}^{d+r / 2} \gamma\left(\boldsymbol{u}_i, y_k\right) \mathrm{d} y_k $

(11)

$ A=\int_0^a \gamma\left(\boldsymbol{u}_i, y_k\right) \mathrm{d} y_k $

(12)

式(11)和(12)中：d是目标尺寸；r是目标尺寸区间宽度；a是CSD曲线主要分布的区间宽度。即产品收率可以计算为：

$ Y=\frac{A_0}{A} $

(13)

实际工程应用中通常期望晶体产品有较小的变异系数，这样有利于下游操作。产品CSD曲线的方差可计算为：

$ \sigma^2=\frac{\sum\limits_{k=1}^K y_k^2 \gamma\left(\boldsymbol{u}_i, y_k\right)}{\sum\limits_{k=1}^K \gamma\left(\boldsymbol{u}_i, y_k\right)}-(\bar{y})^2 $

(14)

式(14)中y可以表示为：

$ \bar{y}=\frac{\sum\limits_{k=1}^K y_k \gamma\left(\boldsymbol{u}_i, y_k\right)}{\sum\limits_{k=1}^K \gamma\left(\boldsymbol{u}_i, y_k\right)} $

(15)

综上分析，这里提出一个综合评价指标：

$ Q=\frac{1}{Y}+\alpha \sigma^2 $

(16)

因此，最优的过程操作条件可以由下面的优化问题确定：

$ \min _u(Q)=\min _u\left(\frac{1}{Y}+\alpha \sigma^2\right) $

(17)

式(16)和(17)中α是用于标称化2项数值量级的加权因子。

4 实验结果与分析 4.1 L-谷氨酸连续结晶实验 4.1.1 实验设备

采用DN15 COBC实验装置，该结晶器使用内径为15 mm、管段总长14 m的COBC，其直管尺寸如图 3所示。夹套温度区段设置如图 1所示。

溶质是β-LGA(C₅H₉NO₄, 质量分数99%, Sigma公司生产)，溶剂是蒸馏水。实验使用具有千分之一精确度的高精度天平(Mettler-Toledo)来称量β-LGA的质量。

采用FBRM (ParticleTrack G400, Mettler-Toledo生产)测量结晶过程采集产品的晶体弦长分布(CLD)。探测器设置为精细模式，在实验前激光聚焦于0 μm，扫描速度设置为2 m·s^-1。CLD直方图由1~1 000 μm范围内100个方格组成，并且每30 s采集1次数据。

4.1.2 晶种制备

晶种制备采用研磨、筛分和洗涤的方法。原料晶体首先在标准研钵中研磨20 min，然后用筛分振动器(AS 200 digit，由Retsch GmbH生产)筛分90 min。筛网的尺寸从上到下分别为120、100、90和70 μm。收集在70和90 μm之间的晶体，最后使用蒸馏水在25 ℃下清洗3 min，以去除细小的颗粒。

4.1.3 实验步骤

首先，将原料晶体溶解在蒸馏水中，配制45 ℃的饱和浓度。其次，将溶液加热并保持在55 ℃约30 min，以确保溶质的完全溶解。然后以1 ℃·min^-1的速度将溶液冷却到45 ℃。最后根据不同的操作条件进行实验。总流速设定为80 g·min^-1，停留时间为62.5 min。使用FBRM得到实验产品的CLD。

4.2 实验设计

基于表 1的操作条件范围与表 2的实验设计方案，选取合适的输入变量如表 3所示。

其他实验条件设定为常值，包括振荡器振幅25 mm，主料液流速60 g·min^-1，晶种悬浊液流速20 g·min^-1，晶种装载率20%。

4.3 双层基函数预测建模 4.3.1 基函数选取

由于Daubechies小波具有正交性和紧支撑性等优点，本论文应用DB4作为第一类基函数。基于表 3的过程操作条件，为了更好的描述输入u与CSD曲线的关系，采用主成分分析法选取三元多项式基函数作为第二类基函数，即：

$ \begin{aligned} & \left\{1, u_1, u_2, u_3, u_1{ }^2, u_2{ }^2, u_3{ }^2, u_1 u_2, \right. \\ & \left.u_2 u_3, u_1 u_3, u_1{ }^2 u_2, u_1 u_2{ }^2, u_1{ }^2 u_3\right\} \end{aligned} $

(18)

4.3.2 模型训练结果

为了避免输入向量的完全共线性，表 2的中心点重复实验C只做1次。通过训练表 3的过程操作条件，得到的双层基函数拟合结果如图 4所示, 其中N为实验批次号。可见模型预测对实验测量达到较好的拟合效果。

4.4 过程操作条件优化 4.4.1 目标函数求解

式(11)中d为160, r为80。式(12)中的a为400。式(15)中的α为10^-3。使用粒子群算法(Particle swarm optimization, PSO)对式(16)的优化问题进行求解，得到的结果如表 4所示。

4.4.2 模拟与实验验证

采用上述双层基函数预测模型得出的结果如图 5所示。

可以看到，采用优化的过程操作条件保证CLD方差的同时，由式(13)计算得到收率提高8.2%。实验验证结果如图 6所示。

5 结论

提出一种基于响应面实验设计与基于双层基函数的晶体产品尺寸分布预测建模方法，由此给出过程操作条件优化方法。基于不完全因子数的实验设计能有效减少实验次数，相应建立的数据驱动预测模型实现准确地预测晶体产品CSD。通过引入一个关于目标晶体尺寸分布和产品收率的目标函数，利用CSD预测模型可以优化COBC结晶器的振荡频率以及两个区段的降温速率，通过对β型L-谷氨酸冷却结晶过程的仿真和实验，验证了该方法的有效性和优点。