蛋氨酸(D, L-Methionine)又称甲硫氨酸，是人体必须的氨基酸之一，参与蛋白质的合成，被广泛应用于饲料、生物医药和化妆品等行业^[1]。各种蛋氨酸生产工艺都需将蛋氨酸进行结晶精制^[2]。目前已有较多学者对蛋氨酸的结晶动力学、热力学性质进行研究^{[3, 4]}，但对蛋氨酸的连续结晶设备研究较少。在蛋氨酸结晶过程中，晶体难以各向均匀生长，容易造成形状不规则，因此在工业中蛋氨酸的生长需要分布均匀的过饱和度；另一方面，蛋氨酸晶体析出速率过快会导致晶体颗粒细小，因此需要避免生产过程中因过饱和度过高而产生大量晶核。DTB(Draft tube baffled)结晶器带有搅拌装置，可以为蛋氨酸的生产提供分布均匀的生长环境，且作为一种典型的晶浆内循环型结晶器，DTB结晶器可以不断地消耗过饱和度，使之维持在较低水平，因此，工业上常选用DTB结晶器作为蛋氨酸的生产设备。

随着数值模拟的高速发展，计算流体力学(CFD)被广泛应用于DTB结晶器的研究，以充分了解流体在DTB结晶器中的流动特性和颗粒悬浮状态，有利于对DTB结晶器进行结构和操作条件的优化^[5-7]。除此之外，利用CFD对结晶过程^[8-10]和DTB结晶器的放大^{[11, 12]}进行研究，可以大幅度降低时间成本，对结晶器的设计与结晶过程的控制具有重要意义。

内循环是DTB结晶器的重要特征，在蒸发冷却结晶过程中，搅拌桨推动流体沿导流筒到达沸腾液面，是DTB结晶器内循环的动力来源，因此许多学者针对结晶器的搅拌桨进行研究^{[6, 12-14]}，希望通过改进搅拌桨的型式提高晶浆密度，降低搅拌速度，减少过剩晶核的形成。除此之外，导流筒是保证结晶器在较低的压头损失(约100~200 mm H₂O)下就能实现良好内循环的关键结构^[15]，尤其在蒸发结晶和真空蒸发冷却结晶过程中，导流筒将大量晶浆送到沸腾液面，逐渐消耗溶液产生的过饱和度，避免了爆发成核，也降低了结晶疤的产生速率，有利于连续结晶的稳定运行；环形挡板将结晶器分为晶体生长区和澄清区，可以有效控制产品的粒度和粒径分布，对工业生产有重要作用。

目前DTB结晶器的研究主要集中在实验室规模，少有学者关注工业级DTB结晶器的结构特性和流场分布情况，因此，本研究通过CFD技术探究增加导叶结构提高内循环的可行性，研究导流筒和挡板长度对工业级DTB结晶器的流场特性和颗粒悬浮状态的影响，为工业结晶器的设计提供参考。

1 DTB结晶器的模型建立与计算方法 1.1 DTB结晶器模型的简化

在实际工业生产过程中，为方便操作，DTB结晶器常设有支架、人孔等辅助结构，这些结构对DTB结晶器的内部流场不会造成明显影响，因此本研究在建模过程中将这些辅助结构忽略^[16]。模拟的DTB结晶器具体结构如图 1所示，该结晶器总体积为26.05 m³，结晶器直径为3.2 m，高度为5.4 m，进料口、物料循环出口和产品采出口的直径都为0.1 m，搅拌桨和导叶直径为0.76 m，导流筒直径为0.8 m，挡板直径为2.4 m。

1.2 湍流模型

标准k-ε模型忽略分子黏性对湍流流动的影响，通过湍动黏度描述湍流动能k和湍流耗散率ε，适用于湍流充分发展的流动，适用范围广，计算结果精确，计算规模相对较小^[17-19]。

式(1)和式(2)分别表示湍流动能k方程和湍流耗散率ε方程^[20]。

$ \begin{gathered} \frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_i}\left(\rho k u_i\right)=\frac{\partial}{\partial x_j}\left[\left(\mu+\frac{\mu_{\mathrm{t}}}{\sigma_k}\right) \frac{\partial k}{\partial x_j}\right]+ \\ G_{\mathrm{k}}+G_{\mathrm{b}}-\rho \varepsilon-Y_{\mathrm{M}}+S_k \end{gathered} $

(1)

