间歇冷却结晶过程的CFD-PBM数值模拟

2023-05-26覃渫兰,张磊,都健

大连理工大学学报 2023年3期

覃渫兰, 张磊, 都健

(大连理工大学化工学院化工系统工程研究所, 辽宁大连 116024)

0 引言

溶液结晶是化学、制药和食品工业中应用较为广泛的分离和纯化技术之一[1].在结晶过程中,由于对物理化学和机械特性的强烈依赖性,需要控制结晶产品质量[2].粒度分布(crystal size distribution,CSD)是评价最终晶体质量的重要指标之一[3-4].

自20世纪60年代中期以来,粒数衡算模型(population balance model,PBM)已被广泛用于结晶过程的模拟[1-2].粒数衡算模型研究的是既包括外部坐标(空间)又包括内部坐标(颗粒属性)的分散系统[1].离散方法是应用较为广泛的数值方法,该方法首先将整个连续范围划分为若干个相邻的小区间,然后将粒数衡算方程转换为一定数量的离散方程进行求解[2].尽管该方法计算成本较高,但其优点在于数值鲁棒性好[5].

实际结晶系统多为固液两相流体系,流场分布和结晶动力学等各个过程共同决定最终产品质量,进而影响晶体粒度分布[6-7].计算机模拟是研究结晶过程的一种有效方法,将计算流体动力学(computational fluid dynamics,CFD)和粒数衡算模型相结合,以解释流场和物质浓度场中强烈的空间不均匀性[2].

关于CFD-PBM耦合模型描述结晶过程的研究已取得了一定成果.2018年,Fu等[8]建立了CFD-PBM耦合模型来模拟对乙酰氨基酚的连续冷却结晶过程,结合生长动力学和成核动力学,探究混合效应和壁面温度对颗粒粒度分布的影响.2019年,Mousavi等[9]对间歇式搅拌结晶器中磷酸铵镁晶体的结晶过程,建立了CFD-PBM耦合模型,仅考虑晶体生长动力学,模拟计算的磷酸铵镁晶体的粒度分布和实验值相差不大.2021年,De Souza等[6]通过CFD-PBM耦合模型模拟了间歇式搅拌结晶器中硫酸铝钾的冷却结晶过程,仅考虑晶体生长动力学,探究实际流动条件对晶体生长动力学的影响.

关于CFD-PBM耦合模型模拟结晶过程的研究中,大多对结晶动力学模型进行简化,关于破碎过程的研究工作有限.然而,结晶过程是一个复杂的动态相平衡过程,晶体的破碎行为时有发生,特别是在搅拌结晶器中,颗粒的破碎会影响离散相的混合和流体流动行为[10].为了准确分析结晶过程,需要更全面地考虑晶体行为[11].

综上所述,本文针对搅拌结晶器中扑热息痛-乙醇体系的间歇冷却结晶过程,采用CFD-PBM耦合模型进行模拟与优化,充分考虑生长、初级成核和二次成核过程,同时采用Luo破碎模型补充结晶动力学模型,以探究晶种加入、破碎过程、搅拌速度和降温速率对晶体粒度分布的影响.

1 数学模型

1.1 粒数衡算模型

粒数衡算方程(population balance equation,PBE)描述了多相体系中离散相状态在时间和物理空间上的分布[12-13].基于颗粒特征尺寸L的粒数衡算模型的表达式如下所示:

∇[Γeff∇n(L,y,t)]+ρl∂[G(L)n(L,y,t)]/∂L=

(1)

晶体的破碎速率为

g(L)n(L,t)

(2)

式中:g为破碎频率;β(L,L′)为从颗粒特征尺寸L′到L的破碎颗粒概率密度函数;P为一个颗粒破碎产生子颗粒的数量.

1.2 CFD-PBM耦合模型

CFD是一种利用流体基本控制方程求解流体流动的数值方法[14].CFD-PBM耦合模型结合计算流体力学模型和粒数衡算模型,即将空间尺度的流场分布和时间尺度的颗粒状态演变相结合,从而得到产品晶体的粒度分布.CFD模拟计算得到的颗粒流动特性和流场分布,用于求解PBM的成核速率、生长速率和破碎速率等,同时求解PBM方程得到颗粒粒度分布和颗粒直径,进而修正和计算相间作用力和湍动能等,改进CFD模型的预测能力,因此CFD-PBM耦合是双向的[7].

为了将离散相的PBM与CFD结合,采用Sauter平均直径表示离散相的颗粒直径.对于离散方法,Sauter平均直径定义为

(3)

如图1所示,在本文研究中,CFD-PBM耦合模型主要考虑并结合以下建模层次[7,15]:

图1 CFD-PBM耦合模型

(1)计算流体动力学模型,描述颗粒在结晶器中的运动和分布.

(2)粒数衡算模型,描述所有颗粒相关性质在时间和空间上的演化.

(3)结晶动力学模型,描述结晶过程中晶体的行为.

2 结晶动力学模型

2.1 结晶过程

(1)晶体成核过程

成核速率是指单位时间单位体积溶液中新生成的微小晶体粒子数目[16].成核过程通常分为初级成核(primary nucleation)和二次成核(secondary nucleation)[3,17].

初级成核是在无晶种加入时自发成核[3,17].初级成核通常被经典成核理论(classical nucleation theory,CNT)所描述,初级成核速率b1的表达式如下[18]:

(4)

式中:kb1为初级成核常数;V为一个溶质分子的体积;k为玻尔兹曼常数;σ为晶体和溶液的界面能;T为温度;s为相对过饱和度,s=c/cs,c为扑热息痛在乙醇中的溶解度,cs为扑热息痛在乙醇中的饱和溶解度.

二次成核是在拥有晶体时形成晶体[3,17].结晶器中二次成核的主要原因是晶体与搅拌器的接触和碰撞,碰撞的概率与搅拌器的转速成正比[18].二次成核速率的经验方程为[19]

(5)

式中:kb2为二次成核常数;ms为晶体质量;α为常数2.08;β为常数0.713.

(2)晶体生长过程

晶体生长的机理有表面能理论和吸附层理论等[3].由于晶体生长的复杂性,至今仍未建立通用的生长理论[3].在工业结晶中,生长速率常采用幂函数型的经验方程,表达式如下[3,19]:

G=kgexp(-Ea/RT)(Δc)γ

(6)

式中:kg为生长速率常数;Ea为晶体活化能;Δc为绝对过饱和度,Δc=c-cs;γ为动力学常数;R为摩尔气体常数.

晶体破碎过程与颗粒性质和扰动情况有关[20].Luo破碎模型是常用的破碎模型,适用于湍流扩散中颗粒破碎的预测.Luo破碎模型是基于各向同性和概率理论建立的,可以预测给定颗粒尺寸的破碎速率.在湍流中,由于涡流在颗粒表面的轰击,离散相表面产生相对速度的波动,当涡流的能量达到某个阈值时,可以满足破碎所需要的表面能量,从而产生破碎过程[21].Luo破碎模型进行了一些假设[21]:湍流被认为是各向同性的;只考虑颗粒的二元破碎;破碎体积分数是一个随机变量;破碎的发生是由涡旋能量等级决定的;只有长度小于或等于颗粒直径的涡流才能引起粒子振荡.

Luo破碎模型的破碎速率为

Bbre-Dbre=Ωbre(L,L′)=(0.923 8ε1/3L-2/3ω)×

2.047-1}ξ-11/3}dξ

(7)

式中:λ为涡流大小;ξ为量纲一涡流大小,ξ=λ/L;ε为离散相的局部分数;f为破碎常数;ρp为颗粒密度;ω为常数1.5.

2.2 扑热息痛冷却结晶的动力学模型

Li等[19]利用实验数据回归了扑热息痛-乙醇结晶过程的生长和成核动力学.结晶实验的装置如图2所示,其中结晶系统包括一个温度探头、搅拌系统、1 L玻璃结晶器和加热冷却金属夹套.搅拌器为4个桨叶搅拌器.

图2 实验装置示意图[19]

实验中,初始溶液体积为0.5 L,加入的晶种质量为0.03 g,晶种粒径为75～106 μm,结晶器的搅拌速度为400 r/min,采用式(8)的线性降温曲线.表1为涉及的物性参数.

表1 物性参数表

(8)

扑热息痛在乙醇中的饱和溶解度为

在大一新生入学时，重视宣传学科竞赛的重要性，组织大一新生赴现场观摩赛事，通过课堂表现和课外接触发现可造之材，鼓励有理想有潜力的新生参与到赛事之中，感受赛事的紧张程度，体验成就感和失落感，为日后成为参赛骨干打下基础。

cs=7.915×10-7T3-6.439×10-4T2+

1.765×10-1T-16.17

(9)

扑热息痛结晶过程的生长速率G、初级成核速率b1和二次成核速率b2公式分别如式(6)、(4)和(5)所示.式中动力学参数由Li等[19]的实验拟合得到,即生长速率常数kg为45.5 m/s,指数γ为1.24,晶体活化能Ea为41.3 kJ/mol,初级成核常数kb1为0.192 1 s-1/kg,晶体和溶液的界面能σ为4.25 mJ/m2,玻尔兹曼常数k为1.380 6×10-23J/K,二次成核常数kb2、α和β分别为1×105s-1/kg、2.08和0.713.

