趙新穎，黃溫赟，黃文超，呂續艦

(1 中國水產科學研究院漁業機械儀器研究所，農業農村部遠洋漁船與裝備重點實驗室，上海 200092;2 海洋國家實驗室深藍漁業工程裝備技術聯合實驗室，山東 青島 266237;3 南京理工大學能源與動力工程學院，江蘇 南京 210094)

冰漿是一種由平均直徑約0.10 mm的冰粒與水溶液組成的混合物[1]，具有良好的流動性和較大的蓄冷密度，是目前廣泛應用的蓄冷二次循環制冷介質之一[2-3]。為降低建造成本，國內漁船大部分未安裝壓縮式制冷裝置，一般攜帶碎冰給漁獲物保鮮，效果不佳，影響漁獲物品質。利用冰漿的流動傳熱對漁獲物進行降溫冷藏，確保漁獲物品質安全，是一種切實可行的方法。研究冰漿的流動和傳熱特性，掌握船載冰漿循環理論，具有重要現實意義。

國內學者在冰漿領域開展了大量頗具成效的研究工作。江煥寶[4]通過對某極地運輸船海水系統中作為滑油冷卻器的管殼式換熱器的數值模擬，實現了換熱器的結構參數優化。劉圣春等[5]總結了換熱器的流動換熱特性研究現狀，認為確定兩相顆粒流動的流變模型是解決目前換熱問題的前提。文力等[6]對冰漿制取裝置的過冷器中流體流動進行了模擬分析。楊帆等[7]利用Bingham模型分析了冰漿在層流和湍流狀態下的摩擦因子,得到了冰漿在水平直管內流動的壓降曲線。熊庭等[8]對水平管道中泥漿輸送的固液兩相流進行了CFD數值模擬，并與Durand 模型的計算值和文獻試驗數據進行了對比分析。徐愛祥等[9]通過建立 CFD-PBM 耦合模型，研究含冰率與流速對水平直管內冰漿的冰晶體積分數與粒徑的分布規律。鄢新[10]采用Fluent對冰水固流轉化的傳熱過程進行了仿真模擬，根據模擬結果初步給出了冰水界面熱傳導問題的解決辦法。葉健[11]通過試驗研究了質分數小于19.2%的TBAB水合物漿體在內徑為2.0 mm和4.5 mm圓直管中的流動特性，通過測量水合物漿體的流速和壓降得到了相應的流動曲線。

此外，何國庚等[12]早在2005年就注意到了冰漿的應用前景，并對冰漿流體的應用研究提出了一些有價值的建議。王繼紅等[13-14]采用計算流體力學方法研究了水平管道內冰漿流體的無相變流動過程，給出了管道內冰漿流體阻力特性；鄧義斌等[15]基于商用CFD軟件Fluent對水平管內冰水兩相流動進行了計算研究，發現管內兩相流動阻力特性受顆粒碰撞彈型恢復系數、水流流速、顆粒直徑、顆粒質量流量等特性的影響較為顯著。張鵬等[16]同樣考慮水平管內的冰漿流動情況，不同之處在于考慮了冰漿在三相點處溫度下水平管內的流動情況；白銀等[17]通過試驗對不同時間、不同加熱量、不同流速下傳熱器進出口的冰漿的溫度、含冰率進行了測試；徐立等[18-19]研究了冰漿在極地船上的應用，包括極地船換熱器中海水-冰晶兩相流的流動及傳熱特性，以及冰晶顆粒在極地船殼管式換熱器海水管內的分布和融化特性[20-21]。

相較于碎冰等傳統方式，采用冰漿對漁獲物進行保鮮能夠提供更為均勻的冷量分布，但相關研究工作鮮有報道。本研究建立了冰漿流動的計算模型，首先采用文獻數據對計算模型進行了驗證，隨后研究了入口冰顆粒體積分數對漁獲物上表面與側面的熱通量特性的影響，最后研究了流動速度對漁獲物表面傳熱量的影響，力圖為冰漿在漁獲物保鮮中的應用提供具有一定實用價值的參考。

1 數學模型

1.1 模型概述