$ \begin{gathered} \frac{\partial}{\partial t}(\rho \varepsilon)+\frac{\partial}{\partial x_i}\left(\rho \varepsilon u_i\right)=\frac{\partial}{\partial x_j}\left[\left(\mu+\frac{\mu_{\mathrm{t}}}{\sigma_k}\right) \frac{\partial \varepsilon}{\partial x_j}\right]+ \\ C_{1 \varepsilon} \frac{\varepsilon}{k}\left(G_{\mathrm{k}}+C_{3 \varepsilon} G_{\mathrm{b}}\right)-C_{2 \varepsilon} \rho \frac{\varepsilon^2}{k}+S_{\varepsilon} \end{gathered} $

(2)

$ \mu_{\mathrm{t}}=\rho C_\mu \frac{k^2}{\varepsilon} $

(3)

$ G_{\mathrm{k}}=\mu_{\mathrm{t}} \frac{\partial u_i}{\partial u_j}\left(\frac{\partial u_i}{\partial u_j}+\frac{\partial u_j}{\partial u_i}\right) $

(4)

式(1)~式(4)中：t为时间；x为空间位置；ρ为密度；μ_t为湍动黏度，可以用式(3)表示；G_k是由平均速度梯度产生的湍流动能，如式(4)所示；G_b为浮力产生的湍流动能；Y_M是可压缩流体湍流中脉动扩张的湍流动能；σ_k和σ_ε分别为k和ε的湍流Prandtl数；C_1ε、C_2ε、C_μ为经验常数；S_k和S_ε表示用户自定义源项。

在Fluent中，C_1ε=1.44，C_2ε=1.92，C_μ=0.09，σ_k=1.0，σ_ε=1.3，当流体不可压缩且无自定义源项时，G_b=0，Y_M=0，S_k=0，S_ε=0。

1.3 欧拉模型的基本理论

在研究晶体颗粒在液体中的悬浮状态时无法忽略固、液两相之间的相互作用，因此选择欧拉-欧拉模型对固、液两相进行模拟分析。在欧拉-欧拉模型中，每一相都被看作是连续的，每一相的体积分数都可以由连续性方程求得：

$ \frac{\partial}{\partial t}\left(\alpha_{\mathrm{q}} \rho_{\mathrm{q}}\right)+\nabla\left(\alpha_{\mathrm{q}} \rho_{\mathrm{q}} \overrightarrow{\boldsymbol{v}}_{\mathrm{q}}\right)=\sum\nolimits_{p=1}^n\left(\dot{m}_{\mathrm{pq}}-\dot{m}_{\mathrm{qp}}\right) $

(5)

式(5)中：α为液相或固相的体积分数；$ \vec{\boldsymbol{v}}$为流体速度矢量；ṁ_pq是p相向q相的传质速率；ṁ_qp是q相向p相的传质速率。由于本研究不涉及到相间传质作用，因此传质速率ṁ_pq和ṁ_qp都为0。

欧拉-欧拉模型中液相和固相的动量方程需要分别定义^[21-23]。液相q的动量方程可以写为：

$ \begin{gathered} \frac{\partial}{\partial t}\left(\alpha_{\mathrm{q}} \rho_{\mathrm{q}} \overrightarrow{\boldsymbol{v}}_{\mathrm{q}}\right)+\nabla \cdot\left(\alpha_{\mathrm{q}} \rho_{\mathrm{q}} \overrightarrow{\boldsymbol{v}}_{\mathrm{q}} \overrightarrow{\boldsymbol{v}}_{\mathrm{q}}\right)=-\alpha_{\mathrm{q}} \nabla \boldsymbol{p}+\nabla \cdot \overline{\overline{\boldsymbol{\tau}}}_{\mathrm{q}}+ \\ \alpha_{\mathrm{q}} \rho_{\mathrm{q}} g+\sum\nolimits_{p=1}^n\left[K_{\mathrm{pq}}\left(\overrightarrow{\boldsymbol{v}}_{\mathrm{p}}-\overrightarrow{\boldsymbol{v}}_{\mathrm{q}}\right)+\dot{m}_{\mathrm{pq}} \overrightarrow{\boldsymbol{v}}_{\mathrm{pq}}-\dot{m}_{\mathrm{qp}} \overrightarrow{\boldsymbol{v}}_{\mathrm{qp}}\right]+ \\ \left(\vec{F}_{\mathrm{q}}+\vec{F}_{\mathrm{lift}, \mathrm{q}}+\vec{F}_{\mathrm{vm}, \mathrm{q}}\right) \end{gathered} $

