梁運(yùn)濤，王樹剛，張騰飛，王繼紅

（1.煤炭科學(xué)研究總院沈陽(yáng)研究院 煤炭安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，沈陽(yáng)110016；2.大連理工大學(xué) 土木學(xué)院，大連116024）

區(qū)域供冷在城市用電結(jié)構(gòu)優(yōu)化中扮演重要角色［1］。冰漿流體作為蓄冷介質(zhì)，因具有良好的熱物性、傳輸性及相變特性正日益受到人們廣泛關(guān)注。作為冰漿介質(zhì)工程推廣的基礎(chǔ)性研究，冰漿流體阻力特性成為當(dāng)前各國(guó)學(xué)者的研究熱點(diǎn)。目前所采用的研究方法多以實(shí)驗(yàn)測(cè)試為基礎(chǔ)，結(jié)合描述冰漿流體的Bingham、Herschel－Bulkley、Powder Law 和Casson等非牛頓流變特性［2－5］，建立冰漿流體管道壓降預(yù)測(cè)模型。此類模型雖然能夠在一定范圍內(nèi)較為便捷地預(yù)測(cè)出冰漿流體阻力特性，但適用條件苛刻，難以同時(shí)兼顧輸送管徑、管速、冰粒子的濃度、粒徑、載流體的種類及添加劑的濃度等多因素對(duì)阻力特性的影響。近年來(lái)，固液兩相流計(jì)算流體力學(xué)（CFD）模擬技術(shù)為研究冰漿流體阻力特性帶來(lái)了新的契機(jī)。筆者以冰粒子與載流體所構(gòu)成的固液混合物為研究對(duì)象，采用無(wú)相變過程的兩相流混合模型，運(yùn)用CFD模擬為工具，對(duì)不同工況下漿體所表現(xiàn)出的牛頓與非牛頓流變特性給予分段考慮，將CFD模擬結(jié)果與既有文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù)及其它模型預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析。

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

基于無(wú)相變過程的兩相流混合模型［6－9］，描述其流動(dòng)特性的數(shù)學(xué)模型可表示如下。

1.1 連續(xù)性方程

式中：ρm為混合相流體的密度，kg·m－3；um為混合相流體的速度，m·s－1；▽為拉普拉斯算子。

1.2 動(dòng)量方程

式中：i＝1和2分別為液體相和冰粒子相；p為靜壓力，Pa；τm為剪切應(yīng)力，Pa；αL和αS分別為液體相和冰粒子相的體積濃度，%；uDi為拖曳速度，m·s－1。

1.3 冰粒子相濃度方程

基于無(wú)相變過程的冰粒子相濃度方程，如式（3）所示。

1.4 湍流方程

以求解雷諾平均的納維爾 斯托克斯方程（RANS）為基礎(chǔ)，描述混合相流體湍動(dòng)的k－ε方程，如式（4）和（5）所示。

1.5 封閉條件

1.5.1 拖曳速度模型 多相流中，拖曳速度uDi表示i相與混合相間速度差。據(jù)此定義可得冰粒子相與液體相的拖曳速度，如式（6）和（7）所示［6］。

式中：uLS為液固兩相間的滑移速度，m·s－1，其具體計(jì)算關(guān)系如式（8）和（9）所示［6］。

式中：dS為冰粒子直徑，m；ReS為冰粒子雷諾數(shù)；μt，m為湍流粘度，kg·m－1·s－1；σD為普朗特?cái)U(kuò)散系數(shù)；fdrag為拖曳系數(shù)。

1.5.2 粘度模型 水平管道內(nèi)冰漿流體的輸送過程中，當(dāng)漿體輸送速度較高時(shí)，冰粒子沿管道截面近似呈均勻分布，此時(shí)Thomas方程［10］可很好地表示混合相流體的粘度特性，如式（10）所示。

式中：μm為混合相流體粘度，kg·m－1·s－1；μL為液體相粘度，kg·m－1·s－1。

隨著漿體速度降低，冰粒子沿管道截面梯級(jí)分布，混合相流體非牛頓特性逐漸顯現(xiàn)。本研究發(fā)現(xiàn)，若系統(tǒng)在低速運(yùn)行時(shí)致使管道頂部冰粒子的濃度αS≥30%，此時(shí)Bingham流變模型［2］可以很好地描述出低速工況下冰漿流體的流變特性，如式（11）所示。

式中：μp為混合相流體塑性粘度，kg·m－1·s－1；τp為混合相流體屈服應(yīng)力，Pa；為混合相流體速度變形率，s－1。其中，μp與τp取值見文獻(xiàn)［11］。

采用CFD方法求解上述微分方程組，首先沿著CFD網(wǎng)格單元對(duì)各微分方程實(shí)施積分，然后基于散度定理離散得到相應(yīng)離散方程，最后應(yīng)用SIMPLE算法耦合連續(xù)性與動(dòng)量方程。方程求解的邊界條件類型在管道入口為質(zhì)量流率入口，出口為壓力出口。湍流流動(dòng)時(shí)壁面處應(yīng)用標(biāo)準(zhǔn)壁面函數(shù)。流動(dòng)計(jì)算區(qū)域采用六面體網(wǎng)格單元，平均尺寸為1.7mm。通過對(duì)所獲得的數(shù)值計(jì)算結(jié)果進(jìn)行網(wǎng)格獨(dú)立性檢驗(yàn)，網(wǎng)格數(shù)量符合獨(dú)立性要求。數(shù)值計(jì)算精度滿足前后兩次迭代值相對(duì)誤差小于10－6。

