梁坤峰，阮春蕾，王林，袁竹林

(1 河南科技大學制冷與空調(diào)技術研究所，河南 洛陽 471003；2 東南大學能源與環(huán)境學院，江蘇 南京 210096)

引 言

基于流態(tài)化原理實現(xiàn)直接接觸式制取流體冰的方法是一種新型的動態(tài)制冰方式[1-2]，制冰過程中循環(huán)流動載冷液體與液滴的直接接觸換熱，具有較高的熱效率。為了研究流態(tài)化制冰過程中的顆粒運 動、相變與碰撞聚并，以及循環(huán)載冷液體的流場與溫度場分布特征，掌握制冰過程的控制方法，需要探討流化床內(nèi)液滴、冰顆粒和循環(huán)載冷液體所構成的顆粒-流體多相系統(tǒng)的流動與傳熱特性。但由于多相系統(tǒng)相間相互作用的復雜性，存在觀測困難，單純采用實驗手段進行其內(nèi)部規(guī)律性研究，難以完整闡釋。隨著計算機的高速發(fā)展，數(shù)值模擬方法以其較強的再現(xiàn)能力、低耗費、零危險性等優(yōu)勢，在對多相流動與換熱體系的研究中發(fā)展很快。

近年來，一種改進的基于歐拉-拉格朗日方法的CCDM（combined continuum and discrete model）模型開始應用于顆粒-流體多相系統(tǒng)的研究[3-6]，該模型結合了用于計算顆粒的離散相模型（discrete phase model，DPM）和模擬連續(xù)相的計算流體動力學模型，既采用歐拉法對連續(xù)相的宏觀場特征進行求解，又基于歐拉-拉格朗日法跟蹤求解場中的每一個真實顆粒，反映顆粒在流動過程中的微觀細節(jié)，很好地描述顆粒-流體多相系中顆粒運動、傳熱與傳質(zhì)的微觀特性。自從Tsuji 等[7]首先提出CCDM 模型以來，已在顆粒-流體多相系統(tǒng)模擬研究中得到關注；Swasdisevi等[8]對二維噴動床內(nèi)的氣固傳熱特性進行了研究，考慮了氣固相間的對流傳熱，并在噴動區(qū)和下降管中采用不同的傳熱系數(shù)；Li 等[9]考慮了顆粒-流體間的對流傳熱和顆粒-顆粒間的熱傳導，基于CCDM模型研究了水平管中的氣固流動與傳熱問題；Zhou等[10]采用CCDM 模型研究了氣固流化床內(nèi)的煤粉燃燒問題，對氣相湍流采用了大渦模擬方法，相間對流傳熱和顆粒-顆粒的接觸導熱也得以考慮；而采用CCDM 模型在液-固兩相流體系的流動與傳熱方面的研究較為缺乏，Malone 等[11]采用準三維的CCDM模型研究了液-固流化床內(nèi)的流動與換熱問題，比較了兩個不同的相間傳熱系數(shù)模型的應用條件。

本文建立了描述顆粒-流體多相系統(tǒng)的流動、傳熱與相變的CCDM 模型，基于歐拉法求解連續(xù)相的宏觀場特征，由歐拉-拉格朗日法跟蹤多相場中的離散顆粒相的運動、換熱與相變特性；采用兩套不同尺度的計算網(wǎng)格準確描述宏觀場特征與微觀顆粒行為，實現(xiàn)雙向耦合作用；獲得了流態(tài)化床內(nèi)顆粒-流體多相系統(tǒng)的流動與傳熱特性，所獲結論為流態(tài)化制冰過程的運行設計提供了理論指導。

1 CCDM 模型建立

1.1 基本方程

1.1.1 連續(xù)相守恒方程

（1）質(zhì)量守恒方程

（2）動量守恒方程 在流化床內(nèi)，由于離散相液滴與連續(xù)相載冷液體間局部的流體溫度梯度將引起局部的流體密度差，考慮到床內(nèi)多相流動的溫度梯度不太大，采用Boussinesq 假設，在動量守恒方程中的重力項中將流體的密度處理為溫度的線性函數(shù)[11]，控制方程為

（3）能量守恒方程 流場局部平均的流體能量守恒方程為

1.1.2 離散相的運動與熱平衡方程 床內(nèi)離散顆粒相的運動與傳熱由DEM 模型求解。根據(jù)牛頓運動定律，顆粒相的運動方程為

忽略熱輻射傳遞量、顆粒-顆粒碰撞接觸的導熱量，離散顆粒的熱平衡方程為

在較小的時間尺度Δt內(nèi)，假設傳遞給顆粒i的熱量保持不變，則有

1.2 顆粒-流體相間全耦合機制