3 模拟与计算

3.1 三维建模与网格划分

通过Solidworks软件建立搅拌结晶器的几何图形,如图3所示.结晶器元件的相关尺寸见表2.ANSYS meshing生成的网格如图4所示.扭曲度(skewness)用于评价网格质量,即趋于理想网格的程度,取值为0～1,数值越大网格质量越差.基于归一化的正角度,扭曲度被定义为

表2 结晶器的尺寸

图3 简化的三维搅拌结晶器

(a)三维视图

(10)

式中:θmax为网格中最大角度,θmin为网格中最小角度,θe为理想网格角度.本文划分的网格最大扭曲度为0.80,平均扭曲度为0.23,网格质量良好.

3.2 数值模拟设置

采用ANSYS FLUENT进行CFD-PBM模拟仿真.模型设置:选择欧拉(Eulerian)多相流,湍流模型选择RNGk-ε模型,该模型引入了旋转、曲率相关项,能够模拟旋流等中等复杂的流动[22].对于搅拌桨旋转的模拟,采用滑移网格方法(sliding mesh method,SMM).降温曲线通过用户定义函数(UDF)在壁面上进行边界条件设置.PBM采用离散法进行求解,离散成30个区域,成核和生长动力学模型通过用户定义函数(UDF)实现.速度与压力的耦合算法为Coupled,残差设置为1×10-4,在单节点64核内存256 GB超算平台进行并行计算.

4 模拟结果与分析

4.1 模型验证

(1)网格独立性验证

通过网格独立性验证排除网格数量对模拟结果的影响.一共划分了网格数量为120 877、155 210、173 774和231 714的4套网格,对比不同网格数量的速度分布,如图5所示.最终173 774个网格被确定为精度和数值成本权衡之间的最佳选择.

(2)搅拌结晶器速度流场流型验证

图6(a)为ANSYS FLUENT软件的模拟结果,图6(b)为文献[23]中搅拌结晶器的计算结果,通过对比发现,两者的流场分布大致相符.结晶器的流场形成上下两个循环,使得物料在结晶器内得到充分的接触与混合.

(3)模拟粒度分布与实验值对比

在与Li等[19]相同结晶条件下,模拟结果与实验结果的对比如图7所示.图7(a)、(b)分别为体积密度函数V(L,y,t)和数量密度函数n(L,y,t)下的粒度分布.体积密度(volume density)函数是指单位物理空间(体积)、单位颗粒长度的颗粒体积.数量密度(number density)函数是指单位物理空间(体积)、单位颗粒长度的颗粒数.由图可知,粒度分布的模拟结果与实验结果偏差较小,且模拟得到的粒度分布与实验趋势基本一致,因此模拟结果与实验结果有较好的一致性.模拟值和实验值的偏差可能是结晶动力学模型的简化,以及CFD-PBM耦合模型所采用的数值方法带来的不确定性造成的.由于体积密度函数下的粒度分布更能体现颗粒分布的特性,后续研究中粒度分布都将采用体积密度函数下的粒度分布.

(a)体积密度函数下的粒度分布

4.2 搅拌结晶器的流场分析

图8为400 r/min搅拌速度下的颗粒体积分数φ分布.搅拌桨作为上下部分的分界线,各自形成循环流股.因为搅拌桨位置靠下,下部循环流股的流速比上部要快,所以下方颗粒被带起,使得搅拌桨下方的搅拌相对更均匀.

图8 400 r/min搅拌速度下的颗粒体积分数分布

图9为结晶器中不同X坐标下颗粒体积分数的纵向分布曲线图,其中X为距离搅拌桨中心的水平距离.由图可知,颗粒体积分数随着高度增加呈阶梯式分布.除了距离搅拌桨中心较近的区域外,其他区域均存在两个阶梯.从结晶器底部离开,颗粒体积分数会出现骤降,然后进入一个平台,即颗粒分布均匀区域,这是由于距离搅拌桨较近,湍流强度较大,混合效果较好.

图9 不同X坐标下颗粒体积分数纵向分布

4.3 晶种的影响

图10为加入晶种和未加入晶种的模拟结果,加入的晶种质量为0.03 g,晶种粒径为75～106 μm.由图可知,加入晶种的粒度分布峰值比未加入晶种的高,即产生的颗粒数更多.同时加入晶种获得的颗粒平均粒径256.20 μm大于未加入晶种的平均粒径254.95 μm.由于加入晶种数量较少,因此两者偏差较小.为了获得颗粒数多、平均粒径大的粒度分布,可以适当加入晶种.