采用歐拉-歐拉雙流體模型對冰漿流動傳熱進行仿真計算，其基本思想是將液體相和顆粒相等效為可以互相貫穿的連續性介質，針對液體相和顆粒擬流體相分別建立N-S方程，同時液體相和顆粒相之間的耦合是通過相間作用力來實現的。顆粒間的相互作用以碰撞為主，壁面處通過設置鏡面反射系數和彈性恢復系數來表現顆粒與壁面的相互作用。其連續性方程為[16]：

(1)

式中：α、ρ和v分別為兩相的體積分數、密度(kg/m3)及速度矢量(m/s)；i表示l(液相)或者s(固相)，下文中的式(3)和(5)～(7)中下標q表示與i相對的相(i=l，q=s；i=s，q=l)。

1.2 動量守恒方程

動量守恒方程可以表示為[16]：

(2)

式中：τ為應力張量(N/m2)；P為流場壓力分布；g是重力加速度(m/s2)；FD,i和FL,i分別表示冰顆粒在第i相中的拖曳力和升力(N)；為哈密頓算子。對于液相有：

(3)

而對于固相，根據擬流體假設，其應力張量τS由顆粒碰撞產生的隨機粒子運動引起：

(4)

式中：μS、ξS和PS分別表示剪切黏度(N·s/m2)、體積黏度(N·s/m2)以及固相壓強(N·s/m2)；為單位向量。

1.3 相間作用力

在歐拉雙流體模型中，兩相間的耦合關系是通過相間作用力建立起來的，本次計算所考慮的相間作用力包括拖曳力和升力。拖曳力是顆粒在液相中的阻力，是固液兩相間最基本的相互作用形式[16]，其表達式為:

FD,i=Ksl(vq-vi)

(5)

式中：Ksl表示顆粒的拖曳系數，可采用Gidaspow模型進行計算：

(6)

顆粒在液相中運動時還受到由液相速度梯度引起的升力，其表達式為：

FL,i=CLαSρl(|vl-vs|)×(×vl)

(7)

式中：CD和CL分別為顆粒阻力系數和升力系數。

結合上述模型，考慮顆粒動力學理論、標準k-ε湍流模型和能量方程模型，即可實現冰漿流動和傳熱特性分析[16]。

2 水平管內冰漿流動模型驗證

2.1 初始條件及邊界條件

初始時水平管內冰顆粒均勻分布，速度垂直于入口截面，且固相與液相速度相同。管道出口采用壓力出口條件，出口為一個大氣壓。壁面采用絕熱條件，液相無滑移，冰顆粒在壁面處采用Johnson-Jackson模型處理，顆粒與壁面碰撞恢復系數為0.9，鏡面系數為0.015。

2.2 網格劃分及計算方法

首先在ANSYS-ICEM中建立了水平管道三維模型，管道直徑D為23 mm，為了保證流動充分發展，管道長度L設置為2 000 mm。管道網格采取了六面體結構網格，網格數為72萬個，采用非穩態求解模型，求解方法為壓力與速度耦合的SIMPLE求解器，控制方程采用QUICK與二階迎風離散格式。開始計算時時間步長設置為較小的0.000 5，以保證收斂性；每個時間步長內控制方程的殘差小于10-4時認為計算收斂。

2.3 計算結果及驗證

在Fluent軟件中對水平管內冰漿兩相流動進行數值求解，計算了在冰顆粒體積分數為0.1、速度為1.0 m/s時的冰顆粒體積分數分布情況，并與Zhang等[1]計算結果進行驗證，以保證利用這種方法模擬冰漿與漁獲物流動傳熱的可行性與可靠性，如圖1所示，吻合情況較好。由于冰顆粒的密度略小于海水，因此在流動過程中冰顆粒會不斷向上懸浮，使得出口截面上的冰顆粒體積分數自上而下遞減，圓管頂部最高為0.103 5，底部最低為0.09。

圖1 水平管出口截面冰顆粒體積分數驗證(U=1.0 m/s)

3 方管內冰漿與漁獲物的流動傳熱計算

3.1 物理模型及計算條件