1.2.1 力的耦合 與氣-固流動體系相比，液-固兩相流的相間作用力要更復雜，需要考慮的力有：穩(wěn)態(tài)流體曳力、虛擬質(zhì)量力、Basset 力、剪切提升力、旋轉提升力、壓力梯度力、浮力?；谖闹醒芯康念w粒-多相系統(tǒng)特征，可忽略Basset 升力和剪切提升力。相間作用力表達式如下

表1 相間作用力系數(shù)Table 1 Equation of force between liquid and solid phases

式中系數(shù)與相關指數(shù)說明見表1。

1.2.2 熱量的耦合 對流傳熱通常是連續(xù)流體相中主要的傳熱模式，傳熱公式常常是建立在實驗數(shù)據(jù)上的經(jīng)驗或半經(jīng)驗公式。對于顆粒-流體多相流動體系，Haid[16]根據(jù)實驗數(shù)據(jù)擬合出了三系數(shù)模型，預測的相對誤差為34.4%；Bremford 等[17]實驗研究了床內(nèi)黏度相對較大的連續(xù)相流體中的相間傳熱系數(shù)，指出固相顆粒的存在對流體黏度和相間傳熱的影響是不可忽略的；而Agarwal[18]通過引入床層空隙率函數(shù)，將描述單個球形顆粒在均勻流場中強制對流的Ranz-Marshall 關聯(lián)式應用到顆粒-流體多相系統(tǒng)中，實現(xiàn)了流化床體系中顆粒對流傳熱的Nusselt 數(shù)和在相同表觀速度流體速度下單個顆粒的Nusselt 數(shù)的關聯(lián)，并以此計算床內(nèi)液滴的對流傳熱系數(shù)

離散顆粒和周圍流體間的對流換熱由下式確定

1.2.3 顆粒的碰撞 流化床內(nèi)液滴與載冷液體之間的流動換熱，促使其表層快速凍結，凍結一定程度后的液滴發(fā)生碰撞時，難以發(fā)生顆粒的聚并，類似彈性剛體間的碰撞，服從胡克定律。由此，文中通過判定液滴凍結程度的深或淺，采用不同的碰撞模型計算，液滴凍結程度淺時，碰撞采用經(jīng)典的O’Rourke 定律，而凍結程度深時，液滴間的碰撞則采用軟球模型。

當連續(xù)相為黏度較大的液體介質(zhì)時，由于液體介質(zhì)在顆粒表面所產(chǎn)生的流動潤滑作用，顆粒之間的碰撞與在氣固流動體系存在著較大的差異。Joseph 等[19]以及Malone 等[11]通過對顆粒的恢復系數(shù)e進行修正，來考慮黏性液體介質(zhì)對顆粒碰撞產(chǎn)生的影響

式中，Stc為顆粒碰撞產(chǎn)生反彈的關鍵Stokes損傷值，計算中取Stc=10[20]；Sti是與黏性流場相關的量

1.3 求解方法

1.3.1 計算條件 數(shù)值模擬的對象為實驗研究裝置中的測試床段，床體的內(nèi)徑為40 mm，考慮到計算時間，根據(jù)所考察內(nèi)容，將床體高分為2 m 與8 m（實驗裝置中，豎直床體段高2 m，而從液滴霧化形成的床體底部到冰晶與載冷液體混合物出口處的全距離為8m）兩種情況。

整個求解區(qū)域劃分非結構六面體網(wǎng)格，網(wǎng)格總數(shù)分別為8000（2 m）和32000（8 m）。床體的入口為速度入口邊界條件，出口為壓力出口邊界條件，而壁面采用無滑移、絕熱邊界條件處理。顆粒的噴入位置為床體底部0～0.03 m 區(qū)間，采用隨機方法確定噴入顆粒的初始位置，應用Rosin-Rammler 分布函數(shù)確定噴入顆粒的初始粒徑分布。計算中所設計到的相關參數(shù)及條件見表2。

