基于群體平衡模型的冰漿流動(dòng)與傳熱特性數(shù)值研究

鄒良旭，馬 非，孟昭男，張 鵬

(上海交通大學(xué) 制冷與低溫工程研究所，上海 200240)

由于全球化石能源的不斷消耗以及溫室氣體的排放使得各種環(huán)境問題日益凸顯，人們愈加重視節(jié)能減排和提高能源利用率.蓄能技術(shù)可以解決能量供需在時(shí)間和空間上不均衡的問題，可以有效地提高能源利用率.冰漿是一種固-液兩相可流動(dòng)的相變材料，其固體顆粒直徑通常不大于1.0 mm[1]，并具有較好的儲(chǔ)能、流動(dòng)和傳熱特性，廣泛應(yīng)用于建筑空調(diào)、區(qū)域供冷和食品存儲(chǔ)等領(lǐng)域[2-4].在實(shí)際應(yīng)用中，冰漿固相顆粒直徑在流動(dòng)和傳熱過程中不斷變化，且冰漿的流動(dòng)特性和固液兩相之間的傳熱傳質(zhì)過程十分復(fù)雜，因此深入理解冰漿的流動(dòng)與傳熱特性對(duì)其應(yīng)用十分重要.

目前所進(jìn)行的研究大多忽略了冰漿固相顆粒直徑在流動(dòng)與傳熱過程中的變化，采用數(shù)值計(jì)算獲得冰漿的流動(dòng)與傳熱特性.Wang等[5]通過對(duì)比實(shí)驗(yàn)與數(shù)值模擬結(jié)果發(fā)現(xiàn)，使用基于固相顆粒動(dòng)力學(xué)理論的 Euler-Euler方法可以更準(zhǔn)確地模擬出冰漿的流動(dòng)過程，以及彎管處的阻力對(duì)冰漿流動(dòng)產(chǎn)生的影響.Kumano等[6]研究了冰漿存儲(chǔ)時(shí)間對(duì)冰漿流動(dòng)特性的影響，作者利用顯微鏡對(duì)冰漿顆粒進(jìn)行觀察，發(fā)現(xiàn)流動(dòng)工況下其直徑會(huì)隨時(shí)間的增長而變大，當(dāng)儲(chǔ)存時(shí)間超過12小時(shí)，圓管內(nèi)因冰顆粒的聚集而出現(xiàn)堵塞現(xiàn)象.

傳統(tǒng)的兩相流數(shù)值模擬計(jì)算模型忽略了冰漿固相顆粒直徑的變化.群體平衡模型是一種針對(duì)顆粒直徑變化的數(shù)值計(jì)算模型，可較為準(zhǔn)確地描述冰漿內(nèi)固相顆粒的破碎、融化和聚并等過程.Xu等[7]采用基于Euler-Euler方法的計(jì)算流體動(dòng)力學(xué)-群體平衡模型(CFD-PBM)對(duì)冰漿在水平直管內(nèi)的流動(dòng)進(jìn)行了數(shù)值計(jì)算，計(jì)算獲得的壓降與實(shí)驗(yàn)數(shù)據(jù)較為吻合，驗(yàn)證了CFD-PBM模型的可靠性及準(zhǔn)確性，但在該文獻(xiàn)中尚未涉及冰漿的傳熱特性.Jin等[8]利用兩相流模型和群體平衡模型對(duì)氮漿在水平管道中的流動(dòng)與傳熱特性進(jìn)行了模擬計(jì)算，研究結(jié)果表明加入群體平衡模型后的計(jì)算模型更準(zhǔn)確，在不同工況下壓降和傳熱系數(shù)的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的偏差在5%以內(nèi)，但該文獻(xiàn)未對(duì)顆粒直徑變化對(duì)流動(dòng)與傳熱特性的影響進(jìn)行討論.

針對(duì)上述問題，本文采用基于Euler-Euler法的CFD-PBM模型，對(duì)冰漿在水平圓管中的流動(dòng)與傳熱特性及其固相顆粒的直徑變化進(jìn)行研究.在流動(dòng)特性的研究中，分析固相體積分?jǐn)?shù)與流速對(duì)顆粒直徑在流動(dòng)過程中的影響，為傳熱特性的研究提供依據(jù)；在傳熱特性的研究中，分析不同熱流密度對(duì)顆粒直徑變化的影響.同時(shí)在上述兩類工況中，將CFD模型與CFD-PBM模型的計(jì)算結(jié)果進(jìn)行對(duì)比，因CFD-PBM模型考慮了冰漿顆粒的直徑變化，使得CFD-PBM模型對(duì)冰漿流動(dòng)與傳熱過程的數(shù)值模擬更接近實(shí)際情況.

