溶液預(yù)凍方式與冰晶粒度分布之間的關(guān)系

賈穎姣，余云霞，劉志強，羅春

(1.中南大學(xué)能源科學(xué)與工程學(xué)院，湖南長沙，410083；2.新奧泛能網(wǎng)絡(luò)科技有限公司，北京，100176)

真空冷凍干燥是一種穩(wěn)定的干燥過程。其原理如下：通過低溫冷凍使液態(tài)水凍結(jié)為固態(tài)冰，在一定的真空度下，固態(tài)水可以直接變?yōu)樗魵猓谡婵窄h(huán)境中相應(yīng)的溫度和氣壓下對物料進(jìn)行冷凍干燥處理。通過真空冷凍干燥得到的產(chǎn)品具有復(fù)水性好、穩(wěn)定性高以及方便運輸、能長期保存等優(yōu)點，因此，廣泛應(yīng)用于真空冷凍干燥技術(shù)廣泛應(yīng)用于醫(yī)藥、食品、材料等領(lǐng)域，但由于該技術(shù)存在干燥時間長、能耗大等缺點，因而其推廣受到限制。為了減少干燥時間，降低操作成本，得到高質(zhì)量的凍干產(chǎn)品，必須對最佳凍干工藝條件進(jìn)行研究[1-5]，但所需試驗量大、成本高，且結(jié)果對整個過程的優(yōu)化缺乏預(yù)測性。從機(jī)理方面，凍干過程中升華干燥階段耗時最長，而預(yù)凍階段影響凍結(jié)體的結(jié)晶情況即干燥階段的孔隙結(jié)構(gòu)，對凍干過程產(chǎn)生決定性影響。NAKAGAWA 等[6-8]通過模擬得到了凍干速率和成核溫度對冰晶尺寸等參數(shù)的影響；賈穎姣等[9]針對液態(tài)物料提出了3 種凍結(jié)方式，并得出了不同凍結(jié)方式對預(yù)凍溶液冰晶粒度分布規(guī)律的影響。為了提高凍干效率，需要定量地得出不同凍結(jié)方式下形成冰晶粒度的分布情況，以及不同冰晶粒度對凍干速率、能耗和復(fù)水性等參數(shù)的影響，從而獲得凍結(jié)條件與凍干效率和產(chǎn)品質(zhì)量的對應(yīng)情況，再通過逆向推理從機(jī)理上求得最佳凍干曲線。為此，本文作者基于數(shù)值模擬針對不同凍結(jié)方式與結(jié)晶粒度平均值和均勻度之間的函數(shù)關(guān)系進(jìn)行研究，以便根據(jù)凍結(jié)方式預(yù)測其粒度分布，并根據(jù)干燥階段的需要來調(diào)控預(yù)凍條件，為研究凍干過程中預(yù)凍階段和干燥過程之間的關(guān)系提供依據(jù)。

1 數(shù)值模擬

1.1 物理模型

基于凍干箱內(nèi)擱板凍結(jié)的預(yù)凍方式，以凍干機(jī)中2 塊擱板間1 個西林瓶內(nèi)溶液和它所分配的周圍空間作為研究對象，忽略西林瓶液面以上的瓶身和瓶口部分對換熱造成的影響。簡化后的物理模型如圖1所示。

1.2 數(shù)學(xué)模型

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

在溶液的凍結(jié)過程中，固液相變界面隨著時間的變化而移動，稱為三維兩相stefan 問題[10-13]。針對該物理模型進(jìn)行如下假設(shè)：1)只發(fā)生熱傳遞，忽略質(zhì)傳遞，且傳熱方式為導(dǎo)熱，忽略熱輻射和熱對流；2)忽略固相和液相的密度差異性，即冷凍前后溶液體積不變；3)凍結(jié)過程只發(fā)生在固液相界面上，其溫度在固化溫度和液化溫度之間，瓶內(nèi)冰晶為異質(zhì)成核，不存在過冷情況；4)凍結(jié)區(qū)域和液相區(qū)域各自分布均勻且內(nèi)部各處物性分別一致，忽略包含的不凝性氣體；5)預(yù)凍后只有晶態(tài)，沒有形成玻璃態(tài)物質(zhì)。針對相變潛熱的問題，采取多孔-焓法[14-16]來解決，其中，液相的體積分?jǐn)?shù)由熱焓平衡法[17]通過計算得出。能量方程如下：

