王 瑾，李占勇，徐 慶，王瑞芳

(天津科技大學(xué)機(jī)械工程學(xué)院，天津 300222)

單液滴撞擊水平冷板面的模擬研究

王 瑾，李占勇，徐 慶，王瑞芳

(天津科技大學(xué)機(jī)械工程學(xué)院，天津 300222)

針對(duì)噴霧冷凍干燥中液滴在設(shè)備壁面的冷凍沉積問(wèn)題，在前期單液滴撞擊水平放置冷板面的實(shí)驗(yàn)基礎(chǔ)上，建立了單液滴與冷板面撞擊的模型，利用水平集方法與固化模型相結(jié)合，模擬了蒸餾水液滴(直徑為3.2mm，溫度為293K)分別從100、250、500mm高度撞擊不銹鋼板面(溫度為268K和253K)的過(guò)程．模擬的最大鋪展直徑同實(shí)驗(yàn)對(duì)比，誤差小于10%，證明該模型可適用于低黏度、低密度的單液滴撞擊水平冷板面的情況．

單液滴；冷板面；模擬；水平集方法；固化模型

噴霧冷凍干燥是一種新型的工程技術(shù)，主要用于奶粉、藥物蛋白、卵清蛋白、流感疫苗[1–4]等高附加值熱敏性產(chǎn)品的制備，產(chǎn)品形態(tài)也多種多樣，合理地優(yōu)化其過(guò)程工藝有助于制備比表面積大、空隙率高、結(jié)構(gòu)均勻穩(wěn)定的粉體[5]．霧化液滴在設(shè)備壁面的冷凍沉積是在噴霧冷凍干燥技術(shù)應(yīng)用中的常見(jiàn)問(wèn)題．此問(wèn)題可以通過(guò)實(shí)驗(yàn)方法解決，如Xu等[6–7]對(duì)單液滴撞擊不同材質(zhì)的冷板面進(jìn)行實(shí)驗(yàn)，得出不同溫度和不同撞擊速度下單液滴的鋪展直徑和鋪展因子，提出了減少冷凍沉積發(fā)生的方法；還可以通過(guò)模擬設(shè)備流場(chǎng)的方法解決冷凍沉積問(wèn)題，如Anandharamakrishnan等[8]利用3,D計(jì)算流體動(dòng)力學(xué)模型模擬噴霧冷凍干燥設(shè)備的內(nèi)流場(chǎng)，以減小霧化液滴與設(shè)備壁面的接觸幾率；也可直接對(duì)單液滴撞擊冷板面進(jìn)行模擬研究，反映液滴在板面的鋪展凍結(jié)過(guò)程，可以彌補(bǔ)實(shí)驗(yàn)法無(wú)法直接測(cè)得的數(shù)據(jù)，如耿縣如等[9]以單個(gè)霧化液滴為研究對(duì)象，模擬了液滴大小、氣體流速和環(huán)境溫度等參數(shù)對(duì)液滴凍結(jié)過(guò)程的影響，發(fā)現(xiàn)液滴越大凍結(jié)時(shí)所需的形核時(shí)間和完全固化時(shí)間越長(zhǎng)，而且凍結(jié)過(guò)程隨著氣體流速的增大和環(huán)境溫度的降低而縮短，為分析液滴撞擊板面后的凍結(jié)過(guò)程提供了理論依據(jù)．

液滴撞擊冷板面的模擬可以認(rèn)為是存在相變的自由表面流動(dòng)問(wèn)題，涉及移動(dòng)界面．同時(shí)，考慮凝固相變時(shí)，對(duì)流體有效的控制方程對(duì)固體不再適用，自由表面的流動(dòng)也會(huì)變得更加復(fù)雜，需要追蹤相界面的變化．

目前對(duì)于處理液體凝固的計(jì)算方法有2種，其中1種是對(duì)自由表面及流體中的焓值變化量進(jìn)行追蹤．Pasandideh-Fard等[10]采用改進(jìn)的VOF法模擬了錫液滴撞擊溫度為25～240℃的不銹鋼板面的凝固過(guò)程，并且模型中引入了液滴接觸角及液滴底部界面的傳熱系數(shù)，該系數(shù)由實(shí)驗(yàn)測(cè)得并作為模擬的邊界條件．但是VOF方法由于數(shù)值不連續(xù)，所以在計(jì)算界面的位置和平均曲率(確定表面張力)時(shí)精度較低.第2種方法是將流體描述為多孔介質(zhì)[11–13]．