1 CFD模型

采用基于Euler-Euler法的CFD模型對(duì)冰漿在水平圓管中的流動(dòng)與傳熱特性進(jìn)行模擬計(jì)算.假設(shè)固-液兩相均為連續(xù)性流體，且均具有質(zhì)量、動(dòng)量和熱量的交換；冰漿在管內(nèi)流動(dòng)為不可壓縮流體湍流動(dòng)，固相為表面光滑、無彈性的球體[9].采用標(biāo)準(zhǔn)k-ε湍流模型描述冰漿在水平圓管中的湍流過程[10].因冰漿的兩相密度差異較小，故相間作用力只考慮曳力、升力和湍流擴(kuò)散力，顆粒動(dòng)力學(xué)模型詳見文獻(xiàn)[9].

1.1 連續(xù)性方程

(1)

(2)

1.2 動(dòng)量守恒方程

固相動(dòng)量方程如下：

(3)

式中：p為局部壓力；g為重力加速度；F為兩相之間的總作用力；ps為固相壓力；τs為固相應(yīng)力張量.

液相動(dòng)量方程如下：

(4)

式中：τl為液相應(yīng)力張量.

1.3 能量守恒方程

在傳熱的工況下需要考慮能量守恒方程，基于Euler-Euler模型的假設(shè)，冰漿中的固-液兩相都遵循能量守恒方程：

(5)

(6)

式中：H為焓值；λeff為等效熱導(dǎo)率；hV為單位體積內(nèi)固-液兩相間的傳熱系數(shù)；T為溫度.

2 群體平衡模型

2.1 模型描述

在冰漿的固-液兩相流中，固相顆粒在實(shí)際物理過程中存在破碎、融化和聚并等現(xiàn)象，這些現(xiàn)象會(huì)導(dǎo)致固相顆粒直徑分布和尺寸大小發(fā)生變化.群體平衡模型可用來描述固體顆粒在各種現(xiàn)象下的分布狀態(tài)，使用正交矩方法對(duì)群體平衡方程進(jìn)行求解并計(jì)算其Sauter平均直徑，即固相顆粒直徑[11-12].關(guān)于固相顆粒的群體平衡方程如下：

(7)

聚并是多個(gè)顆粒之間的團(tuán)聚現(xiàn)象.聚并會(huì)導(dǎo)致顆粒尺寸增大，其顆粒聚并頻率可用布朗運(yùn)動(dòng)和顆粒碰撞作用下的模型描述[13]：

(8)

式中：kB為Boltzmann常數(shù)；μ為液體黏度；ε為湍流擴(kuò)散率；η為動(dòng)力黏度.

破碎的作用與聚并相反，它表示顆粒之間由于碰撞導(dǎo)致的顆粒尺寸減小的過程.破碎現(xiàn)象會(huì)受到冰漿固相顆粒直徑尺寸的影響，顆粒破碎頻率由下式表示[7,14]：

δ(V)=c1ηc2εc3Vc4

其中：c1～c4均為系數(shù)，其值大小分別為 0.000 6，-1.25，0.75，1.在破碎模型中，固相顆粒破碎的概率密度函數(shù)描述了由破碎導(dǎo)致顆粒的體積由V′變?yōu)閂的概率，并且可以描述顆粒數(shù)量以及由破碎導(dǎo)致的尺寸分布.本文的固相顆粒為球狀顆粒，不同顆粒的大小不同，因此采用拋物型概率密度函數(shù)，其表達(dá)式為

β(V|V′)=

(9)

2.2 CFD與PBM模型的耦合過程

CFD-PBM模型首先利用CFD兩相流模型，通過求解質(zhì)量守恒、動(dòng)量守恒、能量守恒等方程來獲得固-液兩相的流速、溫度以及固相體積分?jǐn)?shù)；再采用PBM模型通過正交矩方法求解固相顆粒的聚并、破碎以及融化過程，并且得到顆粒的尺寸及其分布；最后通過獲得的固相顆粒直徑，重新計(jì)算固相的流速、溫度以及體積分?jǐn)?shù)，如此循環(huán)以實(shí)現(xiàn)兩種模型的耦合[15].

