基于群體平衡模型的冰漿流動與傳熱特性數值研究

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

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

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

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

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

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

1 CFD模型

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

1.1 連續性方程

(1)

(2)

1.2 動量守恒方程

固相動量方程如下：

(3)

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

液相動量方程如下：

(4)

式中：τl為液相應力張量.

1.3 能量守恒方程

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

(5)

(6)

式中：H為焓值；λeff為等效熱導率；hV為單位體積內固-液兩相間的傳熱系數；T為溫度.

2 群體平衡模型

2.1 模型描述

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

(7)

聚并是多個顆粒之間的團聚現象.聚并會導致顆粒尺寸增大，其顆粒聚并頻率可用布朗運動和顆粒碰撞作用下的模型描述[13]：

(8)

式中：kB為Boltzmann常數；μ為液體黏度；ε為湍流擴散率；η為動力黏度.

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

δ(V)=c1ηc2εc3Vc4

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

β(V|V′)=

(9)

2.2 CFD與PBM模型的耦合過程

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

3 數值計算

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

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

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

表1 不同條件下的流動與傳熱工況

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

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

4 模型驗證及數值計算結果

4.1 模型驗證

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

4.2 CFD-PBM流動工況結果

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

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

各工況下管道出口截面對稱軸上的顆粒直徑及CFD-PBM1工況出口截面的顆粒分布如圖2所示，其中R為對稱軸上任意位置到軸線的距離與半徑之比.由圖2可知，CFD1模型中的固體顆粒直徑無任何變化，而CFD-PBM模型中近壁面處的顆粒直徑比主流區域的顆粒直徑大.這是由于近壁面處的湍流擴散率較大，導致了聚并效果增強，最終使近壁面處的顆粒直徑比主流區的顆粒直徑大.對比圖2(a)中CFD-PBM1工況和CFD-PBM2工況的結果發現，入口流速越大，顆粒直徑就越較大，這是由于較大的流速有利于顆粒的碰撞，從而強化了顆粒的聚并作用.同時，流速越大，冰漿顆粒的分布越均勻，進而導致管內頂部與底部顆粒直徑大小差異減小.通過上述分析可知，管道內顆粒的最大直徑出現在出口頂部，在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工況下出口截面固相體積分數分布Fig.3 Solid volume fraction distributions of outlet for CFD-PBM1

4.3 CFD-PBM傳熱工況結果

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

圖4 沿流動方向軸線的固相體積分數變化Fig.4 Solid volume fraction along the axis of the flow direction

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

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

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

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

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

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

5 結論

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

(1) 使用CFD-PBM模型能夠更準確地反映出冰漿在水平圓管中的流動與傳熱特性.研究發現，冰漿的固相體積分數與顆粒直徑大小呈正相關，高固相體積分數區域的顆粒直徑比低固相體積分數區域的顆粒直徑大，同時顆粒直徑的增長也導致了固相體積分數的增加.而流速越快會增加聚并作用，最終加速顆粒直徑的增長.

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

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