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

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

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

0 引 言

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

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

1 幾何模型

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

2 網(wǎng)格模型

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

圖2 U形管的網(wǎng)格模型Fig.2 Mesh model of U-shaped tube

圖3 直管的網(wǎng)格模型Fig.3 Grid model of straight pipe

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

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

3.1 連續(xù)性方程

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

3.2 動量守恒方程

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

3.2.1 海水動量守恒方程

3.2.2 冰晶的動量守恒方程

3.3 能量守恒方程

3.4 相間傳熱傳質(zhì)

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

冰晶融化向海水傳遞的質(zhì)量

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

3.5 相間作用力

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

CD為拖曳力系數(shù)，可表示為

3.6 顆粒流動力學(xué)理論

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

從動力學(xué)理論導(dǎo)出的方程表示為

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

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

3.7 湍流模型

4 數(shù)值參數(shù)設(shè)置

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

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

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

5 仿真分析

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

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

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

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

圖6 當速度為2 m/s時，直管出口處冰晶的體積分數(shù)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時，直管中冰晶的體積分數(shù)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型管中冰晶顆粒的體積分數(shù)。從圖9和圖10比較中可看出，由于湍流的影響，冰顆粒的分布更加均勻。從圖11可以看出，由于離心力的作用，更多的冰晶顆粒集中在U形管的拐彎處外側(cè)，隨著海水速度的增加，離心力的影響更加明顯，更多的冰粒子移動到外側(cè)。從圖8、圖9和圖10看出，由于顆粒和管壁之間的碰撞，冰晶顆粒通過U型管的拐角處時，顆粒逐漸移動到流場的主流場。

5.2 冰晶顆粒的熔化特性

圖8 當速度為1 m/s，1.5 m/s和2 m/s時，U管中冰晶的體積分數(shù)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管出口處的冰晶體積分數(shù)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管出口處的冰晶體積分數(shù)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管拐角處冰晶體積分數(shù)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管拐角處冰晶體積分數(shù)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管拐角處冰晶體積分數(shù)Fig.13 Ice flow fraction at the corner of the U tube when the flow rate is 2 m/s

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

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

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

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

6 結(jié) 語

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

[1]沈春, 施偉來.基于Nimbus-7和DMSP衛(wèi)星資料的北極海冰變化特征分析[J].海洋預(yù)報, 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]錢頌文.換熱器設(shè)計手冊[M].北京: 化學(xué)工業(yè)出版社, 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].廣州：華南理工大學(xué), 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]王冶.船用換熱器流動與傳熱的三維流場數(shù)值模擬研究[D].上海: 上海交通大學(xué), 2013: 23–14.