海水-冰晶在极地船换热器中的相变模拟研究

徐 立，邹祥岩，黄振飞，刘志平，孙 强

(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.