3 數(shù)值計(jì)算

考慮兩類計(jì)算工況，一種僅計(jì)算冰漿在非均質(zhì)狀態(tài)下的等溫流動(dòng)工況，稱為流動(dòng)工況；另一種計(jì)算管道壁面有熱流時(shí)冰漿的非均質(zhì)、非等溫流動(dòng)的情況，稱為傳熱工況.上述兩類工況均采用質(zhì)量分?jǐn)?shù)為10.3%的乙醇-水溶液作為冰漿的載流體.在數(shù)值計(jì)算中，對(duì)于壁面邊界條件，液相采用無滑移壁面條件，固相采用Johnson-Jackson壁面滑移條件[9].在傳熱工況中，還需在管壁添加熱流密度邊界，假定圓管入口處的固-液兩相為相同流速和溫度.

在流動(dòng)工況中，圓管的長度為 1.3 m，直徑為23 mm.工質(zhì)的物性考慮為常數(shù)[9].固體顆粒的初始直徑為125 μm.在傳熱工況中，為了觀察到明顯的冰漿融化過程，采用長度為 2.0 m的圓管，直徑為16 mm，固-液兩相的入口流速為1.0 m/s，固相體積分?jǐn)?shù)為10%.在傳熱工況下，工質(zhì)的物性隨溫度發(fā)生變化[16]，固相顆粒的初始直徑設(shè)置為270 μm.除了考慮工況，還需對(duì)比CFD模型與CFD-PMB模型在相同工況下的結(jié)果，數(shù)值模擬計(jì)算的工況及條件如表1所示.其中：CFD表示不加PBM模型的模擬計(jì)算；CFD-PBM表示加入PBM模型的模擬計(jì)算；工況名后的數(shù)字為編號(hào).

利用FLUENT 14.5軟件進(jìn)行數(shù)值模擬計(jì)算，計(jì)算時(shí)間步長為0.01 s.圓管網(wǎng)格劃分采用六面體結(jié)構(gòu)化網(wǎng)格，在近壁面處加密處理，網(wǎng)格均進(jìn)行網(wǎng)格無關(guān)性分析以保證計(jì)算的準(zhǔn)確性.流動(dòng)工況下的網(wǎng) 格數(shù)量為 1 347 312，傳熱工況下的網(wǎng)格數(shù)量為2 739 380.在數(shù)值計(jì)算過程中，所有的控制方程均采用一階迎風(fēng)離散格式進(jìn)行計(jì)算.每個(gè)時(shí)間步長的最大迭代次數(shù)為50次，當(dāng)每一步迭代計(jì)算的速度項(xiàng)殘差低于1.0×10-4、能量項(xiàng)殘差低于1.0×10-6時(shí)，視為計(jì)算收斂并進(jìn)行下一步計(jì)算.

表1 不同條件下的流動(dòng)與傳熱工況

Tab.1 Flow and heat transfer of ice slurry under different conditions

工況入口流速/(m·s-1)固相體積分?jǐn)?shù)/%壁面熱流密度/(kW·m-2)CFD11.010-CFD-PBM11.010-CFD-PBM22.010-CFD-PBM31.015-CFD21.01050CFD-PBM41.01015CFD-PBM51.01025CFD-PBM61.01050

4 模型驗(yàn)證及數(shù)值計(jì)算結(jié)果

4.1 模型驗(yàn)證

CFD-PBM模型對(duì)冰漿流動(dòng)的數(shù)值模擬計(jì)算的準(zhǔn)確性驗(yàn)證見文獻(xiàn)[7].由于冰漿固相顆粒直徑的測量在實(shí)驗(yàn)中難以捕捉測量，所以采用壓降作為模型準(zhǔn)確性的衡量手段.在其實(shí)驗(yàn)中，冰漿載流體為質(zhì)量分?jǐn)?shù)6.5%的乙二醇-水溶液，在不同流動(dòng)工況下，實(shí)驗(yàn)得到的冰漿在水平圓管中流動(dòng)的壓降與CFD-PBM模型的計(jì)算結(jié)果最大偏差為10.32%，驗(yàn)證了CFD-PBM計(jì)算冰漿流動(dòng)工況的準(zhǔn)確性.

4.2 CFD-PBM流動(dòng)工況結(jié)果

