賈穎姣，劉志強(qiáng)，羅春

(中南大學(xué)能源科學(xué)與工程學(xué)院，湖南長沙，410083)

真空冷凍干燥技術(shù)是目前最優(yōu)良的干燥技術(shù)之一，凍干產(chǎn)品具有復(fù)水性好、穩(wěn)定性高以及方便運(yùn)輸、長期保存等優(yōu)點，因而在醫(yī)藥、食品、材料等領(lǐng)域得到廣泛應(yīng)用。但冷凍干燥耗時長、能耗大、設(shè)備初投資大，這些缺點制約了冷凍干燥技術(shù)的工業(yè)化發(fā)展。因此，優(yōu)化凍干過程的工藝參數(shù)以提高凍干效率很有必要。真空冷凍干燥過程主要分為預(yù)凍階段、升華干燥階段和解析干燥階段。其中，升華干燥階段耗時最長，而預(yù)凍階段的凍結(jié)方式是通過改變物料內(nèi)部冰晶的尺寸來決定干燥層內(nèi)孔隙的尺寸，對后期的干燥階段尤其是升華干燥階段產(chǎn)生影響，且冰晶分布的均勻度與干燥過程出現(xiàn)的空洞萎縮等現(xiàn)象有關(guān)。因此，人們對如何減少干燥時間、降低操作成本、得到高質(zhì)量的凍干產(chǎn)品進(jìn)行了研究。徐立沖等[1]研究了真空冷凍條件下停留時間、過冷度等對人工海水中冰晶生長速度的影響及作用機(jī)理；黃鴻兵等[2]研究了凍藏過程中-10～-20 ℃之間溫度波動對豬肉肌間冰晶、顏色和新鮮度的影響，得出溫度波動越劇烈，冰晶體積增大幅度越大，且-10 ℃恒溫凍藏樣品具有最小的直徑；關(guān)志強(qiáng)等[3]對文蛤和波紋巴非蛤肉在-196，-78和-18 ℃中的凍結(jié)過程進(jìn)行了研究，發(fā)現(xiàn)-78 ℃凍結(jié)使組織結(jié)構(gòu)產(chǎn)生的冰晶較均勻，變形程度次之；駱麗君等[4]研究了不同凍結(jié)溫度下面條內(nèi)部水分分布、微觀結(jié)構(gòu)、冰晶形態(tài)等的變化，發(fā)現(xiàn)凍結(jié)溫度越低，凍結(jié)速度越快，凍結(jié)時間和通過最大冰晶形成區(qū)域時間越少，面條內(nèi)部冰晶體積越小，分布越均勻；LIN等[5]通過實驗發(fā)現(xiàn)了使凍干腳手架的空隙分布與網(wǎng)狀結(jié)構(gòu)都好的最佳預(yù)凍溫度；NAKAGAWA 等[6-8]通過模擬和實驗證明了凍結(jié)過程中所建模型的正確性，并模擬出凍干速率和成核溫度對冰晶尺寸以及干燥層滲透率的影響。上述研究大部分都是基于實驗進(jìn)行的，為了更深入地研究冷凍干燥技術(shù)，本文作者基于數(shù)值模擬進(jìn)行系統(tǒng)研究，使生產(chǎn)者能根據(jù)具體需要制定最佳的凍結(jié)方式。

1 物理模型和控制方程

1.1 物理模型