图10 加入晶种和未加入晶种的粒度分布对比

图11为加入晶种和未加入晶种条件下,生长和成核过程产生的源项Sφ随时间的变化曲线,即纵坐标体现了晶体生长速率和成核速率的大小.由图可知,在刚开始结晶过程中,加入晶种的源项高于未加入晶种的源项,这是由于自发成核的过饱和度高于颗粒生长所需的过饱和度.在5 000 s之前,源项曲线呈现快速增长趋势,是颗粒生长和成核过程的关键阶段,该阶段中加入晶种的源项始终大于未加入晶种的源项,即加入晶种的生长速率和成核速率更大,产生的颗粒数更多、平均粒径更大.图12为在3 000 s时颗粒体积分数分布图,其中图12(c)为加入晶种和未加入晶种的颗粒体积分数差值的分布.图中大部分区域的数值为正,说明结晶器内加入晶种的颗粒体积分数大于未加入晶种的.所以在结晶器内,加入晶种条件下,颗粒数更多,颗粒的体积分数更大.

图11 加入晶种和未加入晶种条件下结晶动力学源项随时间的变化曲线

(a)加入晶种

因此,在结晶过程中,加入晶种可以避免新晶体成核造成的高过饱和度,从而获得具有颗粒数多、平均粒径大的粒度分布晶体.

4.4 破碎过程的影响

Luo破碎模型是常用的破碎模型,适用于湍流中颗粒破碎的预测.图13为在考虑破碎和不考虑破碎条件下颗粒的粒度分布以及与实验结果的对比.由图可知,考虑破碎的粒度分布峰值更大,颗粒数更多,更接近于实验结果.因此,补充Luo破碎模型可以减少由于结晶动力学模型简化带来的偏差.

图14为结晶动力学源项Sφ随时间的变化曲线,其中两者差值在一定程度上体现了Luo破碎模型破碎速率的大小.由图可知,考虑破碎的结晶过程粒度分布峰值更大,主要原因是前期考虑破碎的结晶动力学源项明显高于不考虑破碎的,即由于补充Luo破碎模型,颗粒属性受到了破碎速率的影响,颗粒数更多.结晶后期,破碎速率的影响逐渐减小,两者源项相差较小.因此采用合适的破碎模型,考虑更为全面的结晶动力学模型,能够获得提高最终颗粒粒度分布预测的能力.

图14 不同破碎条件下结晶动力学源项随时间的变化曲线

4.5 搅拌速度的影响

图15为不同搅拌速度下颗粒体积分数的分布.由图可知,随着搅拌速度的增大,传递给流体的机械能越来越大,结晶器内的固相体积分数分布趋于均匀.在搅拌速度为300、400 r/min时,结晶器中产品颗粒的分布在上下两部分有明显差异,搅拌桨作为上下部分的分界线,各自形成循环流股.在搅拌速度为500 r/min时,结晶器内的颗粒体积分数分布基本均匀.

(a)200 r/min

图16为不同搅拌速度下粒度分布的对比.由图可知,随着搅拌速度增大,粒度分布峰值增加,颗粒数增加,其中300、400 r/min的峰值相差不大.由图17可知,随着搅拌速度Nr增大,颗粒的最终平均粒径减小.主要是因为随着搅拌速度的增大,结晶器中的流体逐渐混合均匀,有利于结晶器内的传热过程,进而使得溶液的平均温度下降,过饱和度略有增加,导致成核速率增大;同时搅拌速度较大时,湍流耗散率较高,进一步加强了二次成核,产生更多细小晶体,所以颗粒数增加,平均晶体尺寸减小.因此搅拌程度对于粒度分布、平均粒径等性质有较大的影响,为了得到理想的颗粒产品,必须对结晶工艺的搅拌速度进行合理的设计和控制.

图16 不同搅拌速度下粒度分布的对比

图17 平均粒径随搅拌速度的变化曲线

4.6 降温速率的影响

线性降温曲线公式为

(11)

采用不同降温速率k对结晶过程进行模拟,如图18所示.当降温速率k增大时,粒度分布峰值增加,即颗粒数增加,并且粒度分布变窄,其中降温速率k为-1/150 ℃/s和-1/200 ℃/s的峰值相近.由图19可知,随着降温速率k的增大,颗粒的平均粒径逐渐减小,虽然降温速率k从-1/150 ℃/s到-1/200 ℃/s的平均粒径略有增大,但是相差很小.因此适当采用较为平缓的线性降温曲线,能够获得更窄的粒度分布以及更多的颗粒数,但是产品最终的平均粒径会更小.所以在结晶工艺中控制降温速率是很有必要的.