在冰漿的管道輸送中，外界的傳熱容易使冰漿在傳輸過程中發生相變。根據漁獲物尺寸，在ANSYS ICEM中建立方管模型并劃分網格。如圖2所示，方管截面為邊長4 m的正方形，長度20 m；管內漁獲物的長、寬、高分別為3.2 m、2.36 m、2.1 m，網格數為48萬個。

圖2 方管內冰漿和漁獲物流動傳熱模型

由圖2可知，漁獲物在方管內沿x軸向的位置為x=10～13.2 m，漁獲物頂部在截面y=0.1上，兩側面分別位于z=-1.18和z=1.18截面上。

采用歐拉-歐拉雙流體模型對模型進行計算，在初始時刻方管內冰顆粒均勻分布，固相與液相溫度均為273.15 K；冰漿速度垂直于入口截面，且固相與液相速度相同；管道出口設置為壓力出口，條件為一個大氣壓；管內漁獲物邊界條件設置為wall，給定壁溫298.15 K，且壁面處考慮冰顆粒的碰撞與反射。采用耦合SIMPLE求解器進行穩態計算，控制方程采用二階迎風方法進行離散。

3.2 計算結果

3.2.1 不同冰顆粒體積分數對傳熱效果的影響

熱通量(Heat flux)又稱為熱流，表示單位時間通過單位面積的熱能，可以用傅里葉定律來充分描述。與溫度相關的熱通量可以描述為：

(8)

式中：k為表面熱導率；τ為溫度分布；x為沿表面的某一坐標，負號表示熱通量從高溫區傳向低溫區。在對漁獲物與流體傳熱進行計算時，利用了冰顆粒與漁獲物表面之間的熱通量來表征冰漿對漁獲物的冷卻情況，在CFD-Post 19.0中讀取漁獲物表面傳熱量，全管如圖3所示。

圖3 熱通量分布情況(U=0.5 m/s，αin=0.2)

冰顆粒的分布是影響傳熱的主要原因之一，圖4展示了流速為0.5 m/s，冰顆粒入口體積分數αin=0.2時漁獲物附近截面上冰顆粒分布情況。

圖4 漁獲物附近冰顆粒分布情況

從圖4中可以看出，在冰漿流動充分發展以后，由于冰的密度較小，在重力作用下冰顆粒大量聚集在方管頂部，在方管出口截面(x=20)處冰顆粒體積分數最高為0.397；管道中間部分則較為均勻，體積分數在0.2～0.22；而在漁獲物附近，方管底部幾乎沒有冰顆粒分布；在冰漿繞流漁獲物時，由于渦的產生以及壁面處速度梯度較大，湍流脈動較強，導致漁獲物附近冰顆粒質量濃度較低，且分布不均勻，繞流結束后冰顆粒分布有重新趨于穩定的趨勢。

冰漿質量濃度主要由冰顆粒體積分數Vf體現，傳熱效果與冰顆粒的分布情況密切相關。當一定流速的冰漿在方管中沖擊漁獲物迎流面后會發生繞流，冰顆粒與漁獲物表面之間發生碰撞反射，與顆粒自身的彈性碰撞以及冰水兩相相互作用，必然導致冰顆粒分布不均勻，使得漁獲物各表面傳熱效果不同，當冰漿繞流漁獲物的頂部與側面時，這種情況尤為明顯。故分別計算了αin=0.1、0.2、0.3時，漁獲物上表面(x=10)與側面(z=1.18)冰顆粒體積分數Vf與熱通量[q(W/m2)]的分布情況，并讀取上表面和側面中線上的數據，其變化趨勢曲線如圖5所示。

圖5 不同αin下漁獲物上表面與側面熱通量分布情況

圖6中漁獲物上表面冰顆粒質量濃度沿z軸的變化較為明顯，其總體呈中心低、兩側高的規律分布，且隨著入口冰顆粒體積分數αin的提高，表面熱通量呈現出整體增大的趨勢。圖7展示出漁獲物表面附近冰漿的流線分布情況。可以看出：

圖6 不同αin情況下漁獲物上表面和側面熱通量和冰顆粒體積分數變化趨勢