目前，預(yù)凍方法有很多種，如凍干機(jī)內(nèi)擱板凍結(jié)、冷庫凍結(jié)、旋轉(zhuǎn)凍結(jié)器凍結(jié)等，而凍干液態(tài)原料的生產(chǎn)用醫(yī)藥凍干機(jī)多采用干燥箱內(nèi)擱板凍結(jié)[9]，且液態(tài)藥品通常用西林瓶裝排列置于擱板之上。因此，本文選取凍干機(jī)中兩擱板間的1個西林瓶和它所分配的周圍空間作為研究對象，具體物理模型及尺寸如圖1所示。常溫的藥品溶液被裝入西林瓶中放在擱板上，由上下擱板輸送冷量使其凍結(jié)。在建模過程中進(jìn)行如下假設(shè)：1)凍結(jié)過程只發(fā)生熱傳遞，忽略質(zhì)傳遞；2)傳熱方式主要為導(dǎo)熱，忽略熱輻射和熱對流；3)西林瓶為圓柱形，且不考慮液面以上的瓶身和瓶口部分對換熱造成的影響；4)凍結(jié)過程只發(fā)生在固液相界面上，其溫度在固化溫度和液化溫度之間；5)凍結(jié)區(qū)域和液相區(qū)域各相物料分布且內(nèi)部各處物性分別一致，忽略包含的不凝性氣體；6)物料關(guān)于半徑r=0呈軸對稱，忽略物性的方向性差別；7)預(yù)凍后只有晶態(tài)，沒有形成玻璃態(tài)物質(zhì)；8)預(yù)凍階段為常壓；9)瓶內(nèi)冰晶為異質(zhì)成核，且溫度一旦到達(dá)液化溫度即開始結(jié)晶生長，即成核溫度為271.28 K；10)忽略固相和液相的密度差異性，即冷凍前后溶液體積不變。

圖1 西林瓶裝物料預(yù)凍模型示意圖Fig.1 Schematic diagram of pre-freezing material model in a vial

1.2 控制方程

預(yù)凍過程是一個相變界面隨著時間移動的非穩(wěn)態(tài)過程，即三維兩相Stefan 問題[10-11]。本文運(yùn)用多孔-焓法[12-14]處理相變潛熱的問題。固液模糊區(qū)域被看作1個多孔區(qū)域，其中，某個小區(qū)域的孔隙率等于該處的液相體積分?jǐn)?shù)，而液相的體積分?jǐn)?shù)是基于熱焓平衡法[15-16]計算的。能量方程如下：

式中：ρ為材料密度，kg/m3；H為焓，J/kg；u為流體的流速，m/s；k為導(dǎo)熱系數(shù)，W/(m·K)；T為物料溫度，K；φ為熱源，W/m3。本文中u和φ都為0，液相導(dǎo)熱系數(shù)k1和固相導(dǎo)熱系數(shù)ks不同，材料的焓H為顯焓h與潛熱ΔH之和：

式中：href為參考焓值，J/kg；Tref為參考溫度，K；cp為恒壓比熱容，J/(kg·K)。

式中：cp-s和cp-l分別為固相和液相的比熱容，J/(kg·K)；Tsolidus和Tliquidus分別為固化溫度和液化溫度，K。糊化區(qū)域介于兩者之間，某處已釋放的潛熱ΔH可按下式計算：

式中：L為潛熱，J/kg。β為溶液的質(zhì)量分?jǐn)?shù)。

顯然，ΔH在0～L之間變化。

耦合壁面邊界條件如下：

1)當(dāng)z=18 mm時，

2)當(dāng)r=11 mm時，

式中：kg為空氣導(dǎo)熱系數(shù)，W/(m·K)；ρw為西林瓶壁面的密度，kg/m3；cp-w為西林瓶壁面的比熱容，J/(kg·K)；Tw為西林瓶壁的溫度，K；Tg為空氣的溫度，K；。

2 建模對象和計算方法

2.1 研究對象和內(nèi)容

甘露醇作為一種保護(hù)劑，常被應(yīng)用于凍干過程中[17-19]，因此，本文以體積分?jǐn)?shù)為10%甘露醇作為研究對象，其用于模型計算的相關(guān)參數(shù)如表1所示。

為了研究不同凍結(jié)方式對冷凍之后物料內(nèi)部冰晶顆粒粒度的影響，以生產(chǎn)中方便控制的擱板溫度作為預(yù)凍過程工藝曲線的控制變量，設(shè)計3組不同的凍結(jié)方式：1)把物料放在擱板上冷卻，擱板溫度保持恒定(-10，-15，-20 和-25 ℃)；2)物料放在擱板上冷卻，擱板從室溫開始冷卻且降溫速率保持恒定(-0.5，-1.0，-1.5 和-2 ℃/min)；3)先使擱板溫度由常溫按照-1 ℃/min 的速度下降至一定養(yǎng)晶溫度(0，-2，-4和-6 ℃)，并在該溫度維持15 min，然后以-1 ℃/min的速度繼續(xù)降溫。

表1 10%體積分?jǐn)?shù)甘露醇的物性參數(shù)Table 1 Physical properties mannitol with 10%volume fraction