2 算例驗(yàn)證與對(duì)比分析

2.1 算例概況

以文獻(xiàn)［12］所報(bào)道的實(shí)驗(yàn)結(jié)果為參考，利用前文所介紹的兩相流CFD模擬方法研究水平管道內(nèi)冰漿流體等溫流動(dòng)過程。其中漿體輸送管徑D＝0.016m，管長(zhǎng)L＝6m，輸送平均速度um＝0.22～2.02m·s－1，液體相為質(zhì)量濃度CA＝10.6%的乙醇溶液。各相基本物理參數(shù)如表1所示。

表1 流體的基本物理參數(shù)

2.2 數(shù)值結(jié)果分析

圖1給出了不同濃度下水平管道內(nèi)冰漿流體管道壓降與平均輸送速度間變化關(guān)系。總體上漿體的管道壓降與平均輸送速度間呈單調(diào)遞增關(guān)系，但在低速流動(dòng)區(qū)域內(nèi)（本例中，輸送速度小于0.78m·s－1），壓降隨輸送速度的遞變速率略有降低。由圖1（a）可知，當(dāng)采用Thomas粘度方程時(shí)，CFD模擬結(jié)果與實(shí)驗(yàn)值整體上體現(xiàn)出較好的一致性，但隨著漿體速度的降低，CFD預(yù)測(cè)值明顯低于實(shí)驗(yàn)測(cè)試結(jié)果，在某些低速工況點(diǎn)處最大偏差可達(dá)48.9%。這是因?yàn)殡S著漿體速度降低，管道頂部冰粒子濃度升高，此時(shí)冰漿流體逐漸表現(xiàn)出非牛頓流體特性。在所驗(yàn)證的4組工況中，當(dāng)漿體輸送速度低于0.78m·s－1、管道頂部冰粒子濃度高于30%時(shí)，冰粒子沿管道截面分布不均勻性加劇，粒子間相互作用增強(qiáng)，此時(shí)Thomas粘度方程的高階項(xiàng)不足以反映冰粒子間相互作用，以及由此造成的混合相流體粘度特性的變化。但是，借助于Bingham流變模型，則可較為理想地刻畫出低速工況下水平管道內(nèi)冰漿流體的非牛頓特性。基于Bingham流變特性并采用兩相流無(wú)滑移混合模型，對(duì)圖1（a）中虛線內(nèi)工況點(diǎn)進(jìn)行修正處理。由圖1（b）可知，在修正后低速工況點(diǎn)處，CFD模擬獲得了與實(shí)驗(yàn)測(cè)試非常一致的結(jié)果，最大偏差不足15%。

圖1 水平管道內(nèi)不同冰漿濃度下冰漿流體的管道壓降

2.3 預(yù)測(cè)模型對(duì)比