1.3.2 分割式搜索法 目前，對于離散顆粒受力以及顆粒對流場反作用的求解，較為普遍的是采用一種簡單近似、結構網(wǎng)格統(tǒng)計平均的方法，但對于具有復雜外形幾何體的非結構網(wǎng)格而言，難免會產(chǎn)生誤差，甚至極大地偏離真實情形。文中采用分割式搜索法確定顆粒所在的當?shù)胤墙Y構網(wǎng)格，顆粒-流體間相互作用均在該非結構網(wǎng)格中進行計算，實現(xiàn)一對一的雙向耦合作用與反饋。以二維網(wǎng)格為例說明如下：將不規(guī)則外形幾何體，用較大的規(guī)則幾何體覆蓋，如圖1(a)所示，不規(guī)則橢圓區(qū)域為求解區(qū)域，規(guī)則四邊形為覆蓋區(qū)域。將橢圓區(qū)域內(nèi)的非結構網(wǎng)格采用一維編號技術進行編號[21]，并計算其所在結構網(wǎng)格的編號。在求解區(qū)域內(nèi)，每個結構網(wǎng)格中包含至少1 個不規(guī)則非結構網(wǎng)格，然后對其中的非結構網(wǎng)格依據(jù)掃描順序重新進行編號，如圖1(b)所示，并建立數(shù)組進行存儲。在每個求解時間步，根據(jù)顆粒物理位置，確定其所在的結構網(wǎng)格。搜索僅在該結構網(wǎng)格及其鄰近結構網(wǎng)格中進行，若滿足lc=min(l1,l2,…,ln)，即為顆粒所在的當?shù)胤墙Y構網(wǎng)格，其中l(wèi)為非結構網(wǎng)格質(zhì)心到顆粒的距離。采用該方法，能準確地計算顆粒-流體間的耦合作用與反饋，且搜索效率顯著提高。

表2 模型的計算條件和參數(shù)Table 2 Calculation conditions and calculated parameters

圖1 非結構網(wǎng)格顆粒搜索示意Fig.1 Non-structural grid used for searching particles

1.3.3 計算流程 采用先分離后耦合的算法，即顆粒相與連續(xù)相場分開計算，而耦合計算固相顆粒場 與連續(xù)相流場之間的相互作用。計算步驟為：首先，獲得一個穩(wěn)定的初始連續(xù)相流場；接著將離散顆粒噴入連續(xù)相流場中，根據(jù)顆粒所在的位置計算整個流場區(qū)域中每個控制體積的空隙率（即連續(xù)相所占的體積分數(shù)）；然后計算離散顆粒所受的作用力和顆粒-流體之間的對流傳熱量，顆粒所受力包括流體作用力、顆粒間碰撞接觸力、重力等；隨后依據(jù)牛頓第二定律，在多個固相顆粒計算時間步內(nèi)，更新顆粒的速度、位置、溫度等屬性，并在固相顆粒場計算結束時，再次計算流場中每個控制體積的空隙率，獲得固相顆粒對流體的反作用力和總的換熱量；最后以源項的形式將其加入連續(xù)相場的求解控制方程中，迭代求解連續(xù)相流場，進入下一循環(huán)。

2 實驗裝置

流態(tài)化制冰實驗系統(tǒng)原理圖如圖2所示，整個實驗系統(tǒng)由流化床、制冷機和測試裝置3 部分組成。流化床的主體由有機玻璃管制成，流化床底部中心軸線處設置1 個噴射制冰液體的噴頭，噴頭以下 300 mm 處為另一種非相溶的循環(huán)流動的載冷液體入口。制冷機（5 hp 壓縮機，1 hp=745.699 W）通過蒸發(fā)器給吸收制冰液體相變熱的載冷液體提供 足夠的冷量。測試裝置包括溫度[溫度傳感器采 用Pt100 型鎧裝熱電阻，精度等級為A，允許偏 差±(0.15+0.002|t|)]、壓力與流量（壓力測量選用絕壓變送器，靜態(tài)精度為0.5%FS；載冷液體的流量采用LC-40 型齒輪流速計，精度等級為0.5，制冰液體流量由LZJ-10 型玻璃轉子流量計測量）的測試與采集。

制冰過程在流化床內(nèi)部進行，床的底部形成兩個并列兩路液體入口：循環(huán)載冷液體流動的油路與制冰液體噴射霧化的水路。循環(huán)載冷液體流路中，經(jīng)制冷裝置蒸發(fā)器降溫后的載冷液體在離心泵的作用下，自蒸發(fā)器流經(jīng)過濾器，最后流入流化床內(nèi)。當載冷液體的溫度、流速穩(wěn)定后，打開增壓泵將水強制噴入載冷液體中，水在非相溶的載冷液體中液-液霧化形成細小液滴，并與其一起向流化床上部同向流動。此時，可通過測試裝置進行實驗測量，如通過設置在床高0、0.5、1.0、1.5、2.0 m 處的溫度傳感器進行溫度測量。

3 結果與分析

3.1 顆粒碰撞模型的判別準則

當液滴在霧化形成的初始階段，由于液滴凍結程度淺，液滴大部分都處于液態(tài)，此時，液滴-液滴的碰撞極有可能引起液滴聚并現(xiàn)象的發(fā)生，可采用O’Rourke 經(jīng)典的液滴碰撞模型進行求解；而當液滴凍結到一定程度后，液滴外殼堅硬，即使顆粒之間發(fā)生碰撞，但聚并的概率極小，此時可遵循剛性球體間碰撞機理進行計算。如何確定液滴凍結的程度，并以此作為兩種碰撞模型使用的判別準則？文中提出冰晶凍結程度關鍵值ξ的定義：ξ=rc/R，式中，rc為等效半徑，即在某一時刻顆粒表面的對流換熱量等效于其由外向里凍結所需冷量時的顆粒半徑；R為顆粒的半徑。

