海水-冰晶在極地船換熱器中的相變模擬研究

徐 立，鄒祥巖，黃振飛，劉志平，孫 強

(1.武漢理工大學 高性能船舶技術教育部重點實驗室，湖北 武漢 430063；2.武漢理工大學 能源與動力工程學院，可靠性工程研究所，湖北 武漢 430063；3.武漢理工大學 物流工程學院，湖北 武漢 430063)

0 引 言

北極航道是溝通東亞、北美和歐洲的最方便的路線。與通過蘇伊士運河或巴拿馬運河的傳統商業路線相比，北極航道具有明顯的優勢，因為北極航道具有縮短航線，減少二氧化碳排放和降低運輸成本等優點。近年來，由于全球氣溫的上升，北極地區海冰面積陸續減少。以往的衛星觀測數據表明，過去幾十年北極海冰每10年下降約3%，其中夏季海冰融化最為顯著[1]。每年夏天，北極航道都有無冰期，這為北極商船提供了機會，與傳統的航道不同，商船在北極航道受到北極冰區的影響。當船舶在冰區航行時，船舶海底吸入口過濾器將被大塊冰覆蓋，并且冰晶顆粒將與海水一起進入海水系統，海水和冰晶粒子將形成固-液兩相流。兩相流的特性比單相流復雜得多，當海水和冰晶兩相流進入極地船舶的熱交換器時，冰晶會熔化，這對海水的傳熱特性有一定的影響。目前，大多數對極地船的研究主要集中在海冰的融化，航行的可行性分析，法律問題和對極地船舶運行的建議。對海水和冰晶兩相流的流動和相變特性研究的很少。

由于計算流體動力學（CFD）方法的經濟性和有效性以及實驗測量技術的局限性，CFD方法用于研究冰晶粒子在極地船殼管式換熱器的分布和相變特性。歐拉-歐拉模型是描述固液兩相流的有效模型[2]。Zhang成功用相間傳熱傳質模型研究了水平管道中氮漿熔融特性[3]。本文采用歐拉-歐拉模型與界面傳熱傳質模型研究冰晶的融化。隨著北極航道的開通，對極地船換熱器中冰晶的分布和熔化特性的研究具有重要意義，這也為極地船海水冷卻系統的穩定性提供了保障。

1 幾何模型

船舶經常使用的管殼式換熱器如圖1所示。本文研究了殼管式換熱器單個海水管中冰晶的熔化和分布，對整個熱交換器的研究起到了指導作用。常用的海水管是直管和U管。根據GB151-1999行業標準[4]，建立了直管和U管的三維幾何模型。管的直徑為24 mm，其長度為1 000 mm。

2 網格模型

在本研究中，使用商業計算流體動力學軟件Ansys ICEM15.0來研究冰晶-海水兩相流的流動和相變特性。首先，創建直管和U形管的三維模型。然后在模型中進行網格生成，直管和U形管的計算域分為六面體結構的網格單元。將直管分成114 915 6個網格單元，將U形管分成210 296 5個網格單元，模型進行網格劃分后網格質量大于0.65，符合仿真的網格質量要求。如圖2和圖3所示，為U型管和直管的網格模型。

圖2 U形管的網格模型Fig.2 Mesh model of U-shaped tube

圖3 直管的網格模型Fig.3 Grid model of straight pipe

3 數學模型

本研究中使用的CFD模型利用粒狀流動力學理論來描述粒子相互作用。為了簡化數學模型，冰晶顆粒被假定為球形，非彈性和光滑的，體積分數為15%，可以認為是牛頓流體[5]。在本研究中，用歐拉-歐拉模型研究管道中冰晶的分布，相間傳熱傳質模型模擬冰晶的相變。歐拉-歐拉雙流體模型中固相和液相的守恒方程包括連續性方程、動量守衡方程等。

3.1 連續性方程

式中：i表示固相或液相，當i=l時為液相，i=s時為固相，下標p表示i的相對相。在式（1）中，ai為每個相的體積分數；為每個相的速度；為每個相的密度；為相間傳質系數。每相的體積分數的關系可以表示為

3.2 動量守恒方程

固液兩相流，每個相都有自己的動量守恒方程。固相和液相的動量守恒方程略有不同。因此，本文分別描述了它們的動量守恒方程。

3.2.1 海水動量守恒方程

3.2.2 冰晶的動量守恒方程

3.3 能量守恒方程

3.4 相間傳熱傳質

在本文中，相間傳熱傳質模型用于模擬由冰晶融化引起的相間熱和質量傳遞。在能量方程中，表示由冰晶相變引起的動量交換，表示冰晶相變引起的熱量交換。通過計算特定界面面積的單個顆粒與流體之間的傳熱系數到相間傳熱系數，該系數可由下式表示：

