船舶在极地航行时，换热管在冷却系统运行中具有重要作用，必须使其能够高效稳定的运行。由于海水中存在冰晶介质，一定条件下产生堵塞现象^[1]，使换热器无法正常工作；船舶上机械设备的运行产生固有频率在0~100 Hz振动^[2]，使管道内海水-冰晶两相流的流动状态更为复杂。因此，有必要综合研究振动对直管中海水-冰晶两相流流动特性的影响。

近年来，随着振动强化传动技术的发展，Sun等^[3]通过CFD(计算流体动力学)中的离散相模型、冲蚀磨损模型和动网格技术，研究不同振动工况下海水-冰晶两相流对管道冲蚀磨损特性的影响。Xu^[4]建立了CFD模型和数群平衡模型(PBM)相耦合方法。研究水平管道中冰浆流动，表明冰粒大小对冰浆流动特性有明显影响。刘圣春等^[5]通过以颗粒相动力学为基础的欧拉-欧拉模型，研究不同管道模型内冰浆的压降特性，表明直管、弯管、T形管压力随流速增大而增加，且T形管的压降最大。王继红等^[6]通过CFD研究液固两相流在水平直管中的流动特性，颗粒的密度和湍流强度对液固两相流流动状态起到关键作用。Bordet等^[7]通过实验研究了入口平均轴向速度和冰体积分数对冰浆流动模式和压降的影响，在层流和湍流状态下，体积冰浓度高达18.4%。徐立等^[8]基于热焓多孔介质模型建立数学模型，研究流速对管道换热产生的作用。Man等^[9]对管换热器中插入纽带实验研究，在低雷诺数下提高热交换器的传热效率，但可能造成压降增大。但目前对于海水-冰晶两相流的流动特性研究较少，对添加振动的海水-冰晶两相流还未涉及，振动对于管道的摩擦损失也不容忽略。

基于船舶航行的实际工况中振动复杂及不规则性，实际中实验仿真不易，相关研究较少。本文使用Euler双流体模型和UDF自定义动网格技术，研究不同振动下水平直管中海水-冰晶两相流的流动特性，探讨振动对冰晶颗粒分布及流动阻力的影响。

1 数学模型

海水-冰晶两相流流动具有一定的复杂性，本研究通过假设冰晶为球形、光滑非弹性的颗粒简化数学模型，不考虑流动过程中的粘性耗散^[10]。冰晶颗粒在管道中的流动分布通过欧拉-欧拉双流体模型描述，将海水和冰晶看作可以相互贯穿的连续性介质，分别建立相应的Navier-Stokes方程。

1.1 连续性方程

质量守恒是流体力学中每一个流动介质微元体必须遵循的基本规律，表达式为：

$ \frac{\partial }{{\partial t}}\left( {{\alpha _i}{\rho _i}} \right) + \nabla \cdot \left( {{\alpha _i}{\rho _i}{\mathit{\boldsymbol{v}}_i}} \right) = 0 $

(1)

式中：α_i代表各相体积分数；下标i为l、s时分别代表海水和冰晶颗粒；ρ_i为各相密度；v_i为各相的速度；$ \nabla $ ·为散度。

局部固相与液相的体积分数相关性为：

$ {\alpha _s} + {\alpha _l} = 1 $

(2)

1.2 动量守恒方程

由于固液相的作用力不同，用相间动量传递表达动量守恒。海水动量守恒为：

$ \begin{array}{*{20}{c}} {\frac{\partial }{{\partial t}}\left( {{\alpha _L}{\rho _L}{\mathit{\boldsymbol{v}}_L}} \right) + \nabla \cdot \left( {{\alpha _L}{\rho _L}{\mathit{\boldsymbol{v}}_L}{\mathit{\boldsymbol{v}}_L}} \right) = }\\ { - {\alpha _L}\nabla \cdot p + \nabla \cdot {\tau _L} + {\alpha _L}{\rho _L}g + {M_{SL}}} \end{array} $

(3)

$ {{\tau _L} = {\alpha _L}{\mu _L}\left( {\nabla {\mathit{\boldsymbol{u}}_L} + \nabla \mathit{\boldsymbol{u}}_L^{\rm{T}}} \right) + {\alpha _L}\left( {{\lambda _L} - \frac{2}{3}{\mu _L}} \right)\left( {\nabla \cdot {\mathit{\boldsymbol{u}}_L}} \right)\mathit{\boldsymbol{I}}} $