圖3為4 種凍結程度關鍵值下顆粒粒徑的軸向分布情況。為方便比較不同關鍵值對顆粒粒徑的影響，不同工況下冰晶平均粒徑起始值按該床高處的實驗值，統(tǒng)一起始點。由圖可知，在相同床高處，凍結程度關鍵值越大，顆粒的粒徑則越小，且隨著凍結程度關鍵值的增大，液滴粒徑沿床高而增大的幅度明顯減??；表明凍結程度關鍵值較大時，液滴表面的換熱能力強，液滴相變進行得更快，此時液滴間碰撞發(fā)生聚并的概率就小，由此顆粒的粒徑不會有較大幅度的增加。為了確定計算中凍結程度關鍵值的大小，計算了幾個凍結程度關鍵值下顆粒平均粒徑沿床高的變化曲線，通過與相同工況下的實驗結果對比，發(fā)現(xiàn)采用凍結程度關鍵值ξ=10 時的床內(nèi)顆粒平均粒徑沿床高的變化趨勢與試驗結果比 較吻合，因此，在之后的計算分析中統(tǒng)一采用該顆粒凍結程度關鍵值來計算流化床內(nèi)離散顆粒相與連續(xù)油介質(zhì)之間的流動換熱。

圖3 凍結程度關鍵值對顆粒粒徑沿軸向分布的影響Fig.3 Effect of particles freezing degree on axial size distribution (Tcont=-7 ,℃Gcont=0.4 L·s-1,Gdisp=30 ml·min-1)

圖4 載冷液體軸向溫度分布Fig.4 Axial temperature distribution of coolant

3.2 循環(huán)床內(nèi)多相流動與傳熱特性

3.2.1 載冷液體溫度分布 圖4為流化床從起始到穩(wěn)定的不同時刻載冷液體沿軸向的溫度分布曲線。在霧化液滴進入流化床后的很短時間內(nèi)（t=0.5 s），離散液滴還沒來得及分布于整個流化床內(nèi)，與連續(xù)相間只在床的底部換熱，床內(nèi)載冷液體的溫度呈先增大后減小，繼而保持初溫不變；運行5 s 時間后，受離散液滴在床內(nèi)擴散運動的影響，載冷液體的溫度分布曲線峰值朝床高方向移動；運行50 s 時間后，床內(nèi)液-液-固三相流動與換熱體系基本達到穩(wěn)定，載冷液體沿床高的溫度增長趨勢有所平緩。

此外，圖4還給出了相同工況條件下的實驗測試結果，對比發(fā)現(xiàn)模擬與實驗結果的總體變化趨勢一致，最大誤差在13%左右，基本驗證了CCDM 模型對于流態(tài)化制冰過程的模擬能力。分析局部出現(xiàn)誤差的原因有兩方面。一是實驗數(shù)據(jù)的處理方面，實驗的原始數(shù)據(jù)去掉了流化床跑冷的影響，而分析時采用的跑冷系數(shù)本身是根據(jù)實驗數(shù)據(jù)擬合得到的[22]，在不同的階段和工況下易造成誤差的放大和減小；二是傳熱模型的應用方面，文中采用了Agarwal[18]提出的傳熱關聯(lián)式，雖然該模型引入了床層空隙率函數(shù)將流化床體系中顆粒對流傳熱的Nu和在相同表觀速度流體速度下Ranz-Marshall 關聯(lián)式中單個顆粒的Nu關聯(lián)起來，但應用于低溫低黏度載冷液體與顆粒的表面?zhèn)鳠岜尘?，會對計算結果會產(chǎn)生一定影響。

圖5 載冷液體徑向溫度分布Fig.5 Radial distribution of temperature of coolant

圖6 離散顆粒相軸向溫度分布Fig.6 Axial distribution of temperature of discrete particle

圖5為連續(xù)相載冷液體在不同床高處沿徑向的溫度分布曲線。由圖可知，在床的中心區(qū)域，連續(xù) 相沿床高方向的溫度分布基本呈逐漸升高趨勢；而在近壁區(qū)域，載冷液體的溫度分布則比較紊亂，尤其是在床高1 m 處。原因在于床內(nèi)離散顆粒相的紊亂運動，尤其是在床的上半段，由于顆粒向上運動速度的減慢，顆粒數(shù)密度增大，增加了顆粒間碰撞的概率，而此時大部分顆粒處于未凍結完全或凍結程度很淺的狀態(tài)，顆粒之間的碰撞極易引起顆粒間聚并現(xiàn)象的發(fā)生，顆粒間的頻繁碰撞與聚并將使其在床內(nèi)的運動變得紊亂，顆粒的運動速度以及相間的能量耦合作用，都呈較大幅度的波動，這些因素將導致連續(xù)相溫度場分布的不均勻。