式中：ρ為材料密度，kg/m3；k為導(dǎo)熱系數(shù)，W/(m2?K)；T為物料溫度，K；H為材料的焓，J。H由顯焓h和潛熱ΔH之和所得：

其中：href為參考焓，J/kg；Tref為參考溫度，K；比定壓熱容cp定義為

式中：cp-s和cp-l分別為固相和液相的比定壓熱容，J/(kg?K)；Tsolidus和Tliquidus分別為固化溫度和液化溫度，K，糊化區(qū)域的溫度介于這兩者之間。

已釋放的潛熱ΔH為溶液單位潛熱L的函數(shù)：

其中：ΔH在0～L之間變化；β為液相質(zhì)量分?jǐn)?shù)，

邊界條件為：

式中：kg為空氣導(dǎo)熱系數(shù)，W/(m2?K)；ρw為西林瓶壁材料的密度，kg/m3；cp-w為西林瓶壁的比定壓熱容，J/(kg?K)；Tw為西林瓶壁的溫度，K；Tg為空氣側(cè)的溫度，K。

1.3 數(shù)值解法

運用Ansys ICEM 軟件建模并劃分網(wǎng)格，且進(jìn)行網(wǎng)格無關(guān)性驗證[9]。導(dǎo)入Ansys Fluent 15.0 軟件，選取凝固和熔化(融化)模型，采用有限體積法的二階迎風(fēng)格式對控制方程進(jìn)行空間離散化，殘差設(shè)置為10-6，運用SIMPLE 算法對模型溫度進(jìn)行求解。將所得溫度場變化的數(shù)據(jù)導(dǎo)入Matlab軟件計算冰晶粒度，并針對該方法進(jìn)行模型驗證[9]。粒晶顆粒粒度計算公式[6]為

式中：D為冰晶顆粒粒度，μm；c為常數(shù)；R為相界面的移動速率，mm/s；G為凍結(jié)區(qū)域的溫度梯度，K/mm；R和G為通過處理數(shù)值計算結(jié)果即溫度場間接得到，詳細(xì)計算過程見文獻(xiàn)[9]。

關(guān)于冰晶粒度的分布，本文選擇平均值和標(biāo)準(zhǔn)差2個參數(shù)進(jìn)行研究，計算式如下：

式中：D*為西林瓶中冰晶粒度平均值，μm；n為參考冰晶粒度的位置個數(shù)；s為冰晶粒度分布標(biāo)準(zhǔn)差，μm。

1.4 計算參數(shù)

1.4.1 幾何模型

甘露醇學(xué)名為己六醇，常被用作凍干藥品的低溫保護(hù)劑和賦形劑[5,18-19]。本文選取體積分?jǐn)?shù)為10%甘露醇為研究對象，其對應(yīng)圖1中的幾何尺寸如表1所示。

表1 幾何模型尺寸Table1 Geometric model mm

1.4.2 模擬工況

選取擱板溫度作為預(yù)凍過程中的控制變量，分別采用以下3組不同的凍結(jié)方式：1)將常溫物料放入溫度恒定的擱板上進(jìn)行預(yù)凍(以下簡稱恒溫條件)；2)先將物料放入常溫凍干箱內(nèi)，然后控制擱板溫度以一定速率下降直至凍結(jié)(以下簡稱恒溫降條件)；3)先使擱板溫度由常溫按照一定的速度下降至養(yǎng)晶溫度Tm并維持一段時間，而后繼續(xù)以同樣的速度降溫(以下簡稱混養(yǎng)晶條件)。具體模擬工況如表2所示，這3種凍結(jié)方式都在常壓環(huán)境下進(jìn)行。