Voller等[14]使用固定網(wǎng)格來(lái)模擬對(duì)流擴(kuò)散條件下的液體凝固情況，將焓值的變化轉(zhuǎn)化為關(guān)于溫度的函數(shù)，且定義全液相時(shí)孔隙率為1、全固相時(shí)孔隙率為0，將流體速度與孔隙率相乘，當(dāng)流體速度為0時(shí)，則可以認(rèn)為流體固化完全，通過(guò)在Navier-Stokes方程中修正體積力源項(xiàng)來(lái)模擬相變過(guò)程．而Marin[15]利用此方法模擬了直徑2mm的鋁液滴撞擊板面的固化過(guò)程，結(jié)果表明此方法對(duì)模擬高黏度、高密度且具有高熔化潛熱的材料撞擊板面的問(wèn)題具有很好的適用性．

本文結(jié)合Sussman等[16]描述的水平集方法與Marin[15]使用的固化模型，使用基于溫度的熱傳導(dǎo)方程代替原有的焓值方程，并且在熱容源項(xiàng)中加入熔化潛熱，以解決低黏度、低密度液滴(如蒸餾水等)撞擊板面的相變問(wèn)題．

1 控制方程與計(jì)算方法

1.1 流體流動(dòng)控制方程

液滴與空氣的速度場(chǎng)和壓力場(chǎng)控制方程用不可壓縮Navier-Stokes方程表示：

式中：?為哈密爾頓算符；t為時(shí)間；μ為流體黏度；ρ為流體密度；u為速度矢量；ui為速度矢量在i方向上的數(shù)值，i為r或z，r表示液滴鋪展直徑方向，z表示液滴鋪展高度方向；p為流體壓強(qiáng)；iF為體積力在i方向上的數(shù)值．氣液兩相界面用固定網(wǎng)格的歐拉方程處理，由于氣體和液體之間的密度差，所以氣液接觸面存在表面張力的影響．體積力包括重力和水平集法的表面張力，由于本模型中存在液滴固化的情況，因此體積力還包含依賴于凝固過(guò)程的固體比例源項(xiàng)，該源項(xiàng)將在1.3節(jié)詳細(xì)說(shuō)明．

1.2 水平集方法

水平集方法可以用于求解包含表面張力相的不可壓縮兩相流動(dòng)．在水平集方法中，使用水平集函數(shù)?來(lái)描述兩相交界面．在連續(xù)相中水平集函數(shù)始終為正，在分散相中始終為負(fù)，兩相交界的自由表面由水平集函數(shù)為零的點(diǎn)構(gòu)成，即：連續(xù)相?(x, y, t)＞0；交界面?(x, y, t)＝0；分散相?(x, y, t)＜0．

從自由面定義可知，交界面上的單位法向量n由分散相指向連續(xù)相，交界面的曲率κ可以用水平集函數(shù)表示為

交界面的運(yùn)動(dòng)可通過(guò)水平集函數(shù)的對(duì)流表示為

描述物性急劇變化的Heaviside函數(shù)可以用水平集函數(shù)表示為

Heaviside函數(shù)作用在液相與氣相界面間平滑過(guò)渡，液相賦值為1、氣相賦值為0．流體不同屬性的改變可以用Heaviside的函數(shù)分別進(jìn)行計(jì)算，公式見(jiàn)式(7)—式(9)．

式中：k為導(dǎo)熱系數(shù)；下標(biāo)g表示氣體，l表示液體．

表面張力在Navier-Stokes方程中只作用在氣液交界面，由式(10)定義．

式中：σ為表面張力系數(shù)；κ為曲率；n為液–氣界面的單位法向量．

采用平滑Dirac delta函數(shù)擴(kuò)展出給定界面寬度的近似方程：

式中，n2為此函數(shù)的參數(shù)．