(6)

固相p的动量方程可以写为：

$ \begin{gathered} \frac{\partial}{\partial t}\left(\alpha_{\mathrm{p}} \rho_{\mathrm{p}} \overrightarrow{\boldsymbol{v}}_{\mathrm{p}}\right)+\nabla \cdot\left(\alpha_{\mathrm{p}} \rho_{\mathrm{p}} \overrightarrow{\boldsymbol{v}}_{\mathrm{p}} \overrightarrow{\boldsymbol{v}}_{\mathrm{p}}\right)=-\alpha_{\mathrm{p}} \nabla \boldsymbol{p}-\nabla \boldsymbol{p}_{\mathrm{p}}+ \\ \nabla \cdot \overline{\overline{\boldsymbol{\tau}}}_{\mathrm{p}}+\alpha_{\mathrm{p}} \rho_{\mathrm{p}} g+\sum\nolimits_{q=1}^n\left[K_{\mathrm{qp}}\left(\overrightarrow{\boldsymbol{v}}_{\mathrm{q}}-\overrightarrow{\boldsymbol{v}}_{\mathrm{p}}\right)+\right. \\ \left.\dot{m}_{\mathrm{qp}} \overrightarrow{\boldsymbol{v}}_{\mathrm{qp}}-\dot{m}_{\mathrm{pq}} \overrightarrow{\boldsymbol{v}}_{\mathrm{pq}}\right]+\left(\vec{F}_{\mathrm{p}}+\vec{F}_{\mathrm{lift}, \mathrm{p}}+\vec{F}_{\mathrm{vm}, \mathrm{p}}\right) \end{gathered} $

(7)

其中：$ \overline{\overline{\boldsymbol{\tau}}}_q$为各相的应力张量；K_pq为两相之间的动量传递系数；K_pq=K_qp，$ \vec{F}_{\mathrm{p}}, ~ \vec{F}_{\text {lift }, \mathrm{p}} \text { 和 } \vec{F}_{\mathrm{vm}, \mathrm{p}}$分别为外部质量力、升力和虚拟质量力。

1.4 相间作用力

多相流中不同相之间的相互作用力有虚拟质量力、曳力、升力和湍流分散力等，在实际工业结晶过程中，液相与固相之间的曳力占主导作用。因此本研究主要考虑两相之间的曳力作用。

选择精度较高的Syamlal-O′Brien模型^[24]计算固体颗粒在连续液相中受到的曳力，其表达式为：

$ f=\frac{C_{\mathrm{D}} R e_{\mathrm{s}} \alpha_1}{24 v_{\mathrm{r}, \mathrm{~s}}^2} $

(8)

式(8)中C_D为Dalla Valla推导所得的曳力系数^[25]：

$ C_{\mathrm{D}}=\left(0.63+\frac{4.8}{\sqrt{\operatorname{Re}_{\mathrm{s}} / v_{\mathrm{r}, \mathrm{~s}}^2}}\right)^2 $

(9)

v_{r, s}为与分散相相关的最终速度^[26]：