表2 3種不同凍結(jié)方式下的模擬工況Table2 Simulation conditions of three different freezing modes

2 關(guān)系式回歸研究

2.1 恒溫預(yù)凍方式

在恒溫凍結(jié)方式下，模擬得出各設(shè)定溫度下冰晶粒度的平均值和標(biāo)準(zhǔn)差，結(jié)果如圖2所示。為了明確擱板溫度與平均粒度和標(biāo)準(zhǔn)差之間的關(guān)系，進(jìn)行多項式回歸擬合，擬合的多項式為

式中：A0，A1，…，Ap分別為該函數(shù)中的待定系數(shù)；p為待定系數(shù)個數(shù)；?(x)為擬合得出的冰晶粒度分布參數(shù)，即平均值或標(biāo)準(zhǔn)差；x為自變量，即擱板溫度。

圖2 恒溫預(yù)凍方式擬合曲線Fig.2 Fitting curves in constant temperature condition

令

式中：yi為真實的冰晶粒度分布參數(shù)；δi為殘量，表示擬合的誤差。按照最小二乘法準(zhǔn)則[20]確定系數(shù)Ai，即求得使殘量?i平方和最小的Ai。平方和S為

式中：S為殘量平方和。

聯(lián)立并求解式(10)和(12)便可求得擬合多項式。分析對比2 次、3 次和4 次多項式的不確定度，選擇最佳擬合多項式，可得擱板溫度-平均粒度的擬合多項式為

式中：D*1為擱板恒溫凍結(jié)方式下擬合得到的冰晶粒度平均值，μm；Tb為擱板恒定溫度，K。

所得擱板溫度-標(biāo)準(zhǔn)差多項式為

式中：s1為擱板恒溫凍結(jié)方式下得到的冰晶粒度分布標(biāo)準(zhǔn)差，μm。擱板溫度-平均粒度擬合曲線與擱板溫度-標(biāo)準(zhǔn)差擬合曲線如圖2所示。

2.2 恒溫降預(yù)凍方式

在恒溫降條件下，模擬各設(shè)定溫降速率條件下冰晶粒度的平均值和標(biāo)準(zhǔn)差，結(jié)果如圖3所示。根據(jù)同樣步驟進(jìn)行多項式擬合，通過不確定度選擇最佳擬合多項式，得溫降速率-平均粒度的擬合多項式為

式中：D*2為恒溫降凍結(jié)方式下擬合得到的冰晶粒度平均值，μm；v為擱板降溫速率，K/s。所得的溫降速率-標(biāo)準(zhǔn)差擬合多項式為

式中：s2為恒溫降預(yù)凍結(jié)方式下擬合得到的冰晶粒度分布標(biāo)準(zhǔn)差，μm。

圖3 恒溫降預(yù)凍方式擬合曲線Fig.3 Fitting curves in constant temperature condition drop

2.3 混養(yǎng)晶預(yù)凍方式

在混養(yǎng)晶條件下，模擬出各設(shè)定養(yǎng)晶溫度條件下冰晶粒度的平均值和標(biāo)準(zhǔn)差，結(jié)果如圖4所示，可得45 min為粒度分布最均勻的養(yǎng)晶時間。

采用同樣步驟進(jìn)行多項式擬合，通過不確定度選擇最佳擬合多項式，得養(yǎng)晶時間-平均粒度的擬合多項式為

式中：D*3為混養(yǎng)晶預(yù)凍方式下擬合得到的冰晶粒度平均值，μm；t為養(yǎng)晶時間，s。所得的養(yǎng)晶時間-標(biāo)準(zhǔn)差擬合多項式為

式中：s3為混養(yǎng)晶預(yù)凍方式下擬合得到的冰晶粒度分布標(biāo)準(zhǔn)差，μm。

圖4 混養(yǎng)晶預(yù)凍方式擬合曲線Fig.4 Fitting curves in crystal growing condition