1.3 固化控制方程

純凈流體是否發(fā)生凝固由熔點(diǎn)決定．在流體冷卻過(guò)程中，一旦溫度低于熔點(diǎn)，在溫度繼續(xù)降低之前流體開(kāi)始釋放相變潛熱．在多組分體系中，相變發(fā)生在從固體開(kāi)始出現(xiàn)到流體最后凝固的溫度范圍內(nèi)．在此過(guò)程中熔化潛熱不斷釋放．

傳熱方程用系統(tǒng)溫度形式表示為

首先，將相變時(shí)的固體組分FS表示為與溫度相關(guān)的函數(shù)：

式中：T為溫度；Tm為蒸餾水發(fā)生凝固溫度；ε為凍結(jié)溫度范圍中間值．

為了將固化過(guò)程并入體積力中，可以將液相看作是一種多孔介質(zhì)，孔隙率λ依賴于溫度，全液相對(duì)應(yīng)于孔隙率為1，全固相對(duì)應(yīng)于孔隙率為0，孔隙率定義為

對(duì)于固化過(guò)程，固相源項(xiàng)S定義為

式中，A是孔隙率的函數(shù)，根據(jù)式(16)定義．

式中，C、q均為函數(shù)A的參數(shù)．

當(dāng)溫度位于液相線以上時(shí)，系統(tǒng)處于全液態(tài)，固體比例源項(xiàng)取值為0，對(duì)Navier-Stokes方程沒(méi)有任何改動(dòng)．當(dāng)溫度處于液相線和固相線之間，固液共存為多孔區(qū)域，A值增大并影響瞬態(tài)、對(duì)流和擴(kuò)散項(xiàng)，流動(dòng)方程近似相當(dāng)于多孔介質(zhì)中的Darcy定律．當(dāng)溫度降低到固相線以下時(shí)，孔隙率接近于0時(shí)(固相)，該源項(xiàng)使得速度值趨于0(固相)．式中的C和q是能達(dá)到預(yù)期的效果的任意常數(shù)，通常C要足夠大且q要足夠小到避免被0除．

水平集函數(shù)定義了兩相初始態(tài)(該模擬中是氣相與液相)的交界面，當(dāng)固化發(fā)生時(shí)，固相會(huì)出現(xiàn)在液相中，所以需要修改固相源項(xiàng)S，使其嚴(yán)格作用在液相中．FS項(xiàng)在方程中都要乘以Heaviside函數(shù)H(?)(在液態(tài)時(shí)H(?)＝1)，這樣就確保了固化模型是液相獨(dú)有，無(wú)論氣體溫度如何，模型都不會(huì)影響氣相．

式(1)中的體積力最終可以表示為

總體系熱量由兩部分組成：顯熱h和潛熱ΔH．顯熱由式(18)計(jì)算．

式中，cp為比定壓熱容．

潛熱的釋放與溫度的下降呈線性變化，發(fā)生相變的焓值作為溫度的函數(shù)可以寫(xiě)成式(19)形式．

釋放潛熱可以被看作是只在過(guò)渡溫度(從液相到固相)范圍下釋放的顯熱，因此潛熱熱容可以定義為

改進(jìn)后的熱容公式為

2 數(shù)值模擬

2.1 材料屬性

本文模擬初始溫度為293K、直徑為3.2mm的蒸餾水液滴撞擊溫度為268K和253K的不銹鋼板面，空氣為氣相．在模型中，空氣的溫度設(shè)置為與不銹鋼板面的溫度一致，且均低于蒸餾水的凍結(jié)溫度.材料的基本性質(zhì)見(jiàn)表1．

表1 空氣、蒸餾水、不銹鋼板的基本性質(zhì)Tab. 1The basic properties of air，distilled water and stainless steel plate

2.2 參數(shù)與變量

由于模型中液滴存在相變情況，因此需要加入額外的源項(xiàng)來(lái)控制相變．在模型中加入的控制變量如表2所示．

模型建立的過(guò)程中需要進(jìn)行全局定義的與模擬過(guò)程相關(guān)的參數(shù)如表3所示．

表2 相變控制變量Tab. 2 Variables in the phase change control equation

