DOI： 10.11817/j.issn.1672-7207.2019.08.031

溶液预冻方式与冰晶粒度分布之间的关系

贾颖姣^{1, 2},余云霞¹,刘志强¹,罗春¹

（1. 中南大学 能源科学与工程学院，湖南 长沙，410083；

2. 新奥泛能网络科技有限公司，北京，100176）

摘要:考虑到在真空冷冻干燥过程中，预冻阶段会通过影响冻结体冰晶粒度的分布最终对干燥过程的速率、能耗以及冻干产品品质产生重大影响，基于CFD软件对溶液在预冻过程中温度场变化的数值模拟结果，计算冻结后冰晶粒度，并研究控制搁板温度恒定、从常温匀速降温以及含养晶过程中恒温降这3种冻结方式对冰晶粒度分布的影响，分别得出各预冻参数（搁板温度、温降速率和养晶时间）与冰晶粒度分布参数（冰晶粒度平均值和标准差）之间的回归关系式，并且进行不确定度分析和显著性检验。研究结果表明：回归关系拟合优度和参数准确度较高，且多项式分布显著性水平达0.98以上，这对冻干曲线的优化有着一定的参考价值。

关键词:计算流体力学；粒度分布；冻结；均匀度；回归关系式

中图分类号:TB66    文献标志码:A     文章编号:1672-7207（2019）08-2026-07

Relationship between freezing methods of solution and ice grain size distribution

JIA　Yingjiao^{1, 2}, YU　Yunxia¹, LIU　Zhiqiang¹, LUO　Chun¹

(1. School of Energy Science and Engineering, Central South University, Changsha 410083, China；

2. ENN Ubiquitous Energy Network Technology Co. Ltd., Beijing 100176, China)

Abstract: Considering that in the vacuum freeze drying process, the pre-freezing stage will affect the distribution of ice crystal size of frozen material, which will ultimately affect the drying rate, energy consumption and the quality of freeze-dried products, based on the numerical simulation results using CFD of the temperature field changing during the prefreezing process, the size of ice grain size after freezing was calculated, and three kinds of freezing methods, i.e., maintaining shelf constant low temperatures, controlling the shelf to cool from room temperature with constant rates and cooling at constant rate containing crystal growth process, were studied to research the effects on the size distribution of the ice crystals. The formula of the fitting relationship between prefreezing parameter which included the temperature of shelf, the rate of temperature drop, the crystallization time and grain size distribution parameters of ice which included average value and uniforming of ice crystal were obtained, and the uncertainty analysis and the significance test were carried out. The results show that the fitted goodness and parameter accuracy are high and the significance level of polynomial distribution is above 0.98, which has reference for the optimization of freeze-drying curve.

Key words: computational fluid dynamics; particle size distribution; freezing; uniformity; regression formula

真空冷冻干燥是一种稳定的干燥过程。其原理如下：通过低温冷冻使液态水冻结为固态冰，在一定的真空度下，固态水可以直接变为水蒸气，在真空环境中相应的温度和气压下对物料进行冷冻干燥处理。通过真空冷冻干燥得到的产品具有复水性好、稳定性高以及方便运输、能长期保存等优点，因此，广泛应用于真空冷冻干燥技术广泛应用于医药、食品、材料等领域，但由于该技术存在干燥时间长、能耗大等缺点，因而其推广受到限制。为了减少干燥时间，降低操作成本，得到高质量的冻干产品，必须对最佳冻干工艺条件进行研究^[1-5]，但所需试验量大、成本高，且结果对整个过程的优化缺乏预测性。从机理方面，冻干过程中升华干燥阶段耗时最长，而预冻阶段影响冻结体的结晶情况即干燥阶段的孔隙结构，对冻干过程产生决定性影响。NAKAGAWA等^[6-8]通过模拟得到了冻干速率和成核温度对冰晶尺寸等参数的影响；贾颖姣等^[9]针对液态物料提出了3种冻结方式，并得出了不同冻结方式对预冻溶液冰晶粒度分布规律的影响。为了提高冻干效率，需要定量地得出不同冻结方式下形成冰晶粒度的分布情况，以及不同冰晶粒度对冻干速率、能耗和复水性等参数的影响，从而获得冻结条件与冻干效率和产品质量的对应情况，再通过逆向推理从机理上求得最佳冻干曲线。为此，本文作者基于数值模拟针对不同冻结方式与结晶粒度平均值和均匀度之间的函数关系进行研究，以便根据冻结方式预测其粒度分布，并根据干燥阶段的需要来调控预冻条件，为研究冻干过程中预冻阶段和干燥过程之间的关系提供依据。