3 模型檢驗和不確定度分析

3.1 多項式模型的有效檢驗

為了檢驗冰晶粒度和凍結(jié)條件之間確實存在多項式關(guān)系，需要進(jìn)行顯著性檢驗，根據(jù)多元統(tǒng)計理論，若取Xl=xl，則實際模型的一元多項式非線性回歸問題可簡化為多元線性回歸問題。可根據(jù)多元線性回歸的顯著性檢驗統(tǒng)計量檢驗多項式回歸方程的顯著性。

假設(shè)H0：A1=0，A2=0，…，Aj=0。

檢驗統(tǒng)計量：

式中：F為檢驗統(tǒng)計量；為由擬合曲線計算的粒徑，μm；yi為由仿真得到的粒徑，μm；為yi的平均值；p為變量的自由度，它與擬合多項式的階數(shù)有關(guān)，對于多種凍結(jié)方式，分別取3，2和4。

若F＞F1-α(p,n-p-1)，則拒絕H0，所得多項式分布在置信水平α上是顯著的。所得結(jié)果如表3所示。

3.2 不確定度分析

從3種凍結(jié)方式的擬合曲線可以看出，擬合結(jié)果與原始結(jié)果比較接近，擬合效果較好。

3.2.1 均方根誤差

均方根誤差用來衡量預(yù)測值與原始值之間的誤差，在一定程度上反映了所得回歸公式的效果。殘量平方和為

均方根誤差為

式中：SD1為殘差平方和；δi為擬合誤差，μm。

3.2.2 校正決定系數(shù)

雖然均方根誤差可以對擬合進(jìn)行定量判斷，但均方根誤差也存在一定局限。為了獲得最佳的擬合優(yōu)度，參考校正決定系數(shù)，其值在0～1之間變化，當(dāng)校正決定系數(shù)接近1時，表明擬合效果好。

決定系數(shù)R為

校正決定系數(shù)Ra為

3.2.3 回歸參數(shù)標(biāo)準(zhǔn)偏差

標(biāo)準(zhǔn)偏差是用來評定回歸參數(shù)的準(zhǔn)確度的量。標(biāo)準(zhǔn)偏差越小，則所得參數(shù)越準(zhǔn)確。

殘余標(biāo)準(zhǔn)偏差se為

各參數(shù)標(biāo)準(zhǔn)偏差se(Ai)(i=1,2,3,4)為

其中：E=B-1；B為求解待定系數(shù)Ai的系數(shù)矩陣。具體求解過程見文獻(xiàn)[20]；Eii為矩陣E中對應(yīng)的數(shù)值。所得結(jié)果如表4所示。

表3 不同預(yù)凍條件下擬合多項式顯著性檢驗Table3 Significance test of fitting polynomial under different pre-freezing conditions

表4 不同預(yù)凍條件下擬合多項式的不確定度Table4 Uncertainties of fitting polynomials under different pre-freezing conditions

4 結(jié)論

1)在擱板溫度恒定的凍結(jié)方式下，擱板溫度-平均粒徑與擱板溫度-標(biāo)準(zhǔn)差之間的回歸關(guān)系式均為三次項式；在降溫速率恒定的凍結(jié)方式下，溫降溫度-平均粒徑與溫降速率回歸關(guān)系式均為三次多項式。而在含養(yǎng)晶過程的恒溫降過程中，養(yǎng)晶時間-平均粒徑關(guān)系式為二次多項式，而養(yǎng)晶時間-標(biāo)準(zhǔn)差關(guān)系式為四次多項式。

2)擬合優(yōu)度和參數(shù)準(zhǔn)確度比較高，而且多項式分布顯著性水平在0.99以上。

3)后續(xù)研究可根據(jù)結(jié)果推算到其他預(yù)凍條件下的冰晶粒度分布情況，定量研究預(yù)凍條件對凍干效率和產(chǎn)品質(zhì)量的影響，優(yōu)化凍干曲線。