$ \begin{gathered} v_{\mathrm{r}, \mathrm{~s}}=0.5\left(A-0.06 R e_{\mathrm{s}}+\right. \\ \sqrt{\left(0.06 R e_{\mathrm{s}}\right)^2+0.12 R e_{\mathrm{s}}(2 B-A)+A 2} \end{gathered} $

(10)

$ A=\alpha_1^{4.14} $

(11)

$ \left\{\begin{array}{l} \alpha_1 \leqslant 0.85, \quad B=0.8 \alpha_1^{1.28} \\ \alpha_1>0.85, \quad~~ B=\alpha_1^{2.65} \end{array}\right. $

(12)

Re_s为相对雷诺数^[27]，其速度由液固两相的相对瞬时速度表示：

$ R e_{\mathrm{s}}=\frac{\rho_1 d_{\mathrm{s}}\left(v_{\mathrm{s}}-v_1\right)}{\mu_1} $

(13)

式(13)中：下标l和s分别表示液相和固相；d_s为固体颗粒的直径。

1.5 数值模拟方法

DTB结晶器壁面与底面设置为wall，顶部假设为静液面，设置为symmetry，进料口设置为velocity-inlet，产品采出管口和物料循环管口设为pressure-outlet，由于DTB结晶器带有搅拌桨，考虑到搅拌桨与结晶器之间的相对运动，本研究采用多重参考系法(MRF)^[28]将搅拌区域设置为旋转区域，搅拌桨相对搅拌区的速度为0，其它区域设置为静区域。

模拟过程中以水相为连续相，密度为998.2 kg ·m^-3，黏度为0.001 Pa ·s；蛋氨酸颗粒为分散相，密度为1 340 kg ·m^-3，黏度为0.004 4 Pa.s，粒径为0.3 mm，假设流体为不可压缩流体。表 1为其它模拟参数设置。本研究模拟时忽略晶体成核生长过程，因此结晶器内部颗粒全部源于进口管，进口物料为3 m ·s^-1、蛋氨酸体积分数为12%的固-液混合物，搅拌桨转速设为400 r ·min^-1。

1.6 网格划分及无关性检验

利用ICEM软件对DTB结晶器模型进行四面体非结构网格划分，图 2为结晶器导入ICEM后的结构简图，以结晶器底部中心为坐标原点，在确定结晶器主体网格后，对变化梯度较大的搅拌桨、进料管、产品采出管和物料循环管进行网格加密处理，从而得到最优网格数量。

将DTB结晶器的模型按照合适的尺寸划分为4种不同的网格数：316万、626万、806万和1 026万，然后利用Fluent 2022 R1软件计算相同物性参数和操作条件下流体的运动状态，待流场稳定后选择靠近搅拌桨附近流动剧烈的区域进行监控，监测区如图 3所示，读取该处的速度值并绘制成图，结果如图 4所示。由图 4可知，806万网格与1 026万相比，在轴向和径向方向的速度几乎重合，可以保证所需的计算精度，因此本论文选择网格数量为806万的计算模型，最终划分的网格如图 5所示。

2 结果与分析 2.1 导叶结构对流动特性的影响

DTB结晶器中导流筒下端搅拌桨以较低转速旋转，可以看作内循环轴流泵装置^[29]，而导叶作为轴流泵中的一个重要部件，可将搅拌桨叶轮流出的流体由旋转运动转为轴向运动，对轴流泵的高效运行起着重要的作用，因此借鉴轴流泵的结构特点，本研究在导流筒内增加如图 6所示的导叶结构(导叶叶片数为6，倾斜角度为45°)，希望提高DTB结晶器的轴流效果，增大内循环。

图 7为DTB结晶器纵剖面的局部速度矢量图，图 8为纵剖面的体积分数云图。由图 7和图 8可见，常规DTB结晶器在搅拌轴附近的流场较为混乱且流速很慢，颗粒分布较少，部分流体在重力作用下向下流动，不利于晶浆的内循环；增加导叶结构后，通过导叶的流体在到达导流筒液面之前全部往上运动，不仅提高了结晶器的轴流效果，还增强了流场的稳定性。表 2的计算结果也证明导叶结构确实增加了DTB结晶器的内循环流量，更多的晶浆被送到沸腾液面，表层随时存在大量晶体，可以消耗不断产生的过饱和度进行晶体的生长，从而避免过饱和度过高发生爆发成核而降低产品粒度。

图 9和图 10分别为结晶器横剖面Y=2.2 m和Y=4.8 m处的固体体积分数分布图。

由图 9和图 10可知，导叶明显增强了颗粒的沉降效果、降低了循环出口的颗粒占比，这是由于导叶使内循环速度加快，从而降低了澄清区流体的运动速度，使颗粒未达到循环出口就发生沉降，但更多的晶体从产品出口采出，可能会加宽粒径分布。此外，晶体与导叶的碰撞可能会引发二次成核，影响产品粒径，因此实际工业过程中要根据生产目标，选择是否增加导叶结构。

2.2 导流筒长度对流动特性的影响

导流筒形成的内循环通道可以使DTB结晶器的各个流动截面维持较高的流动速度，因此对结晶器的流场分布具有重要的影响，本节以导流筒长度为研究对象(导流筒长度h为1.6、2.0、2.4、2.8和3.2 m)，探究DTB结晶器内的流场和颗粒悬浮状态。

图 11为不同导流筒长度时DTB结晶器的纵剖面矢量图。从图 11中可以看出，当导流筒长度较短时，内循环作用无法充分发挥，大量晶浆难以到达沸腾液面，无法消耗蒸发浓缩所产生的过饱和度；当导流筒长度增加时，循环路程变长，越来越多的流体被送到沸腾液面，为结晶过程提供较大的生长表面，有利于消耗不断产生的过饱和度，使之维持在较低的水平，避免因晶核生成量过多而影响产品粒度。但如果导流筒长度太长，距离沸腾液面过近，到达沸腾区的流体速度过大(如图 12和图 13所示)，不仅会增加颗粒与颗粒、颗粒与器壁的碰撞频率，导致晶体容易发生破碎现象，而且容易造成沸腾区流场不稳定，不利于过饱和度的稳定和均匀生成。

表 3列出了不同导流筒长度时结晶器内的流动特性。由表 3可知，导流筒长度在1.6~2.8 m范围内，结晶器的内循环流量随着导流筒长度的增加而增加，功率随之下降；但当导流筒长度超过2.8 m后，从导流筒出去的流体在参与漩涡之前先一步接触液面，动能损失较大，导致内循环流量下降，功率消耗增加。由此可以说明，导流筒长度太长不仅影响结晶过程的稳定运行，而且会降低结晶器的使用效率。

从导流筒出来的流体在参与内循环时容易产生漩涡，造成沸腾液面颗粒分布不均匀，如图 14所示，随着导流筒长度加长，导流筒与挡板之间的漩涡中心逐渐上移，颗粒分布不均现象得到改善。结晶器的横剖面图(图 15)显示，导流筒长度增加，澄清区颗粒的体积分数越来越小，这是由于循环阻力增大，使流体到达澄清区的速度减小，导致循环母液量减少，颗粒不能被带出，不利于细晶消除，这在工业上对生产粒度较大的产品是不利的，因此工业过程中要根据生产需要合理调整导流筒长度。

2.3 挡板长度对流动特性的影响

DTB结晶器中的挡板能控制流体向上的流速，使不同粒度的晶体颗粒实现分级，沉降速度较大的晶体颗粒发生沉降，从母液中分离出来，沉降速度较小的细晶随母液从循环出口排出，从而实现对微晶量的控制。因此本工作探究了挡板长度(挡板长度分别为1.8、2.0、2.2、2.4和2.6 m)对DTB结晶器内颗粒悬浮状态的影响。

图 16和图 17分别为平面X=0、Y=2.2 m处的固体颗粒体积分数分布云图，由图 16和图 17可知，挡板的长度对澄清区的颗粒分布影响不大，但随着挡板的加长，流体对结晶器锥型壁面的碰撞冲击力增大，颗粒容易在此处分布不均，造成局部沉积，晶体既不能从循环口排出进行微晶消除，也不能从产品出口采出得到晶体产品，这对工业生产非常不利。

表 4结果显示，挡板长度对DTB结晶器的功率消耗、内循环流量流量和外循环流量的影响不大，这是因为挡板位于澄清区，搅拌的影响已经消失，无法改变结晶器的内循环情况，因此只对晶体分级有影响。

图 18和图 19分别为X=0 m, Z=-1.4 m处的轴向液相速度分布图和颗粒体积分数分布图。由图 19可知，随着挡板长度的增加，澄清区循环出口的颗粒分布越来越少，这是因为沉降区域逐渐增大，母液到达循环出口的路径变长，溶液的主体速度下降，无法携带晶体从循环出口流出，更多的晶体产品发生沉降或从产品出口采出，因此不利于通过细晶消除得到粒径较大的产品。

除此之外，挡板长度过长，可能会加剧颗粒与器壁的碰撞，大颗粒在沉降之前可能会发生破碎，导致产品粒度变小。因此，在实际工业生产过程中，如果对产品的粒度要求比较高，不适合安装过长的挡板。

3 结论

利用欧拉-欧拉模型模拟获得了不同内部结构时DTB结晶器内流体的流动特性和颗粒的悬浮状况，得到结论如下：(1)在DTB结晶器的导流筒内增加导叶结构可以提高结晶器的轴流效果，增加内循环量，有利于晶体的均匀生长。(2)导流筒控制着DTB结晶器的内循环。增加导流筒长度可以为结晶过程提供更大的生长表面，避免产生过量的晶核。但导流筒过长，内循环产生的漩涡容易影响沸腾液面，造成二次成核，因此工业生产过程中要控制导流筒长度，使之维持在合适的长度范围。(3)挡板长度主要影响DTB结晶器的澄清区。挡板越长，晶体越容易在澄清区的锥形壁面局部沉积，不利于结晶器的稳定运行；而且随着挡板长度的增加，沉降区域逐渐增大，不利于细晶消除，得到较大粒径的产品。