在流動(dòng)工況中不涉及到熱量的傳遞，冰漿中固相顆粒之間只有聚并和破碎作用，不存在融化作用.固相顆粒直徑D沿管道軸線的分布如圖1所示，其中Lpip為管道軸線距離.由圖1可知，CFD模型固相顆粒直徑保持不變，而CFD-PBM模型的固相顆粒直徑沿流動(dòng)方向逐漸增大.分析CFD-PBM模型的計(jì)算結(jié)果可知，在Lpip=0～0.1 m區(qū)間內(nèi)，固相顆粒的聚并和破碎效果很弱，因此顆粒直徑幾乎沒有任何變化；Lpip=0.1～0.6 m區(qū)間是顆粒直徑增長的發(fā)展段，聚并作用比破碎作用強(qiáng)，顆粒直徑開始增長；Lpip=0.6～1.3 m區(qū)間顆粒直徑以穩(wěn)定的增長速率逐漸增大.在CFD-PBM1、CFD-PBM2、CFD-PBM3工況中，沿管道流動(dòng)方向，顆粒直徑分別從初始狀態(tài)的125.0 μm增長到127.6、127.8和128.9 μm，說明增大固相體積分?jǐn)?shù)及流速均可以促進(jìn)顆粒直徑的增長.

圖1 顆粒直徑沿管道軸線的分布Fig.1 Mean diameter distributions along central axial direction

各工況下管道出口截面對(duì)稱軸上的顆粒直徑及CFD-PBM1工況出口截面的顆粒分布如圖2所示，其中R為對(duì)稱軸上任意位置到軸線的距離與半徑之比.由圖2可知，CFD1模型中的固體顆粒直徑無任何變化，而CFD-PBM模型中近壁面處的顆粒直徑比主流區(qū)域的顆粒直徑大.這是由于近壁面處的湍流擴(kuò)散率較大，導(dǎo)致了聚并效果增強(qiáng)，最終使近壁面處的顆粒直徑比主流區(qū)的顆粒直徑大.對(duì)比圖2(a)中CFD-PBM1工況和CFD-PBM2工況的結(jié)果發(fā)現(xiàn)，入口流速越大，顆粒直徑就越較大，這是由于較大的流速有利于顆粒的碰撞，從而強(qiáng)化了顆粒的聚并作用.同時(shí)，流速越大，冰漿顆粒的分布越均勻，進(jìn)而導(dǎo)致管內(nèi)頂部與底部顆粒直徑大小差異減小.通過上述分析可知，管道內(nèi)顆粒的最大直徑出現(xiàn)在出口頂部，在CFD-PBM1、CFD-PBM2和CFD-PBM33種工況的最大顆粒直徑分別為139，142，146 μm.

圖2 管道出口截面顆粒直徑的分布Fig.2 Mean diameter distributions of ice slurry on the cross-section at the pipe outlet

圖3 CFD-PBM1工況下出口截面固相體積分?jǐn)?shù)分布Fig.3 Solid volume fraction distributions of outlet for CFD-PBM1

4.3 CFD-PBM傳熱工況結(jié)果

由于傳熱工況中熱量的傳遞會(huì)導(dǎo)致冰漿固相顆粒融化，而融化會(huì)使冰漿固相體積分?jǐn)?shù)降低、顆粒直徑減小，因此冰漿中固體顆粒之間同時(shí)具有聚并、破碎和融化作用.不同工況下冰漿沿管道軸線的固相體積分?jǐn)?shù)如圖4所示.在冰漿的傳熱中，由于管壁的熱量先傳遞到液相，使液相溫度升高，隨后液相將熱量再傳遞給固相，使固相顆粒隨著流動(dòng)不斷融化，固相體積分?jǐn)?shù)不斷降低.在CFD2、CFD-PBM4、CFD-PBM5及CFD-PBM6工況中，固相體積分?jǐn)?shù)分別從初始狀態(tài)的10%下降至出口截面處的1.21%、6.37%、4.33%及0.99%.

圖4 沿流動(dòng)方向軸線的固相體積分?jǐn)?shù)變化Fig.4 Solid volume fraction along the axis of the flow direction