表3 全局定義參數(shù)Tab. 3 Globally defined parameters

2.3 模擬過(guò)程

利用He等[17]改進(jìn)的經(jīng)驗(yàn)公式，計(jì)算20℃的蒸餾水液滴在板面上發(fā)生附著或破碎與反彈之間的臨界韋伯?dāng)?shù)，可以得到臨界撞擊速度為2.08m/s，若撞擊速度大于臨界速度，液滴發(fā)生附著或破碎，反之液滴發(fā)生反彈．噴霧冷凍干燥箱中的溫度可以控制在269～233K．故設(shè)置液滴分別從H＝100、250、500mm高度下落，撞擊溫度為253K和268K的不銹鋼板面，對(duì)應(yīng)撞擊速度分別為1.4、2.2、3.1m/s．

模擬采用美國(guó)COMSOL公司的COMSOL Multiphysics 4.2軟件中“流體流動(dòng)”模塊下的“兩相流水平集物理場(chǎng)”與“傳熱”模塊下的“流體傳熱物理場(chǎng)”耦合．考察的是撞擊后的液滴鋪展直徑以及發(fā)生相變的過(guò)程，因此計(jì)算區(qū)域選擇為包含1個(gè)液滴的流體區(qū)域．設(shè)置環(huán)境與板面溫度相同，因?yàn)榭諝鈧鳠嵯禂?shù)較小且下落時(shí)間較短，為簡(jiǎn)化計(jì)算假設(shè)液滴下落過(guò)程中液滴溫度不受影響．由于液滴呈球形結(jié)構(gòu)并且板面相對(duì)液滴無(wú)限大，故將整個(gè)液滴與板面簡(jiǎn)化為二維軸對(duì)稱模型．計(jì)算區(qū)域的初始設(shè)置見(jiàn)圖1．

圖1 液滴在板面鋪展的計(jì)算區(qū)域設(shè)置Fig. 1Computational domain of droplet spreading on the plate

由于相變過(guò)程主要發(fā)生在液滴內(nèi)部以及板面上，且初始界面為液滴邊界，因此在劃分網(wǎng)格的過(guò)程中需要在這些關(guān)鍵界面處將網(wǎng)格細(xì)化，如圖2所示．在劃分網(wǎng)格后，進(jìn)行求解器設(shè)置，對(duì)模型進(jìn)行“相初始化”求解和“瞬態(tài)求解器”求解．

圖2 液滴在板面鋪展的網(wǎng)格劃分Fig. 2 Meshing of droplet spreading on the plate

3 模擬結(jié)果與討論

3.1 板面溫度對(duì)液滴撞擊板面的影響

圖3(a)和圖3(b)為從100mm高度下落液滴分別撞擊溫度為268,、253K板面時(shí)的液滴鋪展過(guò)程，圖4為相應(yīng)的液滴鋪展直徑的對(duì)比．從圖4可以看出：對(duì)于同一下落高度，268K對(duì)應(yīng)液滴的最大鋪展直徑要大于253K的最大鋪展直徑，且在碰撞板面后，268K對(duì)應(yīng)的液滴與253K對(duì)應(yīng)的液滴相比回縮更快，說(shuō)明268K不足以使液滴在板面上快速凍結(jié)，從而液滴鋪展的直徑更大，回縮也更明顯；而從253K的鋪展直徑可以看出，液滴達(dá)到最大鋪展直徑后，回縮速度比較慢，說(shuō)明板面的低溫使液滴與板面接觸部分產(chǎn)生了凝固現(xiàn)象，減緩了液滴由于表面張力造成的回縮現(xiàn)象．

圖3 液滴鋪展過(guò)程模擬圖Fig. 3 Simulated droplet spreading process on the cold plate

圖4 不同冷板溫度下液滴鋪展直徑(H＝100mm)Fig. 4 Droplet spreading diameter under varying cold plate temperatures(H＝100mm)

3.2 撞擊高度對(duì)液滴撞擊板面的影響

