張雪，劉圣春，宋麗瑩，徐智明

（天津商業(yè)大學(xué)機(jī)械工程學(xué)院，天津市制冷技術(shù)重點(diǎn)實(shí)驗(yàn)室，天津 300134）

0 引言

隨著全球電力需求的增長(zhǎng)，許多環(huán)境問(wèn)題如全球變暖、臭氧層損耗等受到了廣泛關(guān)注，并嚴(yán)重影響全球能源政策[1]。為了緩解能源緊張，提高能源利用效率，儲(chǔ)能技術(shù)發(fā)揮著關(guān)鍵作用。冰漿是由冰晶、液態(tài)水和防凍劑組成的絮狀混合物，因其具有良好的流動(dòng)特性和較大的能量密度，在同一運(yùn)行工況下相較于傳統(tǒng)冷凍水冷卻能力高4 倍至5 倍[2-3]，已成為一種新型相變材料在蓄冷系統(tǒng)中被廣泛使用[4-8]。

目前，國(guó)內(nèi)外學(xué)者對(duì)冰漿系統(tǒng)的研究較多。在換熱特性方面，JEAN-PIERRE 等[9]實(shí)驗(yàn)研究了雙管換熱器中冰漿的換熱情況，結(jié)果表明，隨著含冰率的減小，傳熱系數(shù)先減小，在低含冰率時(shí)基本保持不變。趙騰磊等[10]以含冰率在20%以下的冰漿為模型進(jìn)行實(shí)驗(yàn)研究，分析冰漿的密度、導(dǎo)熱系數(shù)、比熱和動(dòng)力黏度，結(jié)果表明冰漿的傳熱系數(shù)與含冰率成正比。張曼等[11]建立了動(dòng)態(tài)冰漿傳熱特性的數(shù)理模型，數(shù)值研究了流速、管徑和冰漿體積分?jǐn)?shù)對(duì)傳熱系數(shù)的影響。LEE 等[12]實(shí)驗(yàn)研究了以冰漿作為冷卻介質(zhì)的雙管換熱器，結(jié)果表明換熱率隨質(zhì)量流量和含冰率的增加而增加，流速低時(shí)，冰粒對(duì)換熱的影響更為明顯。BELLAS 等[13]實(shí)驗(yàn)研究了帶冰漿的板式換熱器的性能，得出當(dāng)含冰率在5%～20%時(shí)，傳熱系數(shù)基本保持不變。LIU 等[14]數(shù)值研究了水平直管直徑、長(zhǎng)度、含冰率和流速對(duì)冰漿流動(dòng)壓降的影響，驗(yàn)證了仿真模型的可靠性，進(jìn)一步發(fā)展了黏度與流速的數(shù)學(xué)模型。簡(jiǎn)夕忠等[15]數(shù)值研究了含冰率、流速、管徑和添加劑種類等因素對(duì)冰漿換熱效果的影響。由于冰漿系統(tǒng)涉及到相對(duì)比較復(fù)雜的流動(dòng)及傳熱問(wèn)題，目前，對(duì)冰漿傳熱特性的研究在理論分析、數(shù)值計(jì)算和推廣應(yīng)用方面仍需完善。

本文利用顆粒動(dòng)力學(xué)理論，通過(guò)數(shù)值模擬的方法得出冰漿在不同類型管道內(nèi)的流動(dòng)換熱特性，并對(duì)直管內(nèi)冰堵現(xiàn)象進(jìn)行了分析。

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

在Fluent 軟件中采用雙歐拉流體模型，結(jié)合顆粒動(dòng)力學(xué)理論建立數(shù)值模型的控制方程。

1.1 質(zhì)量守恒方程

冰漿本質(zhì)上為一種液固兩相混合流體，冰粒固體為直徑較小的顆粒，內(nèi)置顆粒相模型，設(shè)兩相分別為q相和p相。

相位q的體積：

q相的有效密度按式(3)計(jì)算：

q相的連續(xù)性方程為：