不同工況下，冰漿沿流動(dòng)方向軸線的固相顆粒直徑變化如圖5所示.CFD模型下固相顆粒的直徑保持不變，而CFD-PBM模型的固相顆粒直徑沿著流動(dòng)方向的軸線距離逐漸減小.在CFD-PBM4、CFD-PBM5和CFD-PBM6工況中，沿管道流動(dòng)的軸線方向，顆粒直徑分別從初始狀態(tài)的270 μm融化至255、249和154 μm.結(jié)合圖4和圖5可知，在熱流密度較大的條件下，固相顆粒直徑減小的速度更快.顆粒直徑的減小有利于冰漿的融化，使得冰漿固相體積分?jǐn)?shù)減小得更多，因此隨著顆粒直徑的減小，CFD-PBM6工況的固相體積分?jǐn)?shù)比同條件下CFD2工況的固相體積分?jǐn)?shù)略低.

圖5 沿流動(dòng)方向軸線的固相顆粒直徑變化Fig.5 Mean diameter of ice slurry along the axis of the flow direction

CFD-PBM6工況下的固相顆粒平均直徑及其分布如圖6所示.由圖6可知，顆粒直徑沿管道流動(dòng)方向逐漸減小.對(duì)于管道不同位置的截面對(duì)稱軸，顆粒直徑分布如圖7所示.由圖7可知，由于固相體積分?jǐn)?shù)與顆粒直徑的相互影響，管道頂部固相顆粒直徑比底部大.隨著沿程距離的增大，管道近壁面處的顆粒完全融化，管道頂部與底部顆粒直徑的差異逐漸減小，故出口附近的顆粒直徑基本呈對(duì)稱分布.

綜上所述，在傳熱工況下，CFD模型并未考慮管道中的固相顆粒直徑的變化及其分布規(guī)律，與實(shí)際情況不符.CFD-PBM模型同時(shí)考慮了顆粒的聚并、破碎和融化作用(以融化模型為主導(dǎo)作用)，最終導(dǎo)致了顆粒直徑近壁面處的直徑小于主流區(qū)，沿流動(dòng)方向中軸線的固相顆粒直徑逐漸減小.因此，CFD-PBM模型比CFD模型更能合理地反映出冰漿的實(shí)際傳熱特征與變化規(guī)律.

圖6 CFD-PBM6工況下顆粒直徑在管道軸向截面的分布Fig.6 Mean diameter distributions on the vertical symmetrical cross-section of the pipe for CFD-PBM6

圖7 CFD-PBM6工況下顆粒直徑在不同位置的截面對(duì)稱軸上的分布Fig.7 Mean diameter distributions at different cross-sections along the flow direction for CFD-PBM6

5 結(jié)論

本文采用基于Euler-Euler法的CFD-PBM模型，對(duì)冰漿在水平圓管中的流動(dòng)與傳熱特性及其固相顆粒直徑的變化規(guī)律進(jìn)行了數(shù)值模擬計(jì)算.通過對(duì)比不同條件下的CFD模型與CFD-PBM模型的計(jì)算結(jié)果可以獲得以下結(jié)論：

(1) 使用CFD-PBM模型能夠更準(zhǔn)確地反映出冰漿在水平圓管中的流動(dòng)與傳熱特性.研究發(fā)現(xiàn)，冰漿的固相體積分?jǐn)?shù)與顆粒直徑大小呈正相關(guān)，高固相體積分?jǐn)?shù)區(qū)域的顆粒直徑比低固相體積分?jǐn)?shù)區(qū)域的顆粒直徑大，同時(shí)顆粒直徑的增長也導(dǎo)致了固相體積分?jǐn)?shù)的增加.而流速越快會(huì)增加聚并作用，最終加速顆粒直徑的增長.

(2) 在流動(dòng)工況中，隨著在圓管中的流動(dòng)，固相顆粒直徑逐漸增大.在入口流速為1.0 m/s與固相體積分?jǐn)?shù)為10%的情況下，固相顆粒的直徑最大可從初始狀態(tài)的125 μm增長到139 μm.其中，管道頂部的顆粒直徑大于底部，近壁面處的顆粒直徑大于主流區(qū).

(3) 在傳熱工況中，固相顆粒直徑隨著固相的融化不斷減小.在入口流速為1.0 m/s、固相體積分?jǐn)?shù)為10%與壁面熱流密度為50 kW/m2的工況下，固相顆粒直徑從初始狀態(tài)的270 μm直至完全融化.其中，管道頂部的顆粒直徑大于底部，近壁面處的顆粒直徑小于主流區(qū).出口處附近顆粒直徑呈對(duì)稱分布.熱流密度的增大會(huì)加速固相的融化，進(jìn)而加快顆粒直徑的減小.