3.2.2 顆粒相軸向溫度分布 圖6為離散顆粒相軸向溫度分布曲線，反映了不同床層高度內(nèi)顆粒平均溫度的變化特征。由圖可知，床內(nèi)離散顆粒相溫度分布的波動較為紊亂，但總的變化趨勢仍與水的結晶過程曲線較為一致。顆粒相進入流化床后的很短時間內(nèi)其平均溫度便迅速下降，由初溫5℃降至 0.5℃，這是由于顆粒入床時，顆粒速度與溫度均處于最大值，相間較大的溫度差以及較大的流速比，加劇了相間的傳熱，使顆粒溫度大幅度降低。到達床高2 m 時，床內(nèi)部分顆粒已達到相變溫度。在床高2～7 m 區(qū)間，床內(nèi)大部分顆粒均處于相變傳熱階段，在床高7 m 以后，床內(nèi)大部分顆粒基本上已完全凍結，此時，床層離散顆粒相的溫度降至冰點溫度以下，溫度下降的幅度也比較平緩。

3.2.3 顆粒相軸向速度分布 圖7為不同粒徑時顆粒軸向速度分布曲線，選取了6 種粒徑檔的顆粒作為示蹤顆粒，考察了其在床內(nèi)流動速度的變化情況。此時制冰液體的噴射流率較小，為3.3 ml·min-1，由于噴射流率較小，床內(nèi)顆粒數(shù)密度低，顆粒間發(fā)生碰撞聚并的概率較低。由圖可知，在粒徑檔小于3.12 mm 時，顆粒在床內(nèi)的運動速度呈先急劇減小，后逐漸增大，最終達到其在床內(nèi)流動的終端速度；隨著顆粒粒徑的增大，這種變化趨勢逐漸減弱，當顆粒粒徑檔大于3.62 mm 以后，顆粒的軸向速度則呈持續(xù)增大的變化趨勢；但同時需注意不同粒徑檔顆粒的終端速度差異極小。出現(xiàn)這種情況的原因主要在于顆粒的流化過程中，噴射初始動量、自身重力、載冷液體的攜帶作用之間的相互作用所致；顆粒在床內(nèi)霧化形成后，其運動速度在很短的時間內(nèi)急劇減小，隨后由于其速度的降低，其受到周圍載流液體的流動攜帶作用，運動速度又逐漸增大，最后達到顆粒的終端速度。

圖8為離散顆粒相軸向平均速度沿床高的變化曲線，此時噴射流率較大，為35 ml·min-1。由圖可知，離散相噴射流率較大時，由于床層的固相顆粒數(shù)密度分布增大，床內(nèi)顆粒間碰撞、聚并與結團等現(xiàn)象的發(fā)生概率增加，致使顆粒相在床內(nèi)的運動變得紊亂，其軸向平均速度沿床高方向的分布相應 地也變得紊亂；尤其是在床高2 m 以下，由于顆粒運動速度的減慢，離散顆粒相軸向平均速度變化幅度較大。但從總體變化趨勢看，離散顆粒相的軸向平均速度呈現(xiàn)先增大，隨后在同一水平線上下波動的趨勢。

圖7 不同粒徑時顆粒軸向速度Fig.7 Axial velocity of particles of different radius

圖8 顆粒相軸向平均速度Fig.8 Axial average velocity distribution of particles phase

圖9 顆粒相數(shù)密度軸向分布Fig.9 Axial number density distribution of particle phase

3.2.4 顆粒相濃度分布 圖9為顆粒相數(shù)密度的軸向分布，圖中曲線代表離散數(shù)據(jù)點的變化趨勢。由圖可知，床內(nèi)固相顆粒的數(shù)密度沿床高而降低，在床高0.5～2.5 m 區(qū)間內(nèi)，顆粒數(shù)密度變化幅度較大；達到床高6 m 以后，顆粒相數(shù)密度逐漸趨于平緩，基本保持不變。原因在于該高度區(qū)間內(nèi)液滴未凍結完全，且運動速度減小，大量顆粒匯集在該區(qū)域，增加了顆粒間碰撞并發(fā)生聚并的概率，從而使該區(qū)域離散顆粒相數(shù)密度波動幅度較大。在床高2.5 m 以上，根據(jù)圖8顆粒相軸向平均速度的變化趨勢看，已基本趨于穩(wěn)定，且已有部分顆粒的凍結程度較深，顆粒相的運動相對比較穩(wěn)定，數(shù)密度也趨于穩(wěn)定。