由質(zhì)量守恒方程，得：

式中，αq為q相的體積分?jǐn)?shù)；ρq為q相的密度，kg/m3；mpq為p相到q相的質(zhì)量傳遞速率，kg/(m3·s)；vq為q相的速度，m/s。

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

將冰漿看作固液兩相流體，用液相動(dòng)力參數(shù)描述動(dòng)量方程[16]，可得：

式中，下標(biāo)s 表示固相；l表示液相；vsl為固液兩相間的相對(duì)速度，若msl＞0，則vsl=vs；若msl＜0，則vsl=vl，且vsl=vls；τl為液相切應(yīng)力，Pa；Fl為外部體積力，N；Flift,l為升力，N；FVm,l為虛擬質(zhì)量力，N；Ksl為固液相間的動(dòng)量交換系數(shù)。

固體顆粒相動(dòng)量方程為：

液相與固體顆粒相之間的作用力可采用GIDASPOW[17]模型確定。

當(dāng)αl＞0.8 時(shí)，固-液動(dòng)量交換系數(shù)Ksl為：

當(dāng)αl≤0.8 時(shí)，Ksl有如下形式：

式中，Re為相對(duì)雷諾數(shù)；ds為顆粒直徑，m。

q相和p相的相對(duì)雷諾數(shù)為：

1.3 能量方程

能量控制方程[18]如下：

式中，λeff為冰漿的有效導(dǎo)熱系數(shù)，W/(m·K)；ST為源項(xiàng)，W/m3。

1.4 湍流方程

選擇k-ε湍流模型，考慮固液兩相間的相互作用[16]：

式中，C1ε=1.44；C2ε=1.92；σε=1.3；σk=1；μt,m=ρmCμk2/ε；Cμ=0.09；下標(biāo)m為固液混合相；k為湍流脈動(dòng)動(dòng)能，J；ε為湍流脈動(dòng)動(dòng)能耗散率；μt,m為混合相湍流動(dòng)力黏度，kg/(m·s)；σk為湍流脈動(dòng)動(dòng)能普朗特?cái)?shù)。

1.5 黏度模型

冰漿是一種固液兩相混合物，物理密度、動(dòng)力黏度以及結(jié)冰點(diǎn)不僅與受添加物體積分?jǐn)?shù)影響，還與含冰率有關(guān)，THOMAS[19]通過(guò)反復(fù)實(shí)驗(yàn)和理論推算，得出黏度關(guān)于含冰率的公式：

式中，φi為含冰率，%。

1.6 相變模型

在冰漿換熱特性模擬過(guò)程中，開(kāi)啟能量方程，將冰-水換熱模擬寫(xiě)入相變UDF 程序中，收斂精度設(shè)為10-4。

2 數(shù)值模擬模型圖及邊界條件設(shè)置

2.1 管道模型

90°彎管和T 型管模型分別如圖1(a)和圖1(b)所示。管道直徑為0.02 m，管總長(zhǎng)為1 m，彎管和T 型管均在中間位置處發(fā)生改變，彎管曲率半徑為0.1 m，網(wǎng)格質(zhì)量在0.7 以上，管道分流處采用結(jié)構(gòu)化網(wǎng)格形式。溶液進(jìn)出口位置及重力方向見(jiàn)圖1。

圖1 90°彎管和T 型管管道網(wǎng)格劃分模型

2.2 邊界條件設(shè)置

對(duì)冰漿在90°彎管、T 型管內(nèi)的換熱特性進(jìn)行模擬。具體參數(shù)設(shè)置為：含冰率IPF 為5%～25%，流速為0.2～1.5 m/s。冰粒直徑約為10-4m。邊界條件設(shè)為速度入口，出口設(shè)為出流，壁面設(shè)為無(wú)滑移，采用定壁溫條件，溫度為1 ℃；入口處冰相和水相的溫度均為0 ℃；湍流模型選用k-ε模型。冰漿兩相流體的物性參數(shù)如表1所示。

表1 冰漿兩相流體的物性參數(shù)[20]