圖3(b)、圖3(c)和圖3(d)是不同高度下落的液滴撞擊溫度為253K板面的液滴鋪展過(guò)程．從圖中可以看出：液滴在碰撞板面后開(kāi)始鋪展，液滴直徑逐漸增大，在此過(guò)程中凝固阻礙液滴充分展開(kāi)；撞擊高度大時(shí)，液滴的鋪展直徑也相應(yīng)增大．

從圖5可以看出：對(duì)于同一板面溫度(253K)，不同撞擊高度對(duì)應(yīng)的最大鋪展直徑也不相同，對(duì)于從100mm高度下落的液滴，其最大鋪展直徑為9.1mm；對(duì)于250mm高度下落的液滴，其最大鋪展直徑為12.5mm，鋪展到最大直徑的時(shí)間最短；而對(duì)于從500mm高度下落的液滴，其最大鋪展直徑為14.6mm，鋪展到最大直徑的時(shí)間較長(zhǎng)．這是因?yàn)橐旱蔚南侣涓叨仍礁?，撞擊速度增大，液滴具有了更高的撞擊能量．?duì)于253K的板面溫度，液滴在回縮過(guò)程中也在不斷固化，雖然液滴隨下落高度增加具有更大的動(dòng)能，但是固化過(guò)程阻礙其回縮到較小直徑；液滴回縮的趨勢(shì)表明，更大高度的撞擊會(huì)最終造成液滴沉積直徑的增大．

圖5 不同撞擊高度下液滴鋪展直徑(T＝253K)Fig. 5Droplet spreading diameter under varying impact heights(T＝253K)

3.3 最大鋪展直徑的模擬結(jié)果與實(shí)驗(yàn)值的對(duì)比

液滴在板面上的鋪展面積(最大鋪展直徑)是研究冷凍沉積時(shí)的重要參數(shù)，將模擬結(jié)果的最大鋪展直徑值與Xu等[5]的實(shí)驗(yàn)值進(jìn)行對(duì)比，結(jié)果見(jiàn)表4．可以看出，本文采用的模型與實(shí)驗(yàn)值的吻合性較好．

表4 最大鋪展直徑模擬值與實(shí)驗(yàn)值Tab. 4Simulated and experimental value of maximum spreading diameter

4 結(jié) 語(yǔ)

本文利用COMSOL Multiphysics 4.2軟件的多物理場(chǎng)耦合功能，使用兩相流水平集與傳熱物理場(chǎng)相結(jié)合的方法，建立了水平集方法的固化模型，對(duì)蒸餾水單液滴撞擊水平冷板面的鋪展和凍結(jié)固化進(jìn)行了模擬．使用基于溫度的熱傳導(dǎo)方程代替焓值方程，加入熔化潛熱修正了熱容源項(xiàng)，可以解決低黏度、低密度單液滴(蒸餾水)撞擊水平冷板面的固化問(wèn)題．

模擬對(duì)比了蒸餾水液滴在不同撞擊高度以及不同溫度板面產(chǎn)生的相變情況，發(fā)現(xiàn)隨著撞擊高度的增大液滴的最大鋪展直徑增大，最終造成液滴沉積直徑的增大；而板面溫度越低，液滴產(chǎn)生固化現(xiàn)象越明顯，其回縮過(guò)程越緩慢，回縮程度越?。?/p>

最大鋪展直徑的模擬值與實(shí)驗(yàn)值吻合較好，誤差小于10%；撞擊過(guò)程中，撞擊高度對(duì)液滴的鋪展和固化的影響較明顯．

［1］ Rogers S，Wu W D，Saunders J，et al. Characteristics of milk powders produced by spray freeze drying[J]. Drying Technology，2008，26(4)：404–412.

［2］ Leuenberger H，Plitzko M，Puchkov M. Spray freeze drying in a fluidized bed at normal and low pressure[J]. Drying Technology，2006，24(6)：711–719.

［3］ Yeom G S，Song C S. Experimental and numerical investigation of the characteristics of spray-freeze drying for various parameters：Effects of product height，heating plate temperature，and wall temperature[J]. Drying Technology，2010，28(2)：165–179.

［4］ Maa Y F，Ameri M，Shu C，et al. Influenza vaccine powder formulation development：Spray-freeze-drying and stability evaluation[J]. Journal of Pharmaceutical Sciences，2004，93(7)：1912–1923.