圖10為顆粒相體積分數(shù)軸向分布，圖中曲線代表離散數(shù)據(jù)點的變化趨勢。與圖9相對應，在床高0.5～2.5 m 區(qū)間內(nèi)，由于顆粒的匯集，床內(nèi)顆粒相體積分數(shù)達到最大，這一點與圖8的顆粒相平均速度最小的表現(xiàn)一致；而在床高2.0 m 以上，由于離散顆粒流動的穩(wěn)定，床層固相顆粒的體積份額也相應地達到穩(wěn)定。

圖10 顆粒相體積分數(shù)軸向分布Fig.10 Axial volume fraction distribution of particle phase

圖11 不同粒徑檔顆粒的停留時間Fig.11 Residence time for particles of different radius

3.2.5 顆粒停留時間 顆粒在流化床內(nèi)的停留時間是循環(huán)流化床制冰運行的重要參數(shù)。不同粒徑檔的顆粒在床內(nèi)的運動特性存在差異，決定了其停留時間的差異；而在同一工況的相同操作參數(shù)下，必然會影響顆粒的凍結程度，并進而影響顆粒的碰撞聚并。因對不同粒徑檔的顆粒在床內(nèi)停留時間的分析，有利于對操作參數(shù)及系統(tǒng)的優(yōu)化設計。圖11為不同粒徑檔的顆粒在流化床內(nèi)停留時間的統(tǒng)計平均值，計算時間為100 s，其中，圖中曲線代表離散數(shù)據(jù)點的變化趨勢。由圖可知，顆粒在床內(nèi)的停留時間隨 其粒徑的變化雖較為紊亂，但從總的變化趨勢上看，仍呈現(xiàn)出隨粒徑的增大，顆粒停留時間先增大后減小的變化趨勢，并以粒徑檔3 mm 為分界點。粒徑檔小于3 mm 的顆粒，由于其體積較小，凍結過程所需冷量也少，在流化床下半段的流動過程中與載冷液體換熱即可凍結到較深的程度，顆粒在隨后的流動過程中，即使碰撞，發(fā)生聚并的概率很小，從而在一定程度上減小了小粒徑顆粒流動紊亂性；在此粒徑檔范圍內(nèi)，顆粒的停留時間主要取決于顆粒的大小，大顆粒較小顆粒而言，不易被攜帶，故其停留時間也就較長。然而當顆粒的粒徑檔大于3 mm時，其凍結時間增加很多，使得液滴流過其初始動量較強的流化床下半段，其凍結程度仍較淺；越過該高度后，由于液滴流速的降低，床內(nèi)局部液滴濃度增大，增大了顆粒間的碰撞概率，從而對于凍結程度相對較淺的大顆粒而言，發(fā)生聚并、粘連、分散等現(xiàn)象的概率較大，使其流動變得較為紊亂與復雜。

3.2.6 相間平均傳熱系數(shù) 圖12為相間平均傳熱系數(shù)的軸向分布，通過統(tǒng)計床高0.1 m 高度內(nèi)所有顆粒的數(shù)量和表面?zhèn)鳠嵯禂?shù)，并對這段高度內(nèi)的所有顆粒的表面?zhèn)鳠嵯禂?shù)進行加和平均獲得，體現(xiàn)了流化床內(nèi)相間傳熱系數(shù)隨時間由非平衡而逐漸達到平衡的過程。在流化床運行4.2 s 時，床內(nèi)離散顆粒相流動未穩(wěn)定，床內(nèi)相間平均傳熱系數(shù)沿床高呈先增大后減小的變化趨勢，在床高1 m 處最??；隨后平均傳熱系數(shù)在床高1 m 以上的床層間，沿床高而增大的幅度逐漸減小，最終在t=7.2 s 時達到穩(wěn)定，床內(nèi)平均傳熱系數(shù)沿床高基本不變。分析認為，計算的初始時刻，由于床內(nèi)固相顆粒的數(shù)密度較小，初始噴入顆粒的運動受到顆粒間碰撞聚并現(xiàn)象的影響極小，不存在速度突增或突減的情形；在床的上 半段，由于床內(nèi)連續(xù)相載冷流體與離散顆粒相間的換熱量也很少，其溫度處于較低的水平，相間溫差隨著床高而增大，由此表現(xiàn)為初始時刻平均傳熱系數(shù)在床的上半段沿床高增大。然而隨著時間的變化，床內(nèi)不同高度處的離散相顆粒濃度分布增大且達到穩(wěn)定，顆粒相與連續(xù)相流動穩(wěn)定，使床內(nèi)連續(xù)相載冷液體沿床高的溫度分布達到穩(wěn)定，而相間的平均傳熱系數(shù)也相應地達到穩(wěn)定。