2.2 結(jié)構(gòu)模型設(shè)置和網(wǎng)格無關(guān)性驗證

利用Ansys ICEM 軟件建立幾何模型。由于圖1中該模型關(guān)于xz平面和yz平面對稱，因此，截取整個模型的1/4劃分網(wǎng)格，并設(shè)置對稱邊界條件進(jìn)行計算。表2所示為4種計算網(wǎng)格數(shù)。在表2所示4種網(wǎng)格數(shù)量下，給出相同的邊界條件，分別計算300 s 內(nèi)甘露醇內(nèi)部的溫度變化，并選取中軸線(r=0)上z=1 mm這一點進(jìn)行研究，其溫度隨時間的變化關(guān)系如圖2所示。由圖2可看出：當(dāng)網(wǎng)格數(shù)為34.8萬時，溫度曲線出現(xiàn)明顯波動；當(dāng)網(wǎng)格數(shù)增加1倍時，波動趨勢減小但仍能看到有幾處凸點；當(dāng)網(wǎng)格數(shù)增大到70 萬時，曲線趨于平滑且網(wǎng)格數(shù)目對模擬結(jié)果的影響很小，可認(rèn)為此時已達(dá)到網(wǎng)格無關(guān)。因此，考慮到計算速度和模擬結(jié)果準(zhǔn)確性，本文后續(xù)的模擬工作均采用70 萬網(wǎng)格進(jìn)行計算，且都為六邊形結(jié)構(gòu)化網(wǎng)格，平均網(wǎng)格連長為0.2 mm，并按照需求在相變區(qū)域局部加密。

表2 4種計算網(wǎng)格數(shù)Table 2 Four kinds of grid numbers used in calculation 104個

2.3 計算方法和模型驗證

在以上模型的基礎(chǔ)上，基于Ansys Fluent 15.0 軟件，采用有限體積法的二階迎風(fēng)格式進(jìn)行控制方程的空間離散化，運(yùn)用迭代法對相變區(qū)域溫度進(jìn)行求解，其中，殘差設(shè)置為低于10-6，從而得出溫度場。

相變域內(nèi)溫度分布有一定的特點，取預(yù)凍條件為擱板溫度恒定-15 ℃進(jìn)行研究。120 s和2 625 s這2個時刻的3 條徑向線(z分別1，9 和17 mm)溫度分布如圖3所示。由圖3可得：西林瓶內(nèi)徑向溫差很小(最大溫差約為1.12 ℃)，因此，本文以縱向中心線上的點為對象研究每一層溫度以及粒度分布情況。

圖2 點(0,0,1)用4種網(wǎng)格數(shù)得出的溫度-時間圖Fig.2 Temperature-time diagrams of point at(0,0,1)with four kinds of grid numbers

在溶液的中軸線上(即r=0 且0≤z≤18 mm)，每隔2 mm 取1 點觀察其降溫及相變過程中溫度隨時間的變化，結(jié)果如圖4所示。得到溫度變化曲線后，利用每個點凍結(jié)完成的時間計算相界面通過該點時的移動速率R，由圖5中折線段的斜率表示。

相界面經(jīng)過上述觀測點時西林瓶內(nèi)的溫度場如圖6所示。從圖6可見：ti越大，則溫度分布越趨于均勻，對應(yīng)時刻的冷凍區(qū)域溫度梯度Gi由此得出。

一些研究者研究了冰晶顆粒的粒度D與相界面移動速率R和溫度梯度G之間的關(guān)系。例如，在金屬凝固過程中，KURZ等[20]得出D∝R-1G-1或者D∝R0.5；對于蘋果冷凍，BOMBEN 等[21]得出D∝R-0.5G-0.5；對于淀粉凝膠，REID[22]得出D∝R-x；對于明膠冷凍，WOINET 等[23]通過理論推導(dǎo)和實驗驗證得出D∝R-0.5G-0.5；對甘露醇冷凍，NAKAGAWA[6]運(yùn)用關(guān)系式D=cR-0.5G-0.5對冰晶粒度進(jìn)行計算，并用實驗證明了其適用性。本文采用如下關(guān)系式計算冰晶粒度：