1 数值模拟

1.1　物理模型

基于冻干箱内搁板冻结的预冻方式，以冻干机中2块搁板间1个西林瓶内溶液和它所分配的周围空间作为研究对象，忽略西林瓶液面以上的瓶身和瓶口部分对换热造成的影响。简化后的物理模型如图1所示。

图1　西林瓶装物料预冻模型示意图

Fig. 1　Schematic diagram of freezing material model in vial

1.2　数学模型

在溶液的冻结过程中，固液相变界面随着时间的变化而移动，称为三维两相stefan问题^[10-13]。针对该物理模型进行如下假设：1) 只发生热传递，忽略质传递，且传热方式为导热，忽略热辐射和热对流；2) 忽略固相和液相的密度差异性，即冷冻前后溶液体积不变；3) 冻结过程只发生在固液相界面上，其温度在固化温度和液化温度之间，瓶内冰晶为异质成核，不存在过冷情况；4) 冻结区域和液相区域各自分布均匀且内部各处物性分别一致，忽略包含的不凝性气体；5) 预冻后只有晶态，没有形成玻璃态物质。针对相变潜热的问题，采取多孔-焓法^[14-16]来解决，其中，液相的体积分数由热焓平衡法^[17]通过计算得出。能量方程如下：

(1)

式中：为材料密度，kg/m³；k为导热系数，W/(m²K)；T为物料温度，K；H为材料的焓，J。H由显焓h和潜热ΔH之和所得：

H=h+ΔH                                  (2)

(3)

其中：h_ref为参考焓，J/kg；T_ref为参考温度，K；比定压热容c_p定义为

(4)

式中：c_p-s和c_p-l分别为固相和液相的比定压热容，J/(kg·K)；T_solidus和T_liquidus分别为固化温度和液化温度，K，糊化区域的温度介于这两者之间。

已释放的潜热ΔH为溶液单位潜热L的函数：

ΔH=βL                                      (5)

其中：ΔH在0~L之间变化；β为液相质量分数，

(6)

边界条件为：

(7)

(8)

式中：k_g为空气导热系数，W/(m²K)；ρ_w为西林瓶壁材料的密度，kg/m³；c_p-w为西林瓶壁的比定压热容，J/(kgK)；T_w为西林瓶壁的温度，K；T_g为空气侧的温度，K。

1.3　数值解法

运用Ansys ICEM软件建模并划分网格，且进行网格无关性验证^[9]。导入Ansys Fluent 15.0 软件，选取凝固和熔化(融化)模型，采用有限体积法的二阶迎风格式对控制方程进行空间离散化，残差设置为10^-6，运用SIMPLE算法对模型温度进行求解。将所得温度场变化的数据导入Matlab软件计算冰晶粒度，并针对该方法进行模型验证^[9]。粒晶颗粒粒度计算公式^[6]为

(9)

式中：D为冰晶颗粒粒度，μm；c为常数；R为相界面的移动速率，mm/s；G为冻结区域的温度梯度，K/mm；R和G为通过处理数值计算结果即温度场间接得到，详细计算过程见文献[9]。

关于冰晶粒度的分布，本文选择平均值和标准差2个参数进行研究，计算式如下：

(10)

(11)

式中：D^*为西林瓶中冰晶粒度平均值，μm；n为参考冰晶粒度的位置个数；s为冰晶粒度分布标准差，μm。

1.4　计算参数

1.4.1　几何模型　

甘露醇学名为己六醇，常被用作冻干药品的低温保护剂和赋形剂^[5,18-19]。本文选取体积分数为10%甘露醇为研究对象，其对应图1中的几何尺寸如表1所示。

表1　几何模型尺寸

Table 1　Geometric modelmm

注：参数见图1。