圖12 床層相間平均傳熱系數(shù)隨時間的變化Fig.12 Change of average heat transfer coefficient with time

4 結 論

（1）基于所建立的CCDM 模型求解了流化床制冰過程所涉及的顆粒-多相系統(tǒng)的流動和傳熱特性，該模型充分考慮了顆粒碰撞對流化床內(nèi)液固相間的動量與能量交換的影響，并基于一種高效的非結構網(wǎng)格搜索法，快速地搜索每一計算時刻離散顆粒在流場中所在的非結構網(wǎng)格，實現(xiàn)了連續(xù)相與離散相間一對一的相間耦合作用。

（2）由于顆粒間碰撞、聚并、結團等現(xiàn)象的存在，水的噴射流率與顆粒的數(shù)密度和顆粒發(fā)生碰撞形成大顆粒的概率的變化趨勢一致，而顆粒的粒徑對于顆粒的運動速度和停留時間具有重要影響，從而導致相間傳熱系數(shù)沿床高發(fā)生變化；尤其是在流化床的下半段，顆粒的聚并引起了顆粒運動速度呈2～3 倍的增加，而相間傳熱系數(shù)增大了近百倍，不僅如此，此階段的床內(nèi)固相顆粒的數(shù)密度和體積分數(shù)都比較大。由此可見，在流化床制冰研究中難以忽略顆粒間的碰撞與聚并現(xiàn)象對于流化床的多相流動與換熱特性的影響。

（3）在流化床的底部入口處，處于液態(tài)的顆粒被噴入床中時，顆粒與冷媒之間存在較大的溫度差與流速比，強化了相間的傳熱，顆粒相平均溫度在進入流化床后的很短時間內(nèi)就迅速下降；在隨后的相變換熱階段，顆粒相溫度降低的幅度比較平緩。當流化床處于穩(wěn)定制冰工作時，較大的相間平均傳熱系數(shù)出現(xiàn)在床高0.5 m 以下，超過這個高度，沿床高幾乎不再變化。通過對比載冷液體軸向溫度分布的模擬結果與實驗結果，發(fā)現(xiàn)僅有13%左右的相對誤差，對所建立的CCDM 模型在模擬流化床制冰過程的能力方面進行了驗證，同時表明文中提出的凍結程度關鍵值，對于判定顆粒間碰撞模型和顆粒相表面?zhèn)鳠崃康牡刃в嬎愫唵斡行А?/p>

符 號 說 明

Ac——循環(huán)床層截面積，m2

CA,i——虛擬質(zhì)量力系數(shù)

CD,i——流體曳力系數(shù)

CS,i——旋轉提升力系數(shù)

cf——連續(xù)相流體的比定壓熱容，J·kg-1·K-1

cp,l——水的比定壓熱容，J·kg-1·K-1

cp,s——冰的比定壓熱容，J·kg-1·K-1

dd——霧化液滴粒徑，mm

F——顆粒的受力，N

Ffp——流體-顆粒間動量交換項，N

fA,i——虛擬質(zhì)量力，N

fB,i——相間密度差引起的浮力，N

fD,i——流體曳力，N

fP,i——流體壓力梯度力，N

fS,i——旋轉提升力，N

Gdisp——離散相流率，ml·min-1

Hb——循環(huán)床高度，m

h——相間對流傳熱系數(shù)，J·m-2·K-1

Ip——顆粒轉動慣量，J·kg-1·m-2

Mp——顆粒所受合力矩，N·m

mp——顆粒質(zhì)量，kg

Nu——Nusselt 數(shù)

Pr——Prandtl 數(shù)

p——連續(xù)相壓力，Pa

Qfp——流體與顆粒能量交換項，J

q——相間傳熱率，W

R——顆粒半徑，mm

Rep,i——顆粒流動Reynolds 數(shù)

St——Stokes 數(shù)

T——溫度，K

Tc,in——冷油介質(zhì)入床溫度，K

Tf——連續(xù)相流體溫度，K

T*f——參考溫度，K

Tp——離散顆粒溫度，K

Δt——迭代時間步長，s

uf——連續(xù)相流體速度矢量

up——離散顆粒速度矢量

vc,in——載流介質(zhì)入床流速

vp——顆粒運動速度，m·s-1

vp,in——顆粒相入口速度，m·s-1

β——流體熱膨脹系數(shù)，m·K-1

μ——動力黏度，kg·m-1·s-1

γ——亞松弛系數(shù)

ε——網(wǎng)格空隙率