圖3 溫度恒為-15 ℃時不同時間3條徑向線的溫度分布Fig.3 Temperature distribution of three radial lines at different time when the temperature is constantly-15 ℃

圖4 溫度恒定為-15 ℃時的溫度(T)-時間(t)關(guān)系曲線Fig.4 T-t diagrams when temperature is constantly-15 ℃

圖5 溫度恒定為-15 ℃的相變界面位置變化曲線Fig.5 Variation of position of freezing front versus time with shelf temperature of-15℃

圖6 不同時間西林瓶內(nèi)溫度分布Fig.6 Distribution of temperature in a vial at different time

其中：c為常數(shù)；R為相界面移動速率，mm/s；G為溫度梯度，K/mm。

為了驗證模型以及數(shù)值解的正確性，用該方法對文獻(xiàn)[6]中的同類型實驗進(jìn)行模擬。在計算區(qū)域中，各介質(zhì)的相關(guān)物性參數(shù)依據(jù)NAKAGAWA等[6]提供的參數(shù)設(shè)定，邊界條件設(shè)定如下：初始溫度為280 K，下擱板溫度下降控制為-0.72 K/min，上擱板和四周空氣壁面溫度為80 K 恒溫，溫度場數(shù)據(jù)在MATLAB 中進(jìn)行處理得到冰晶粒度分布，并對模擬的不同位置的粒度與其進(jìn)行比較。選取3 個點的冰晶粒度進(jìn)行對比，分別為A(0，0，2)mm，B(0，0，4)mm 和C(0，0，6)mm，結(jié)果如表3所示。

表3 模型驗證結(jié)果Table 3 Validation results of model

運(yùn)用該模型模擬得出粒度隨高度的分布趨勢與參考文獻(xiàn)[6]中的一致，且由表3可知，仿真計算得到的冰晶粒徑與參考值之間的誤差在允許范圍之內(nèi)，表明該仿真計算的結(jié)果具有可靠性。

3 模擬結(jié)果與分析

在第1 種凍結(jié)方式(即擱板溫度恒定)下，西林瓶中粒度分布如圖7所示。由圖7可得：粒度隨著高度增加呈現(xiàn)先增大后減小的變化趨勢，在z=16 mm附近達(dá)到峰值。不同溫度下的冰晶粒度分布見圖8。從圖8可見：擱板溫度Tb越低，則平均粒度D′ 越小，標(biāo)準(zhǔn)差S也越??；隨著溫度降低，平均粒度和標(biāo)準(zhǔn)差曲線都都趨于平緩，即要使西林瓶內(nèi)粒度小且分布均勻，則可以適當(dāng)調(diào)低擱板溫度。但若溫度過低，則會造成能源浪費(fèi)，效果也不明顯，且有可能使冷凍材料轉(zhuǎn)為玻璃態(tài)[24]，影響之后的干燥過程。

圖7 恒溫時西林瓶縱向冰晶粒度分布模擬Fig.7 Simulation of ice crystal size distribution in vial vertical direction with constant shelf temperature

圖8 恒溫時冰晶粒度平均值D′和標(biāo)準(zhǔn)差S與溫度的關(guān)系Fig.8 Mean value and standard deviation of ice crystal size distribution versus temperature

圖9 恒定溫降時西林瓶縱向冰晶粒度分布模擬Fig.9 Simulation of ice crystal size distribution in vial vertical direction with constant shelf temperature drop

在第2 種凍結(jié)方式(即擱板溫度以一定的速度下降)下，恒定溫降時西林瓶縱向冰晶尺寸分布模擬見圖9，恒溫時冰晶粒度平均值和標(biāo)準(zhǔn)差與降溫速率的關(guān)系見圖10。從圖9和圖10可見：西林瓶內(nèi)縱向粒度分布也呈先增大后減小趨勢，在z=16 mm附近達(dá)到峰值。對比擱板恒溫和擱板恒定降溫速率2種凍結(jié)方式，當(dāng)選定相同的目標(biāo)溫度時，第2種凍結(jié)方式比第1種凍結(jié)方式得出的平均粒度更大但分布更均勻。例如，當(dāng)降溫速率v為-0.5 ℃/min 時，溫度達(dá)到-25 ℃時物料已全部凍結(jié)即冰晶粒度分布已形成，此時，冰晶的平均粒度為130.25 μm，標(biāo)準(zhǔn)差為19.44 μm，而在恒溫條件下即溫度為-25 ℃時，冰晶的平均粒度為82.85 μm，標(biāo)準(zhǔn)差為38.19 μm。