3 數(shù)值模擬結(jié)果分析

3.1 90°彎管內(nèi)冰漿換熱特性結(jié)果分析

圖2 和圖3所示為含冰率為10%時(shí)，90°彎管截面上溶液平均溫度分布和焓值分布。由圖2 和圖3 可知，縱截面上，溶液平均溫度變化較小，保持在0～1 ℃，彎管處溶液焓值最大為2.46×104J/kg。

圖2 90°彎管截面上溶液平均溫度分布

圖3 90°彎管截面上溶液的焓值分布

圖4所示為90°彎管截面上流體流動(dòng)的速度分布。由圖4 可知，在90°彎管中，流體流動(dòng)速度從入口處逐漸增大直至彎管處達(dá)到最大，約為1 m/s，經(jīng)過(guò)轉(zhuǎn)彎處，出現(xiàn)明顯分層現(xiàn)象，在離心力的作用下，冰?？拷芡鈧?cè)速度明顯大于內(nèi)側(cè)，逐漸在外側(cè)堆積，這與文獻(xiàn)[21]的結(jié)果相符。

圖4 90°彎管截面上流動(dòng)速度分布

綜合以上分析，在90°彎管處，由于存在摩擦和二次流現(xiàn)象[22]，冰漿流體流動(dòng)阻力大，擾動(dòng)強(qiáng)烈，摩擦損失加大。而流速?gòu)娜肟谔幹饾u增大，在彎曲處流速最大且開(kāi)始出現(xiàn)分層現(xiàn)象，在流速較大區(qū)域，冰-水換熱效果緩慢，冰粒不易完全融化[23]。因此，溶液平均溫度隨流速增大而升高，在彎曲處溫度高于入口處，且冰的密度略小于水，由于受到浮升力作用和彎管處冰漿流速較快的影響，冰粒在流體的帶動(dòng)下，逐漸沿著管壁一側(cè)分散開(kāi)來(lái)。經(jīng)過(guò)彎曲處，平均流速減小，溫度稍降低。溶液溫度分布從大到小依次為彎曲處、出口處和入口處，但總體溫度比較均勻。

3.2 T 型管內(nèi)冰漿換熱特性結(jié)果分析

含冰率為10%時(shí)，T 型管縱截面處冰粒焓值分布，如圖5所示。

圖5 T 型管縱截面冰粒焓值分布

由圖5 可知，入口處溶液的焓值較高，隨著管道流動(dòng)呈梯形下降趨勢(shì)，在分流處?kù)手禍p小。這是由于冰漿在管內(nèi)流動(dòng)中，摩擦損失和紊流效應(yīng)使能量損失。在交叉處，擾動(dòng)強(qiáng)烈，液態(tài)水與冰粒之間交換動(dòng)量。分流之后固液兩相流動(dòng)趨于穩(wěn)定，在出口處壓力再次降低。

3.3 直管內(nèi)冰漿流動(dòng)冰堵分析

為了解決冰漿堵塞造成換熱效果差這一問(wèn)題，對(duì)直管中冰漿的分布進(jìn)行數(shù)值模擬。管道直徑為0.02 m；入口邊界條件為速度型，出口為出流型。重力方向?yàn)閅軸負(fù)方向。當(dāng)含冰率為20%，取管道中間橫截面位置進(jìn)行分析，得到冰粒的體積分?jǐn)?shù)分布結(jié)果如圖6所示。可以分為低體積分?jǐn)?shù)區(qū)、高體積分?jǐn)?shù)區(qū)和過(guò)渡區(qū)。圖7所示為文獻(xiàn)[24]中的3 種冰漿的流動(dòng)流態(tài)。

圖6 不同流速下冰粒在直管橫截面上體積分?jǐn)?shù)分布

圖7 文獻(xiàn)[24]中3 種冰漿流動(dòng)流態(tài)