冰晶融化向海水傳遞的質量

其中ΔH表示冰晶顆粒的潛熱。

3.5 相間作用力

相間作用力包括拖曳力和湍流擴散力。拖曳力可以表示為

CD為拖曳力系數，可表示為

3.6 顆粒流動力學理論

因為離散相被認為是連續相，離散相具有與連續相相似的性質，如體積粘度，剪切粘度, 顆粒壓力等，這些性質也被稱為“假流體性質”。顆粒流動力學理論用于描述冰晶顆粒的假流體性質，該理論基于氣體的動力學理論提出，并且顆粒被認為是致密的空氣分子。顆粒溫度被提出來描述顆粒的波動能量，其被定義為

從動力學理論導出的方程表示為

體積粘度ξs用于表征顆粒抵抗變形的能力，其可以計算為

顆粒剪切粘度由2部分組成。

3.7 湍流模型

4 數值參數設置

海水和冰晶的熱物性與海水的鹽度有很大的關系。海水的平均鹽度為35‰。由于海冰的融化，北極海水的鹽度必須降低。在本研究中，海水的鹽度被分配為15‰。海水和冰晶的熱物性如表2所示，鹽度為15‰[10,11]。

在本研究中，采用速度入口邊界條件，速度設置為1 m/s，1.5 m/s和2 m/s。冰晶的體積分數為15%。根據王治[12]的計算，在取海水管與被冷卻水之間的換熱系數為 759 W·m–2·K–1和被冷卻水的溫度為 50℃。使用自由流出口邊界條件。管道入口處冰粒和海水的溫度分別為272.34 K和273.15 K。為了求解控制方程，使用有限體積方法來離散這些方程，并且使用二階方案。相耦合SIMPLE算法用于求解這些離散方程。時間步長設為0.005 s，當殘差小于1.0×10–3時，數值結果被認為是收斂的。 在本研究中，實現了網格獨立性。

表2 海水和冰晶的熱物性Tab.2 Thermophysical properties of seawater and ice crystals

5 仿真分析

5.1 直管和U管中冰晶顆粒的分布

從圖4～圖7可以看出直管中冰晶顆粒的體積分數特性，冰晶顆粒主要集中在主流區。由于冰晶顆粒的密度較小，在直管的上部有更多的冰顆粒。從圖7中可以看出隨著入口速度的增加，湍流強度更強，分布情況受冰晶顆粒密度影響變小，冰晶顆粒的分布更加均勻。在邊界層區域中，由于冰晶顆粒和壁之間的碰撞和冰晶顆粒的熔化，導致冰晶顆粒減少，這種現象在入口速度增加時更為明顯。

圖4 當速度為1 m/s時，直管出口冰晶的體積分數Fig.4 Volume fraction of ice crystals at straight exit at 1 m/s

圖5 當速度為1.5 m/s時，直管出口處冰晶的體積分數Fig.5 Volume fraction of ice crystals at straight exit at 1.5 m/s

圖6 當速度為2 m/s時，直管出口處冰晶的體積分數Fig.6 Volume fraction of ice crystals at straight exit at 2 m/s

圖7 當速度為1 m/s，1.5 m/s和2 m/s時，直管中冰晶的體積分數Fig.7 Volume fraction of ice crystals in straight tubes at speeds of 1 m/s, 1.5 m/s and 2 m/s

圖8～圖13描述了U型管中冰晶顆粒的體積分數。從圖9和圖10比較中可看出，由于湍流的影響，冰顆粒的分布更加均勻。從圖11可以看出，由于離心力的作用，更多的冰晶顆粒集中在U形管的拐彎處外側，隨著海水速度的增加，離心力的影響更加明顯，更多的冰粒子移動到外側。從圖8、圖9和圖10看出，由于顆粒和管壁之間的碰撞，冰晶顆粒通過U型管的拐角處時，顆粒逐漸移動到流場的主流場。

5.2 冰晶顆粒的熔化特性

圖8 當速度為1 m/s，1.5 m/s和2 m/s時，U管中冰晶的體積分數Fig.8 Volume fraction of ice crystals in U tube at speeds of 1 m/s,1.5 m/s and 2 m/s

圖9 當速度為1 m/s和1.5 m/s時，U管出口處的冰晶體積分數Fig.9 Ice volume fraction at the exit of the U tube when the speed is 1 m/s and 1.5 m/s

圖10 當速度為2 m/s時，U管出口處的冰晶體積分數Fig.10 The volume fraction of ice crystals at the exit of the U tube when the speed is 2 m/s

圖11 當流速為1 m/s時，U管拐角處冰晶體積分數Fig.11 Fractional volume of ice crystals at the corners of the U tube when the flow rate is 1 m/s