1.4.2　模拟工况　

选取搁板温度作为预冻过程中的控制变量，分别采用以下3组不同的冻结方式：1) 将常温物料放入温度恒定的搁板上进行预冻(以下简称恒温条件)；2) 先将物料放入常温冻干箱内，然后控制搁板温度以一定速率下降直至冻结(以下简称恒温降条件)；3) 先使搁板温度由常温按照一定的速度下降至养晶温度 T_m并维持一段时间，而后继续以同样的速度降温(以下简称混养晶条件)。具体模拟工况如表2所示，这3种冻结方式都在常压环境下进行。

表2　3种不同冻结方式下的模拟工况

Table 2　Simulation conditions of three different freezing modes

2 关系式回归研究

2.1　恒温预冻方式

在恒温冻结方式下，模拟得出各设定温度下冰晶粒度的平均值和标准差，结果如图2所示。为了明确搁板温度与平均粒度和标准差之间的关系，进行多项式回归拟合，拟合的多项式为

(12)

式中：A₀，A₁，…，A_p分别为该函数中的待定系数；p为待定系数个数；为拟合得出的冰晶粒度分布参数，即平均值或标准差；x为自变量，即搁板温度。

图2　恒温预冻方式拟合曲线

Fig. 2　Fitting curves in constant temperature condition

令

(13)

式中：y_i为真实的冰晶粒度分布参数；δ_i为残量，表示拟合的误差。按照最小二乘法准则^[20]确定系数A_i，即求得使残量平方和最小的A_i。平方和S为

=

(14)

式中：S为残量平方和。

联立并求解式(10)和(12)便可求得拟合多项式。分析对比2次、3次和4次多项式的不确定度，选择最佳拟合多项式，可得搁板温度-平均粒度的拟合多项式为

(15)

式中：为搁板恒温冻结方式下拟合得到的冰晶粒度平均值，μm；T_b为搁板恒定温度，K。

所得搁板温度-标准差多项式为

(16)

式中：s₁为搁板恒温冻结方式下得到的冰晶粒度分布标准差，μm。搁板温度-平均粒度拟合曲线与搁板温度-标准差拟合曲线如图2所示。

2.2　恒温降预冻方式

在恒温降条件下，模拟各设定温降速率条件下冰晶粒度的平均值和标准差，结果如图3所示。根据同样步骤进行多项式拟合，通过不确定度选择最佳拟合多项式，得温降速率-平均粒度的拟合多项式为

(17)

式中：为恒温降冻结方式下拟合得到的冰晶粒度平均值，μm；v为搁板降温速率，K/s。所得的温降速率-标准差拟合多项式为

(18)

式中：s₂为恒温降预冻结方式下拟合得到的冰晶粒度分布标准差，μm。

图3　恒温降预冻方式拟合曲线

Fig. 3　Fitting curves in constant temperature condition drop

2.3　混养晶预冻方式

在混养晶条件下，模拟出各设定养晶温度条件下冰晶粒度的平均值和标准差，结果如图4所示，可得45 min为粒度分布最均匀的养晶时间。

图4　混养晶预冻方式拟合曲线

Fig. 4　Fitting curves in crystal growing condition

采用同样步骤进行多项式拟合，通过不确定度选择最佳拟合多项式，得养晶时间-平均粒度的拟合多项式为

(19)

式中：为混养晶预冻方式下拟合得到的冰晶粒度平均值，μm；t为养晶时间，s。所得的养晶时间-标准差拟合多项式为

+

(20)

式中：s₃为混养晶预冻方式下拟合得到的冰晶粒度分布标准差，μm。

3 模型检验和不确定度分析

3.1　多项式模型的有效检验

为了检验冰晶粒度和冻结条件之间确实存在多项式关系，需要进行显著性检验，根据多元统计理论，若取，则实际模型的一元多项式非线性回归问题可简化为多元线性回归问题。可根据多元线性回归的显著性检验统计量检验多项式回归方程的显著性。

假设H₀：A₁=0，A₂=0，…，A_j=0。

检验统计量：

  (21)

式中：F为检验统计量；为由拟合曲线计算的粒径，μm；y_i为由仿真得到的粒径，μm；为y_i的平均值；p为变量的自由度，它与拟合多项式的阶数有关，对于多种冻结方式，分别取3，2和4。