［5］ 徐慶，耿縣如，李占勇. 噴霧冷凍干燥對(duì)顆粒產(chǎn)品形態(tài)的影響[J]. 化工進(jìn)展，2013，32(2)：270–275.

［6］ Xu Qing，Li Zhanyong，Wang Jin，et al. Characteristics of single droplet impact on cold plate surfaces[J]. Drying Technology，2012，30(15)：1756–1762.

［7］ 徐慶，李占勇，王瑞芳，等. 單液滴撞擊球形表面的涂覆效果[J]. 天津科技大學(xué)學(xué)報(bào)，2013，28(1)：50–54.

［8］ Anandharamakrishnan C，Gimbun J，Stapley A G F，et al. Application of computational fluid dynamics(CFD) simulations to spray-freezing operations[J]. Drying Technology，2010，28(1)：94–102.

［9］ 耿縣如，徐慶，李占勇，等. 噴霧冷凍法單個(gè)液滴凍結(jié)過(guò)程模擬[J]. 化工進(jìn)展，2012，31(5)：981–986.

［10］ Pasandideh-Fard M，Bhola R，Chandra S，et al. Deposition of tin droplets on a steel plate：Simulation and experiments[J]. International Journal of Heat and Mass Transfer，1998，41：2929–2945.

［11］ Voller V R，Markatos N C，Cross M. Techniques for accounting for the moving interface in convection/diffusion phase change[M]//Lewis R W，Morgan K. Numerical Methods in Thermal Problems. Swansea：Pineridge Press，1985：595–609.

［12］ Voller V R，Markatos N C，Cross M. Solidification in convection and diffusion[M]//Markatos N C，Tatchell D G，Cross M，et al. Numerical Simulations of Fluid Flow and Heat/Mass Transfer Processes. Berlin：Springer，1986：425–432.

［13］ Voller V R，Cross M，Markatos N C. An enthalpy method for convection/diffusion phase changes[J]. International Journal for Numerical Methods in Engineering，1987，24(1)：271–284.

［14］ Voller V R，Prakash C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. Journal of Hat and Mass Transfer，1987，30(8)：1709–1719.

［15］ Marin T L. Solidification of a liquid metal droplet impinging on a cold surface[C]//Proceedings of the COMSOL Users Conference. Trondheim，Norway：COMSOL AS，2006：1492-1–1492-8.

［16］ Sussman M，Smereka P，Osher S. A level set approach for computing solutions to incompressible two-phase flow[J]. Journal of Computational Physics，1994，114：146–159.

［17］ He Zheng，Gao Ye，Gu Xuan，et al. Investigating a droplet-wall collision model[J]. Journal of Harbin Engineering University，2009，30(3)：267–270.

責(zé)任編輯：常濤

Numerical Simulation of the Impact of Single Droplet on Cold Horizontal Plate Surface

WANG Jin，LI Zhanyong，XU Qing，WANG Ruifang
(College of Mechanical Engineering，Tianjin University of Science & Technology，Tianjin 300222，China)

In spray-freeze drying(SFD)，freeze deposition of liquid droplets on the equipment walls has become a new challenge to its applications. Based on the previous experimental results of droplet collision with cold horizontal solid plates，a model of a single droplet collision with cold solid plates was established with the combination of the Level Set method and the solidification model to simulate the impacting process of the distilled water droplet(3.2mm in diameter，at 293K)on the stainless steel plate(268K or 253K) from the height of 100mm，250mm and 500mm，respectively. The simulation results of maximum spreading diameter were compared with the experimental data and the errors were less than 10%，which indicates that the model is suitable for dealing with the impact of droplets of low viscosity and low-density on the cold horizontal plate surface.

single droplet；cold plate surface；simulation；level set method；solidification model

TQ051

A

1672-6510(2014)01-0036-06

10.13364/j.issn.1672-6510.2014.01.008

2013–10–07；

2013–11–15

國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)資助項(xiàng)目(2012AA021303)

王 瑾（1990—），女(滿)，陜西人，碩士研究生；通信作者：李占勇，教授，zyli@tust.edu.cn.