圖12 當流速為1.5 m/s時，U管拐角處冰晶體積分數Fig.12 Frozen volume fraction at the corners of the U tube when the flow rate is 1.5 m/s

圖13 當流速為2 m/s時，U管拐角處冰晶體積分數Fig.13 Ice flow fraction at the corner of the U tube when the flow rate is 2 m/s

圖14 不同的入口速度，直管中離入口處每一固定截面冰晶的平均體積分數變化Fig.14 Volume fraction of fixed section ice crystals at the inlet of the sraight tube at different speeds

在本研究中，冰晶顆粒的熔融特性可用冰晶從入口處到出口處的體積分數來表示。圖14至16描述了直管和U管中的冰晶顆粒的融化特性。在圖14中，由于冰晶顆粒和海水之間的溫度差的增加，冰晶的體積分數沿軸向方向逐漸減小，并且在直管的入口區域中的海水的溫度比在中心區域中增加地更快，因此，在入口區域時，冰晶的體積分數變化的最快，曲線下降的最快，并且隨著入口速度的不同，冰晶的體積分數變化情況不同，速度越大，湍流波動增加，導致更多的冰粒隨著管道入口速度的增加而融化。而在圖15和圖16中，分別描述了U型管下部和上部的冰晶顆粒的體積分數變化情況，并表明U型管扁平部分的冰晶顆粒融化特性與直管類似。但是，在入口速度相同時，直管和U管中的冰晶的融化特性不同，計算結果表明，由于更多的熱量傳入U形管中，U形管中有更多的冰粒融化，并隨著入口速度的增加，冰晶的融化速度加快。

圖15 不同的入口速度，U管中離入口處每一固定截面冰晶的平均體積分數變化Fig.15 Volume fraction of fixed section ice crystals at the inlet of the U tube at different speeds

圖16 不同的入口速度，U管中離出口處每一固定截面冰晶的平均體積分數變化Fig.16 Volume fraction of fixed section ice crystals at the exit of the U tube at different speeds

6 結 語

本次研究調查了極性船舶熱交換器海水管中冰粒的分布和相變特性，可以得出冰晶-海水兩相流的流動模式是懸浮的，隨著入口速度的增加，湍流強度增強，冰晶顆粒的分布更加均勻。冰晶顆粒在U型管的拐角處受到離心力的作用，會從U型管內壁流到U型管角部的外壁。通過比較出口處冰晶的體積分數，發現在進口速度相同的情況下，U型管中有更多的冰粒融化，并且融化的冰晶顆粒量隨入口速度的增加而增加。

[1]沈春, 施偉來.基于Nimbus-7和DMSP衛星資料的北極海冰變化特征分析[J].海洋預報, 2015, 04: 31–43.

[2]GUPTA A, ROY S.Eulerian-eulerian simulation of bubbly flow in a rectangular bubble column: experiment validation with radioactive particle tracking[J].Chem.Eng.J, 2013(225):818–836.

[3]ZHANG P, JIANG Y Y.Forced convective heat transfer of slush nitrogen in a horizontal pipe[J].Int.J.Heat Mass Transfer,2014(71): 158–171.

[4]錢頌文.換熱器設計手冊[M].北京: 化學工業出版社, 2002:32–34.

[5]AYEL V, LOTTIN O, PEERHOSSAINI H.Rheology flow behavior and heat transfer of ice slurries: a view of the state of the art[J].International Journal of Refrigeration, 2003, 26:95–107.

[6]GUNN D J1.Transfer of heat or mass to particles in fixed and fluidized beds[J].Int.J.Heat Mass Transfer, 1978(21):467–476.

[7]BURNS A D B, FRANK TH, HAMILL I, et al.The favre averaged drag model for turbulent dispersion in eulerian multiphase flows[C]//Fifth International Conference on Multiphase Flow, ICMF-2004, Yokohama, Japan.2004.

[8]EKAMBARA K, SANDERS R S, NANDAKUMAR K, et al.Hydrodynamic simulation of horizontal slurry pipeline flow using ANSYS-CFX[J].Ind.Eng.Chem.Res, 2009(48):8159–8171.

[9]GIDASPOW D, BEZBURUAH R, DING J.Hydrodynamics of circulating fluidized beds, kinetic theory approach[C]//In Fluidization VII, Proceedings of the 7th Engineering Foundation Conference on Fluidization.Australia, Brisbane,1992: 75–82.

[10]朱傳芳.海水制冰機制冰傳熱過程特性研究[D].廣州：華南理工大學, 2010: 21–24.

[11]MOSTAFA H S, JOHN H L, SYED M Z.On thermal performance of seawater cooling towers[J].Journal of Engineering for Gas Turbines and Power, 2011.133(4): 1–7.

[12]王冶.船用換熱器流動與傳熱的三維流場數值模擬研究[D].上海: 上海交通大學, 2013: 23–14.