若F＞F_1-α(p, n-p-1)，则拒绝H₀，所得多项式分布在置信水平α上是显著的。所得结果如表3所示。

3.2　不确定度分析

从3种冻结方式的拟合曲线可以看出，拟合结果与原始结果比较接近，拟合效果较好。

3.2.1　均方根误差　

均方根误差用来衡量预测值与原始值之间的误差，在一定程度上反映了所得回归公式的效果。

残量平方和为

(22)

均方根误差为

(23)

式中：S_D1为残差平方和；δ_i为拟合误差，μm。

3.2.2　校正决定系数　

虽然均方根误差可以对拟合进行定量判断，但均方根误差也存在一定局限。为了获得最佳的拟合优度，参考校正决定系数，其值在0~1之间变化，当校正决定系数接近1时，表明拟合效果好。

决定系数R为

(24)

校正决定系数R_a为

(25)

3.2.3　回归参数标准偏差　

标准偏差是用来评定回归参数的准确度的量。标准偏差越小，则所得参数越准确。

残余标准偏差s_e为

(26)

各参数标准偏差s_e(A_i)(i=1, 2, 3, 4)为

(27)

其中：E=B^-1；B为求解待定系数A_i的系数矩阵。具体求解过程见文献[20]；E_ii为矩阵E中对应的数值。所得结果如表4所示。

表3　不同预冻条件下拟合多项式显著性检验

Table 3　Significance test of fitting polynomial under different pre-freezing conditions

表4　不同预冻条件下拟合多项式的不确定度

Table 4　Uncertainties of fitting polynomials under different pre-freezing conditions

4 结论

1) 在搁板温度恒定的冻结方式下，搁板温度-平均粒径与搁板温度-标准差之间的回归关系式均为三次项式；在降温速率恒定的冻结方式下，温降温度-平均粒径与温降速率回归关系式均为三次多项式。而在含养晶过程的恒温降过程中，养晶时间-平均粒径关系式为二次多项式，而养晶时间-标准差关系式为四次多项式。

2) 拟合优度和参数准确度比较高，而且多项式分布显著性水平在0.99以上。

3) 后续研究可根据结果推算到其他预冻条件下的冰晶粒度分布情况，定量研究预冻条件对冻干效率和产品质量的影响，优化冻干曲线。

参考文献：

[1] 赵延强. 具有初始孔隙多孔物料冷冻干燥的实验研究[D]. 大连: 大连理工大学化工机械与安全学院, 2015: 24-28.

ZHAO Yanqiang. Experimental study on freeze-drying of porous materials with initial porosity[D]. Dalian: Dalian University of Technology. School of Chemical Machinery and Safety Engineering, 2015: 24-28.

[2] 黄鸿兵, 徐丰莲, 周光宏. 冷冻贮藏过程中温度波动对猪肉肌间冰晶颜色和新鲜度的影响[J]. 食品科学, 2006, 27(8): 49-52.

HUANG Hongbing, XU Fenglian, ZHOU Guanghong. Effects of temperature fluctuation on ice crystals color and freshness in pork muscle during freezing and storage[J]. Food Science, 2006, 27(8): 49-52.

[3] ARGYROPOULOS D, HEINDL A, MULLER J. Assessment of convection, hot-air combined with microwave-vacuum and freeze-drying methods for mushrooms with regard to product quality[J]. International Journal of Food Science & Technology, 2011, 46(2): 333-342.

[4] LIN Liulan, WANG Zhikun, ZHOU Liping, et al. The influence of prefreezing temperature on pore structure in freeze-dried β-TCP scaffolds[J]. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 2013, 227(1): 50-57.

[5] GROHGANZ H, LEE Y Y, RANTANEN J, et al. The influence of lysozyme on mannitol polymorphism in freeze-dried and spray-dried formulations depends on the selection of the drying process[J]. International Journal of Pharmaceutics, 2013, 447(1/2): 224-230.

[6] NAKAGAWA K, HOTTOT A, VESSOT S, et al. Modeling of freezing step during freeze-drying of drugs in vials[J]. AIChE Journal, 2007, 53(5): 1362-1372.