圖10 恒溫時冰晶粒度平均值D′和標(biāo)準(zhǔn)差S與降溫速率δ的關(guān)系Fig.10 Mean values and standard deviation of crystal size distribution versus temperature drop

在第3 種凍結(jié)方式(即擱板溫度先由常溫以速度-1 ℃/min 下降到某一養(yǎng)晶溫度Tm并保持15 min，然后以相同的溫度下降)下，含不同養(yǎng)晶溫度時西林瓶縱向冰晶粒度分布模擬結(jié)果見圖11，含不同養(yǎng)晶溫度時冰晶粒度的平均值和標(biāo)準(zhǔn)差和與養(yǎng)晶溫度的關(guān)系見圖12。從圖11和圖12可知：當(dāng)養(yǎng)晶溫度為0 ℃和-2 ℃，這時的養(yǎng)晶溫度均高于甘露醇溶液的凝固溫度，西林瓶內(nèi)粒度分布曲線幾乎是重合的，2種養(yǎng)晶溫度下得到的平均粒度以及標(biāo)準(zhǔn)差也近似相等；當(dāng)養(yǎng)晶溫度為-4 ℃和-6 ℃，這2 個養(yǎng)晶溫度都低于甘露醇溶液的凝固溫度，最大冰晶粒度出現(xiàn)的位置并不相同。此外，其他位置的冰晶粒度分布比較均勻，這是因為養(yǎng)晶溫度變小。由此可見，養(yǎng)晶時間是調(diào)控冰晶粒度的1個重要因素。

圖11 含不同養(yǎng)晶溫度時西林瓶縱向冰晶粒度分布模擬Fig.11 Simulation of ice crystal size distribution in vial vertical direction with different crystal growing temperatures

圖12 含不同養(yǎng)晶溫度時冰晶粒度的平均值D和標(biāo)準(zhǔn)差S與養(yǎng)晶溫度的關(guān)系Fig.12 Mean values and standard deviation of ice crystal size distribution as functions of crystal growing temperature

與第2種凍結(jié)方式中降溫速度為-1 ℃/min的情況相比，含養(yǎng)晶階段的恒溫降(即第3 種凍結(jié)方式)會使得平均粒度變大，若養(yǎng)晶溫度高于凝固溫度，則幾乎不影響粒度分布的均勻度；要使均勻度發(fā)生變化，則可適當(dāng)調(diào)節(jié)養(yǎng)晶溫度，使其低于凝固溫度并控制養(yǎng)晶時間。

4 結(jié)論

1) 建立了凍干箱預(yù)凍西林瓶裝體積分?jǐn)?shù)為10%甘露醇溶液三維仿真模型，利用Ansys軟件模擬了該模型在預(yù)凍過程中的溫度場，并通過Matlab 計算得到凍結(jié)后冰晶粒度的分布。

2) 在擱板恒溫條件和勻速降溫條件下，粒度分布沿西林瓶高度方向都呈中間大、兩端小的變化趨勢，峰值出現(xiàn)在略低于物料與空氣界面處；隨著擱板溫度降低以及降溫速度加快，平均粒度減小，分布更加均勻，同時減小的趨勢逐漸變平緩；當(dāng)選定1個目標(biāo)溫度時，勻速降溫的方式比保持恒溫得出的冰晶粒度更大但更均勻；勻速降溫且當(dāng)溫度到達(dá)稍高于凝固點的溫度時，加入一定的養(yǎng)晶階段可增加平均粒度而不影響均勻度，要增大某些位置的粒度或調(diào)整均勻度，可視具體情況從降低養(yǎng)晶溫度以及調(diào)整養(yǎng)晶時間方面著手研究。

3) 本文研究可為工程實際中凍結(jié)方式的選取提供參考和指導(dǎo)，也為深入研究冰晶粒度對凍干過程中干燥過程的時間和凍干產(chǎn)品質(zhì)量的影響從而優(yōu)化凍干曲線提供參考。