冰漿流體管道壓降受多因素制約，而現(xiàn)行的基于宏觀實(shí)驗(yàn)現(xiàn)象所建立起的阻力特性預(yù)測(cè)模型往往難以很好地平衡模型的精度與通用性間關(guān)系。這是因?yàn)榇祟惸Ｐ投嗍墙柚诒鶟{流體非牛頓流變特性，再結(jié)合單相流體的Darcy阻力關(guān)系建立的。在建模過程中，由于實(shí)驗(yàn)條件所限，難以綜合考慮多類影響因素。以圖1中所描述的各工況點(diǎn)為例，分別選取文獻(xiàn)［13］、［14］及［15］所建立的3類具有代表性的管道壓降預(yù)測(cè)模型進(jìn)行計(jì)算，其詳細(xì)的數(shù)學(xué)描述如表2所示。

圖2給出了各類預(yù)測(cè)模型的預(yù)測(cè)結(jié)果與實(shí)驗(yàn)對(duì)比。由圖2中（a）、（b）及（c）不難發(fā)現(xiàn)，文獻(xiàn)［13］所介紹的壓降預(yù)測(cè)模型雖能較好地預(yù)測(cè)出大部分工況點(diǎn)的壓降值，但在個(gè)別工況點(diǎn)處模型預(yù)測(cè)精度并不穩(wěn)定，預(yù)測(cè)效果的階躍性及隨機(jī)性較大。而文獻(xiàn)［14］和［15］所建立的預(yù)測(cè)模型僅在高速、低濃度工況下體現(xiàn)出較好的預(yù)測(cè)效果，當(dāng)漿體濃度升高、輸送速度降低時(shí)，模型預(yù)測(cè)值與實(shí)驗(yàn)值間產(chǎn)生較大偏差，從而大大削弱了預(yù)測(cè)模型的實(shí)用價(jià)值。圖2（d）為所采用的CFD模型預(yù)測(cè)結(jié)果，由對(duì)比可知，模擬值與實(shí)驗(yàn)值間的相對(duì)誤差可控制在±15%內(nèi)，模擬結(jié) 果與實(shí)驗(yàn)測(cè)試結(jié)果吻合良好。

表2 冰漿流體管道壓降預(yù)測(cè)模型數(shù)學(xué)描述

圖2 水平管道內(nèi)冰漿流體管道壓降預(yù)測(cè)模型的預(yù)測(cè)值與實(shí)驗(yàn)值對(duì)比

3 結(jié) 論

采用兩相流混合模型，對(duì)不同工況下漿體所呈現(xiàn)的牛頓與非牛頓流變特性給予分段考慮，通過CFD模擬方法研究了水平管道內(nèi)冰漿流體的等溫流動(dòng)，結(jié)論如下。

1）在水平管道內(nèi)冰漿流體輸送過程中，當(dāng)漿體速度較高時(shí)，冰粒子沿管道截面近似呈均勻分布，此時(shí)Thomas方程可很好地描述混合相流體的粘度特性。隨著漿體速度降低，管道頂部冰粒子濃度升高，冰漿流體逐漸表現(xiàn)出非牛頓流體特性，此時(shí)利用Bingham流體模型來(lái)描述低速工況下混合相流體的流變特性所獲得的效果較為理想。

2）同基于冰漿流體非牛頓流變特性管道壓降預(yù)測(cè)模型相比，筆者所采用的數(shù)值計(jì)算模型因全面考慮到漿體輸送管徑、管速、冰粒子的濃度、粒徑、載流體的種類及添加劑的濃度等多因素影響，模型的精度與通用性可得到較好平衡。在所驗(yàn)證的4組工況中，采用的數(shù)值模型的預(yù)測(cè)值與實(shí)驗(yàn)值之間的相對(duì)誤差均控制在±15%以內(nèi)。

［1］魏琪.等熱流圓管內(nèi)冰漿流體傳熱特性的解析研究［J］.蘇州大學(xué)學(xué)報(bào)，2006，26（3）：1－3.WEI Qi.Analytical study on heat transfer of ice slurry in a circular tube with constant heat flux［J］.Journal of Soochow University：Engineering Science，2006，26（3）：1－3.

［2］Grozdek M， Khodabandeh R， Lundqvist P.Experimental investigation of ice slurry flow pressure drop in horizontal tubes［J］.Experimental Thermal and Fluid Science，2009，33（2）：357－370.

［3］Illan F，Viedma A.Experimental study on pressure drop and heat transfer in pipelines for brine based ice slurry Part II：dimensional analysis and rheological model ［J］.International Journal of Refrigeration，2009，32（5）：1024－1031.

［4］Ayel V，Lottin O，Peerhossaini H.Rheology，flow behaviour and heat transfer of ice slurries：a review of the state of the art ［J］.International Journal of Refrigeration，2003，26（1）：95－107.

［5］Kitanovski A，Vuarnoz D，Caesar D A，et al.The fluid dynamics of ice slurry ［J］.International Journal of Refrigeration，2005，28（1）：37－50.

［6］Ling J，Skudarnov P V，Lin C X，et al.Numerical investigations of liquid－solid slurry flows in a fully developed turbulent flow region ［J］.International Journal of Heat and Fluid Flow，2003，24（3）：389－398.

［7］Lin C X，Ebadian M A.A numerical study of developing slurry flow in the entrance region of a horizontal pipe［J］.Computers ＆Fluids，2008，37（8）：965－974.

［8］Skudarnov P V，Lin C X，Ebadian M A.Doublespecies slurry flow in a horizontal pipeline［J］.Journal of Fluids Engineering，2004，126（1）：125－132.

［9］張宏兵，陳露露，謝榮華.水平圓管固液兩相穩(wěn)態(tài)流動(dòng)流動(dòng)特性數(shù)值模擬［J］.化工學(xué)報(bào)，2009，60（5）：1162－1168.ZHANG Hongbing，CHEN Lulu，XIE Ronghua.Numerical simulation of solid－liquid two－phase steady flow in horizontal pipe ［J］.Journal of Chemical Industry and Engineering，2009，60（5）：1162－1168.

［10］Thomas D G.Transport characteristics of suspension［J］.Journal of Colloid Science，1965，20（3）：267－277.

［11］Beata N Z，Zelasko J.Generalized non－Newtonian flow of ice－slurry［J］.Chemical Engineering and Processing，2007，46（10）：895－904.

［12］Beata N Z，Wojciech Z.Momentum transfer of ice slurry flows in tubes，experimental investigations［J］.International Journal of Refrigeration，2006，29（3）：418－428.

［13］Tomita Y.On the fundamental of non－Newtonian flow［J］.Bulletin of the Japanese Society of Mechanical Engineers，1959，2（7）：469－474.

［14］Kauffeld M，Kawaji M，Egolf P W.Handbook on ice slurries－fundamentals and engineering ［M］.France：International Institute of Refrigeration，2005.

［15］劉永紅.冰漿流動(dòng)特性和傳熱性能的研究［D］.上海：同濟(jì)大學(xué)，1997.