[7] NAKAGAWA K, HOTTOT A, VESSOT S, et al. Modeling of freezing step during vial freeze-drying of pharmaceuticals: influence of nucleation temperature on primary drying rate[J]. Asia-Pacific Journal of Chemical Engineering, 2011, 6(2): 288-293.

[8] NAKAGAWA K, OCHIAI T. A mathematical model of multi-dimensional freeze-drying for food products[J]. Journal of Food Engineering, 2015, 161: 55-67.

[9] 贾颖姣, 刘志强. 冻结方式对预冻甘露醇冰晶粒度分布的影响[J]. 中南大学学报(自然科学版), 2019, 50(7): 1704-1711.

JIA Yingjiao, LIU Zhiqiang. Effects of freezing methods on size distribution of ice crystals during pre-freezing mannitol[J]. Journal of Central South University(Science and Technology), 2019, 50(7): 1704-1711 .

[10] 王新房, 汪春华, 吴光超. 基于有限元的一类 Stefan 问题数值模型研究[J]. 自动化技术与应用, 2007, 26(5): 9-12.

WANG Xinfang, WANG Chunhua, WU Guangchao. Numerical model of a class of Stefan problem based on finite element method[J]. Techniques of Automation and Applications, 2007, 26(5): 9-12.

[11] BARRETT J W, GARCKE H, NURNBERG R. On stable parametric finite element methods for the Stefan problem and the Mullins–Sekerka problem with applications to dendritic growth[J]. Journal of Computational Physics, 2010, 229(18): 6270-6299.

[12] GUPTA S C. The classical Stefan problem: basic concepts, modelling and analysis with quasi-analytical solutions and methods[J]. Contemporary Physics, 2018, 59(4): 428-429.

[13] SLOTA D. The application of the homotopy perturbation method to one-phase inverse Stefan problem[J]. International Communications in Heat and Mass Transfer, 2010, 37(6): 587-592.

[14] ZHAO Weihuan, ELMOZUGHI A F, OZTEKIN A, et al. Heat transfer analysis of encapsulated phase change material for thermal energy storage[J]. International Journal of Heat and Mass Transfer, 2013, 63: 323-335.

[15] ZHANG Y, DU K, HE J P, et al. Impact factors analysis of the enthalpy method and the effective heat capacity method on the transient nonlinear heat transfer in phase change materials (PCMs)[J]. Numerical Heat Transfer, Part A: Applications, 2014, 65(1): 66-83.

[16] BRENT A D, VOLLER V R, REID K T J. Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal[J]. Numerical Heat Transfer, Part A: Applications, 1988, 13(3): 297-318.

[17] KOTAS T J. The exergy method of thermal plant analysis[M]. Stoneham Ma: Butterworth Publishers, 1985: 99-161.

[18] KAIALY W, NOKHODCHI A. Freeze-dried mannitol for superior pulmonary drug delivery via dry powder inhaler[J]. Pharmaceutical Research, 2013, 30(2): 458-477.

[19] MEHTA M, BHARDWAJ S P, SURYANARAYANAN R. Controlling the physical form of mannitol in freeze-dried systems[J]. European Journal of Pharmaceutics and Biopharmaceutics, 2013, 85(2): 207-213.

[20] 李庆扬. 数值分析[M]. 北京: 清华大学出版社有限公司, 2001: 73-76.

LI Qingyang. Numerical analysis[M]. Beijing: Tsinghua University Press Ltd., 2001: 73-76.

(编辑  陈灿华)

收稿日期： 2018 -09 -10; 修回日期： 2018 -11 -12

基金项目(Foundation item)：国家自然科学基金资助项目(51676209)；湖南省科技计划重点研发项目(2016NK2139); 中南大学研究生自由探索创新项目(2018zzts489)(Project(51676209) supported by the National Natural Science Foundation of China; Project(2016NK2139) supported by the Key Science and Technology Plan of Hunan Province; Project(2018zzts489) supported by Independent Innovation and Exploration Programs for Graduate Students of Central South University);

通信作者：刘志强，博士，教授，从事低温工程传热传质与系统节能研究；E-mail: liuzq@csu.edu.cn