渔获物冰浆保鲜流动传热分析

赵新颖，黄温赟，黄文超，吕续舰

(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种流动速度的工况进行对比可以发现，适当提高流动速度可以增大渔获物表面冰颗粒的体积分数，从而增加表面热量交换。当流动速度增大时，渔获物表面的冰颗粒数量一开始明显增大，当速度增大到一定水平时，渔获物表面上冰颗粒体积分数十分接近，此时渔获物上表面传热量增幅很小。

□