(4)

冰晶颗粒动量守恒：

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _S}{\rho _S}{\mathit{\boldsymbol{v}}_S}} \right) + \nabla \cdot \left( {{\alpha _S}{\rho _S}{\mathit{\boldsymbol{v}}_S}{\mathit{\boldsymbol{v}}_S}} \right) = - {\alpha _S}\nabla \cdot p + \nabla \cdot {\tau _S} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\alpha _S}{\rho _S}g - {\alpha _S}\nabla \cdot {p_S} + {M_{LS}}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\tau _S} = - {P_S}I + {\alpha _S}{\mu _S}{\left( {\nabla \cdot {\mathit{\boldsymbol{u}}_S} + \nabla {\mathit{\boldsymbol{u}}_S}} \right)^{\rm{T}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\alpha _S}\left( {{\lambda _S} - \frac{2}{3}{\mu _S}} \right)\left( {\nabla \cdot {\mathit{\boldsymbol{u}}_S}} \right)\mathit{\boldsymbol{I}} \end{array} $

(5)

式中：g表示重力加速度；P表示静压；P_s表示冰晶颗粒相间的正应力；μ表示剪切黏度；τ表示剪切应力；λ表示粘滞系数；M表示海水与冰晶颗粒两相间的动量交换；I表示单位张量。

1.3 能量守恒方程

$ \frac{{\partial t}}{{\partial \theta }} + {\mathit{\boldsymbol{u}}_i}\frac{{\partial t}}{{\partial x}} + {\mathit{\boldsymbol{v}}_i}\frac{{\partial t}}{{\partial y}} + {\mathit{\boldsymbol{w}}_i}\frac{{\partial t}}{{\partial z}} = {a_i}\left( {\frac{{{\partial ^2}t}}{{\partial {x^2}}} + \frac{{{\partial ^2}t}}{{\partial {y^3}}} + \frac{{{\partial ^2}t}}{{\partial {z^2}}}} \right) $

(6)

式中：u_i、v_i、w_i分别表示各相在x、y、z方向上的速度，m/s。a_i为各项的扩散率，m²/s。

1.4 湍流模型

模型以海水-冰晶混合液为研究对象，考虑液体相和颗粒相湍流动能作用对流动过程的影响。海水-冰晶两相流在管道中的流动特性用混合相k-ε湍流模型描述，并通过壁面函数法解决壁面流动问题。k-ε湍流模型^[11]中的湍流动能方程为：

$ \frac{\partial }{{\partial {\rm{t}}}}\left( {{\rho _m}k} \right) + \nabla \cdot \left( {{\rho _m}k{\mathit{\boldsymbol{v}}_m}} \right) = \nabla \cdot \left( {\frac{{{\mu _{t.m}}}}{{{\sigma _k}}}\nabla k} \right) + {G_{k,m}} - {\rho _m}\varepsilon $

(7)

式中σ_k=1。湍流动能的耗散率方程可表示为：

$ \begin{array}{*{20}{c}} {\frac{\partial }{{\partial {\rm{t}}}}\left( {{\rho _m}\varepsilon } \right) + \nabla \cdot \left( {{\rho _m}\varepsilon {\mathit{\boldsymbol{v}}_m}} \right) = }\\ {\nabla \cdot \left( {\frac{{{\mu _{t.m}}}}{{{\sigma _\varepsilon }}}\nabla \varepsilon } \right) + \frac{\varepsilon }{k}\left( {{C_{1\varepsilon }}{G_{k,m}} - {C_{2\varepsilon }}{\rho _m}\varepsilon } \right)} \end{array} $

(8)

式中：

$ \left\{ {\begin{array}{*{20}{l}} {{\rho _m} = {\alpha _l}{\rho _l} + {\alpha _s}{\rho _s}}\\ {{\mathit{\boldsymbol{v}}_m} = {\alpha _l}{\rho _l}{\mathit{\boldsymbol{v}}_l} + {\alpha _s}{\rho _s}{\mathit{\boldsymbol{v}}_s}/{\alpha _l}{\rho _l} + {\alpha _s}{\rho _s}}\\ {{\mu _{t.m}} = 0.09\left( {{\alpha _l}{\rho _l} + {\alpha _s}{\rho _s}} \right){k^2}/\varepsilon }\\ {{G_{k,m}} = {\mu _{t.m}}\left( {\nabla {\mathit{\boldsymbol{v}}_m} + {{\left( {\nabla {\mathit{\boldsymbol{v}}_m}} \right)}^{\rm{T}}}} \right):\nabla {\mathit{\boldsymbol{v}}_m}} \end{array}} \right. $

(9)

1.5 相间作用力

溶液流动时，固相和液相之间的动量交换是由相间力引起，其相间作用力主要考虑湍流扩散力和拖曳力。其中湍流扩散力F_{td, l} 通过Burns等^[12]的模型获得；拖曳力表述为：

$ {F_{{\rm{drag}}{\rm{. }}l}} = {k_{sl}}\left( {{\mathit{\boldsymbol{v}}_s} - {\mathit{\boldsymbol{v}}_l}} \right) $

(10)

式中k_sl 表示动量交换系数:

$ {k_{sl}} = 3{C_D}{\rho _I}{\alpha _s}{\alpha _l}/4{d_S}\mid {\mathit{\boldsymbol{v}}_s} - {\mathit{\boldsymbol{v}}_l}\mid $

(11)

拖曳力系数C_D，表示为：

$ {C_D} = {\left( {0.63 + \frac{{4.8}}{{\sqrt {{\rm{R}}{{\rm{e}}_s}/|{\mathit{\boldsymbol{v}}_s} - {\mathit{\boldsymbol{v}}_l}|} }}} \right)^2} $

(12)

1.6 颗粒流动力学理论

固体颗粒的性质通过颗粒流动力学理论描述，其基本思想为把固体颗粒处理为致密气体分子。在国内外两相流研究中，Gidaspow^[13]提出的基于气体分子运动论的颗粒动力学理论模型被广泛使用，表达式为：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{3}{2}\left[ {\frac{\partial }{{\partial t}}\left( {{\alpha _s}{\rho _s}{\theta _s}} \right) + \nabla \cdot \left( {{\alpha _s}{\rho _s}{\mathit{\boldsymbol{v}}_s}{\theta _s}} \right)} \right] = \\ \left( { - {P_s}\mathit{\boldsymbol{I}} + {\tau _s}} \right):\nabla {\mathit{\boldsymbol{v}}_s} + \nabla \cdot \left( {{k_{{\theta _s}}}\nabla {\theta _s}} \right) + {\delta _{ls}} - {\gamma _{{\theta _s}}} \end{array} $

(13)

式中：$ ( - {P_s}\mathit{\boldsymbol{I}} + {\tau _s})$ : $ \nabla $v_s是粒子应力张量；k_θs是扩散系数；k_θs$ \nabla $θ_s 为扩散能；δ_ls 为相间转换能量；γ_θs是颗粒碰撞所消耗的能量。

式(13)中，颗粒温度θ_s 用于描述颗粒波动能量，其表达式为：

$ {\theta _s} = \frac{1}{3}v_s^\prime v_s^\prime $

(14)

式中v_s^′表示颗粒波动速度。

式(13)中，颗粒压力P_s的方程表示为:

$ {P_s} = {\alpha _s}{\rho _s}{\theta _s} + 2{\rho _s}(1 + {e_{ss}})\alpha _s^2{g_o}{\theta _s} $

(15)

式中：e_ss =0.9为粒子碰撞恢复系数；g_o是颗粒径向分布函数；δ_ls=-3k_slθ_s。

体积粘度ξ_s表征颗粒抵抗变形的能力，表示为：

$ {\xi _s} = \frac{4}{3}\alpha _s^2{\theta _s}{\rho _s}{d_s}\left( {1 + {e_{ss}}} \right)\sqrt {\frac{{{\theta _s}}}{{\rm{ \mathsf{ π} }}}} $

(16)

颗粒剪切粘度由粒间碰撞引起的颗粒粘性μ_{s. c}及动力粘度μ_{s. k}组成:

$ {{\mu _s} = {\mu _{s.c}} + {\mu _{s.k}}} $

(17)

$ {{\mu _{s.c}} = \frac{4}{5}\alpha _s^2{\theta _s}{\rho _s}{d_s}\left( {1 + {e_{ss}}} \right)\sqrt {\frac{{{\theta _s}}}{{\rm{ \mathsf{ π} }}}} } $

(18)

$ {{\mu _{s.k}} = \frac{{10{\rho _s}{d_s}\sqrt {{\theta _s}{\rm{ \mathsf{ π} }}} }}{{6\left( {3 - {e_{ss}}} \right)}}\left[ {1 + \frac{2}{5}\left( {1 + {e_{ss}}} \right)\left( {3{e_{ss}} - 1} \right){\alpha _s}{g_o}} \right]} $

(19)

1.7 壁面振动模型

如图 1为水平管内冰晶流动示意图，沿x轴为方向流动，在管壁外壁面上添加振动。振动具有简谐运动规律，将振动相位角设为0，则圆管振动过程位移、速度及加速度分别表示为^[14]:

$ \left\{ {\begin{array}{*{20}{l}} {y = A\sin (\omega t + \varphi )}\\ {v = {\rm{d}}y/{\rm{d}}t = A\omega \cos (\omega t + \varphi )}\\ {a = {\rm{d}}v/{\rm{d}}t = - A{\omega ^2}\sin (\omega t + \varphi )} \end{array}} \right. $

(20)

2 模型设置及验证 2.1 物理模型

物理模型如图 2所示，选取长度为1 000 mm，外径为25 mm，内24 mm的水平直管。数值模拟水平直管相关参数如表 1所示，材料为碳钢管，密度为7 850 kg/m³, 采用ANSYS ICEM对直管进行网格划分，管道两端采用O-Grid网格，将直管划分至75 411个六面体网格，数值计算结果符合网格无关性要求。通过C++语言编译振动条件UDF自定义动网格程序，将动网格用于壁面进行求解。海水-冰晶热物性参数与海水盐度有关，本研究海水盐度取为15‰。根据在北极实地航行的科学考察船“永盛轮”提供的资料和壳管式换热器设计规范^[15-16]，海水-冰晶热物性参数如表 2所示。

2.2 边界和初始条件

为了研究海水-冰晶两相流流动特性，本研究忽略直径较小的冰晶颗粒，颗粒直径统一为0.5 mm。在管道边界条件方面，计算域入口流体采用速度入口，由于振动形成非稳定流场，造成出口处压强和速度改变，因此采用自由流出条件。管壁面为无滑移速度，近壁面处采用标准壁面函数，忽略颗粒在流动过程中的损耗，管壁壁面热流量为0，绝热情况。采用相耦合SIMPLE算法及二阶迎风精度求解离散方程，将时间步长设为0.02 s，计算残差小于1.0×10^-4时为收敛状态。

2.3 模型验证

为了验证所建立的CFD数值模型的有效性，搭建了海水冰浆流动实验平台，如图 3所示为其实物图，相对应的海水冰浆制取流动实验测试系统原理如图 4。实验系统包括：冰浆的制取装置、冰浆的储存装置、冰浆输送系统、流动实验测试段。本实验通过对整个管道系统外壁包裹2 cm厚度的聚氨酯材料进行保温绝热，直管内海水冰浆流动在绝热条件下测试。测试段为长度1 m，管径为24 mm的直管段，通过无纸记录仪及差压变送器读取压降值，实验测量装置参数如表 3。

本实验保证冰浆运行温度在-1.7 ℃左右，冰浆样品如图 5所示。冰粒直径通过制冰机制取0.5 mm范围内，含冰率IPF控制约为15%，流速控制在1.5 m/s。

实验主要测试振动频率为20 Hz条件下，通过振动调节器改变振幅大小，测量冰浆流动过程中压降的结果，并与CFD双欧拉模型模拟结果相比较，如图 6所示。图 6可以看出，实验结果与仿真计算对于压降影响的变化趋于一致，且两者的最大相对误差在20%范围内，这由于数值模拟忽略粘性耗散与摩擦效应，未考虑实验中存在的相变等因素，在实验中随着振幅的增大，颗粒与海水的能量交换越多。由此，数值模拟方法适用于研究海水-冰晶两相流流动特性。

3 数值模拟结果与讨论 3.1 振动对海水-冰晶颗粒分布影响

直管出口处不同振动工况下冰晶颗粒体积分布如图 7所示，冰晶体积分数分布范围0.03~0.21，沿管道圆形截面垂直方向上，冰晶体积分数呈对称性分布。冰晶体积分布出现明显分层，冰晶颗粒多聚集在上壁面，这主要因为冰晶颗粒密度小于海水密度，在浮升力的作用下颗粒群出现上浮现象。图 7表示振幅为0.5 mm下，上壁面冰晶颗粒体积分布明显减少，由于振动频率增加，导致冰晶颗粒之间碰撞频繁，致使冰晶颗粒受到的浮升力小于冰晶颗粒碰撞产生的力。

图 8在频率为20 Hz下，由于振幅的增大导致冰晶颗粒振动势能增加，颗粒与壁面碰撞增强，冰晶体积分布逐渐减少。

3.2 振动对速度分布影响

如图 9(a)研究振动对轴线上冰晶颗粒速度分布，流速分布先增加再减小，最终趋于平缓。主要因为进口段扰动较大，冰粒与流体间相互作用增强，能量转换增多使得压降增加，造成流速迅速增加。随着冰粒与壁面的磨损作用及冰粒间的碰撞损耗，冰粒速度逐渐降低，在0.5 m处趋于平缓。与静态相比，随着振幅与频率的添加，流速有微小的减弱。主要因为振动对中心区域流速影响较小，流速基本趋于稳定。如图 9(b)，沿直管轴线截面湍流动能先逐渐增加，0.5 m处湍流动能达到较大值，主要由于入口段压降较大，颗粒与海水相互作用增强，导致轴向拖曳力增加。由于管内粘性摩擦与机械摩擦作用，流动损耗增加，湍流动能逐渐降低至稳定值。

图 10为不同振动工况下振动对出口截面y轴颗粒流速分布系数的影响，从近壁面到主流区流速逐渐增加且主流区流速分布系数基本不变。这由于靠近壁面粘性摩擦与机械摩擦作用，同时冰粒在上壁面聚集，导致颗粒间的相互作用力阻碍颗粒的流动。随着振动频率及振幅的增加，近壁面流速分布系数逐渐增大。这由于频率增加导致颗粒与壁面之间碰撞频繁，振幅的增加促使颗粒与壁面碰撞作用增强，导致机械振动能量转换为动能的量增大，颗粒的速度增加。

3.3 流动阻力变化

海水-冰晶流动过程中，由于固液流体之间的相互作用力及进口处扰动，流体流动阻力不容忽略。流体的流动阻力采用摩擦系数f表示，摩擦系数由流体在管内流动过程中压力降表示为^[17-18]：

$ f = 2\frac{{\nabla p \cdot D}}{{L{\rho _c}{u^2}}} $

(21)

式中：D和L分别是直管直径和长度；ρ_c是冰晶颗粒密度；u为流体流速。

流体流动阻力f与流道压降$ \nabla $P成正比关系，流动阻力的变化可通过分析压力降得出。

图 11为不同振动工况下压降变化趋势。随着振动频率及振幅的增加，直管进出口的压降逐渐增大。由于振幅和频率的增大，颗粒的机械振动能量增加，造成颗粒与壁面及颗粒间的相互作用增大，致使更多的振动势能转换为颗粒动能，通过式(21)压降与速度平方成正比，压降增加。

4 结论

1) 冰晶两相流在管道中流动由于浮升力作用产生明显分层现象，上壁面附近冰晶颗粒积聚。随着频率从0增到50 Hz、振幅从0增加到1.2 mm，管壁附近区域流体的扰动增强，造成上壁面附近的冰晶颗粒体积分布逐渐减少。

2) 由于颗粒与壁面间粘性摩擦与机械摩擦作用，壁面附近速度小于主流区速度。主流区流速及湍流动能受振动影响较小，在近壁面振动增强引起振动能量增加，转换成动能量增大，流速与湍流动能增大。

3) 通过数值模拟与实验结果先比较，验证了数值模拟的有效性。在频率为20 Hz条件下，振幅从0增加到1.2 mm，海水冰晶两相流的流速增大，造成管道进出口压降增大。