由圖6(a)和圖6(b)可知，當(dāng)流速為0.2 m/s 和0.3 m/s 時(shí)，冰粒體積分布分層明顯，上層冰粒體積分?jǐn)?shù)最高約為0.4，底部冰粒體積分?jǐn)?shù)最低，幾乎為0，中間為過(guò)渡區(qū)，這與圖7(a)對(duì)應(yīng)，由于流速較小，在高體積分?jǐn)?shù)區(qū)黏滯力較大，使冰層沉積，管道堵塞。過(guò)渡區(qū)流速較高、體積分?jǐn)?shù)區(qū)大，小于臨界沉積速度，仍有可能造成冰堵。當(dāng)流速為1 m/s時(shí)，由圖6(c)可知，冰粒體積分?jǐn)?shù)分布整體較為均勻，上層為移動(dòng)床流動(dòng)，下層為非均質(zhì)流動(dòng)，即懸浮床流動(dòng)，與圖7(b)對(duì)應(yīng)。當(dāng)流速大于1.5 m/s 后，如圖6(d)和圖6(e)所示，此時(shí)，流速大于臨界沉積速度，整個(gè)管道內(nèi)呈現(xiàn)為懸浮床模型，冰粒速度相當(dāng)，體積分?jǐn)?shù)分布均勻，與圖7(c)相對(duì)應(yīng)，發(fā)生冰堵情況較少。

通過(guò)對(duì)不同流速下的冰粒體積分?jǐn)?shù)模擬結(jié)果的比較，選擇1.5 m/s 作為冰漿由移動(dòng)床流動(dòng)變?yōu)閼腋〈擦鲃?dòng)的第二臨界流速，冰粒分布趨于最均勻。

不同含冰率條件下，冰漿體積分?jǐn)?shù)的最值與流速的關(guān)系如圖8所示。由圖8 可知，隨著含冰率的增加，A 點(diǎn)的第一臨界速度逐漸減小，在1～2 m/s之間，C 點(diǎn)的第二臨界速度減小，在10 m/s 左右。其原因是含冰率增加對(duì)湍流效應(yīng)的影響大于對(duì)黏性摩擦的影響，因此，含冰率高的冰漿更有可能在較低流速下轉(zhuǎn)化為均勻流，減少冰堵發(fā)生。根據(jù)固體冰的分布輪廓，冰漿的總體積為：

圖8 不同含冰率條件下冰漿體積分?jǐn)?shù)最大值（最小值）隨流速的變化

4 結(jié)論

本文利用Fluent 軟件對(duì)冰粒在不同管道（90°彎管和T 型管）內(nèi)的換熱特性和直管內(nèi)分布情況進(jìn)行理論分析，得出如下結(jié)論：

1）當(dāng)管徑為0.02 m、管總長(zhǎng)為1 m 時(shí)，在90°彎管處，冰粒出現(xiàn)速度分層現(xiàn)象，管外側(cè)速度大于內(nèi)側(cè)；冰漿溫度在彎曲處最大，在出口處，溫度略有降低，但整體溫度分布均勻；T 型管內(nèi)溶液總焓值持續(xù)降低，經(jīng)過(guò)分流處后，焓值降低更加劇烈，能量損失明顯；

2）對(duì)冰堵風(fēng)險(xiǎn)的預(yù)測(cè)結(jié)果表明，當(dāng)直管管徑為0.02 m、管長(zhǎng)為1 m、流速大于1.5 m/s 時(shí)，流態(tài)為懸浮床流，冰粒體積分?jǐn)?shù)分布均勻，冰堵風(fēng)險(xiǎn)低；

3）當(dāng)直管管徑為0.02 m、管長(zhǎng)為1 m 時(shí)，得出冰漿最大值體積分?jǐn)?shù)、最小值體積分?jǐn)?shù)與流速之間的變化曲線，第一臨界速度和第二臨界速度分別為1～2 m/s 和10 m/s 左右，兩者均隨含冰率的增加而減小，含冰率高的冰漿更有可能在較低流速下轉(zhuǎn)化為均勻流，減少冰堵發(fā)生。