圖7 漁獲物上表面附近流線分布

(1)上表面中心處熱通量較小，原因可能為當冰漿繞流90°棱角的漁獲物時，上表面位于方管中心位置，流速較高，顆粒與壁面碰撞所受的排斥力較大，使得湍流脈動較為強烈，冰漿繞流后在漁獲物表面上方形成渦流，使冰顆粒有遠離壁面的趨勢，進而影響了顆粒與壁面間的傳熱，而在αin=0.3的工況下，發現冰顆粒體積分數與表面熱通量的波動明顯小于另外兩種工況；

(2)在上表面兩側冰顆粒的質量濃度和熱通量都相對較高，αin為0.1、0.2和0.3時熱通量最高分別達到9 221.9 W/m、21 396.5 W/m和21 375.3 W/m。因為在冰漿繞流過程中，由于冰顆粒密度較低，且繞過棱角時速度高，壓力較低，造成漁獲物兩側的流體向上表面翻轉的現象，因此上表面兩側冰顆粒質量濃度較高，熱通量明顯大于中心位置；

(3)這種向上翻轉的情況同樣使得漁獲物側面的熱通量分布不均，側面與上表面接觸的位置熱通量明顯高于其他位置，在αin為0.2時最高達到21 612.3W/m，超過αin為0.1時的2倍；

(4)在重力作用下，冰顆粒隨著流動過程逐漸向管道頂部聚集，使得漁獲物側面的冰顆粒體積分數與熱通量自上而下呈遞減趨勢。

3.2.2 不同流速對傳熱效果的影響

對于固液兩相流動而言，在近壁面區域，冰顆粒所受到的升力和顆粒與壁面之間碰撞所產生的排斥力對其分布情況影響較大；而升力與液相速度梯度成正比，顆粒與壁面的碰撞作用主要受到入口體積分數αin的影響，因此在相同入口體積分數下，排斥力隨著入口速度的增大而增大，故管道上壁面聚集的冰顆粒隨著入口速度的增大呈減小的趨勢，冰顆粒會逐漸趨向于管道中間位置，如圖8所示。對αin=0.2，入口速度分別為0.5 m/s、1.5 m/s和2.5 m/s的情況下漁獲物表面冰顆粒分布情況以及傳熱情況計算發現：

圖8 不同流速下漁獲物表面冰顆粒體積分數與熱通量變化曲線

(1)隨著流動速度的增大，漁獲物表面冰顆粒體積分數和熱通量均有所提高，其中U=1.5 m/s情況下，在靠近側棱約0.2 m的距離內熱通量略高于2.5 m/s的工況，且側棱處冰顆粒體積分數基本相同，其原因可能為在速度較高時，側棱處發生繞流時冰顆粒慣性較大，導致和漁獲物表面距離較大，造成熱通量反轉的情況；

(2)在流速為1.5 m/s和2.5 m/s的工況下，漁獲物側面熱通量在靠近側棱處的冰顆粒質量濃度和熱通量相差無幾；相較于上表面，漁獲物側面冰顆粒體積分數和熱通量的變化較大，在2.5 m/s工況下漁獲物側面最高熱通量約是0.5 m/s工況下的4.5倍，分別為104 822 W/m和23 132 W/m，原因可能是在流速較低的情況下冰顆粒易于向方管頂部聚集，導致繞流漁獲物的下層冰漿冰顆粒體積分數較小，因此傳熱效果相對較差。

4 結論

冰顆粒的分布情況是影響其與漁獲物表面傳熱的主要因素，隨著冰顆粒入口體積分數的提高，漁獲物上表面與側面的熱通量整體有所增大；隨著入口冰顆粒體積分數的增大，其在漁獲物表面的分布趨于均勻。通過對3種流動速度的工況進行對比可以發現，適當提高流動速度可以增大漁獲物表面冰顆粒的體積分數，從而增加表面熱量交換。當流動速度增大時，漁獲物表面的冰顆粒數量一開始明顯增大，當速度增大到一定水平時，漁獲物表面上冰顆粒體積分數十分接近，此時漁獲物上表面傳熱量增幅很小。

□