λf——連續(xù)相流體熱導率，W·m-1·K-1

ξ——凍結程度關鍵值

ρ——密度，kg·m-3

ρ*f——參考溫度下的流體密度，kg·m-3

τ——時間步長，s

ωp——顆粒轉動的角速度，rad·s-1

?——量綱1 角速度

[1]Zheng Keqing (鄭克晴),Zhang Xuejun (張學軍),Tian Xinjian (田新建),et al.Heat transfer performance of a single air bubble in direct contact ice slurry generator [J].CIESC Journal(化工學報),2010,61 (S2):58-61

[2]Peng Zhengbiao,Yuan Zhulin,Liang Kunfeng,et al.Ice slurry formation in a cocurrent liquid-liquid flow [J].Chinese Journal of Chemical Engineering,2008,16 (4):552-557

[3]Li Hongchang,Li Yaoming,Gao Fang,et al.CFD-DEM simulation of material motion in air-and-screen cleaning device [J].Computers and Electronics in Agriculture,2012,88:111-119

[4]Kevin F Malone,Bao H Xu,Michael Fairweather.The combined-continuum-and-discrete-model (CCDM) for simulation of liquid-particle flows [J].Computer Aided Chemical Engineering,2006,21:227-232

[5]Chu K W,Wang B,Xu D L,et al.CFD-DEM simulation of the gas-solid flow in a cyclone separator [J].Chemical Engineering Science,2011,66 (5):834-847

[6]Liu Xiangjun (劉向軍),Shi Lei (石磊),Xu Xuchang (徐旭常).Activities of dense particle-gas two-phase flow modeling in Eulerian-Lagrangian approach [J].Chinese Journal of Computational Mechanics(計算力學學報),2007,24 (2):166-172

[7]Tsuji Y,Tanaka T.Cluster patterns in circulating fluidized beds predicted by numerical simulation (discrete particle model versus two-fluid model) [J].Powder Technology,1998,95 (3):254-264

[8]Swasdisevi T,Tanthapanichakoon W,Charinpanitkul T,et al.Prediction of gas-particle dynamics and heat transfer in a two-dimensional spouted bed [J].Advanced Powder Technology,2005,16 (3):275-293

[9]Li J T,Mason D J,Mujumdar A S.A numerical study of heat transfer mechanisms in gas-solids flows through pipes using a coupled CFD and DEM model [J].Drying Technology,2003,21:1839-1866

[10]Zhou H,Flamant G,Gauthier D.DEM-LES of coal combustion in a bubbling fluidized bed (Ⅰ):Gas-particle turbulent flow structure [J].Chemical Engineering Science,2004,59 (20):4193-4203

[11]Malone K F,Xu B H.Particle-scale simulation of heat transfer in liquid-fluidized beds [J].Powder Technology,2008,184 (2):189-204

[12]Felice R D.The voidage function for fluid-particle interaction systems [J].International Journal of Multiphase Flow,1994,20 (1):153-159

[13]Oesterle B,Dinh T B.Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers [J].Experiments in Fluids,1998,25:16-22

[14]Tsuji Y,Morikawa Y,Mizuno O.Experimental-measurement of the Magnus force on a rotating sphere at low Reynolds numbers [J].Journal of Fluids Engineering Transactions of the ASME,1985,107:484-488

[15]Cai X L,Wallis G B.A more general cell model for added-mass in two-phase flow [J].Chemical Engineering Science,1994,49 (10):1631-1638

[16]Haid M.Correlations for the prediction of heat transfer to liquid-solid fluidized beds [J].Chemical Engineering and Processing,1997,36 (2):143-147

[17]Bremford D J,Müller-Steinhagen H,Duffy G G.Heat transfer to kraft black liquor in a liquid-solid fluidized bed [J].Heat Transfer Engineering,2000,21:57-69

[18]Agarwal P K.Transport phenomena in multi-particle systems (Ⅱ):Particle-fluid heat and mass transfer [J].Chemical Engineering Science,1988,43 (9):2501-2510

[19]Joseph G G,Zenit R,Hunt M L,et al.Particle-wall collisions in a viscous fluid [J].Journal of Fluid Mechanics,2001,433:329-346

[20]Gondret P,Lance M,Petit L.Bouncing motion of spherical particles in fluids [J].Physics of Fluids,2002,14 (2):643-652

[21]Cai Jie (蔡杰),Yuan Zhulin (袁竹林).Numerical simulation on solids fraction distribution in straw burned circulating fluidized beds [J].Proceedings of the CSEE(中國電機工程學報),2007,27 (20):71-75

[22]Peng Zhengbiao (彭正標).Experiment and modeling of ice slurry production in liquid-liquid calculating fluidization [D].Nanjing:Southeast University,2009