1 前言

缠绕管式换热器属于一种管壳式换热器，其换热管束以螺旋缠绕的形式排列在壳体内部，具有结构紧凑、耐高压、热应力小、换热系数高、流体适应性强、可实现多种介质同时换热等优点，广泛应用于制药、食品、石油、化工、冶金、电力和纺织等行业^[1-3]。很多学者对圆形光管缠绕管换热器进行了研究。LU等^[4-5]实验研究并模拟了恒定功率电加热条件下壳侧空气流动及换热。NEERAAS等^[6-7]对多种实验工质单相换热和两相混合换热进行了实验研究，拟合了努塞尔数Nu和摩擦系数f的经验关联式。ZENG等^[8]模拟了多层缠绕管换热器壳侧流动和换热，使用田口算法评估了主要结构参数的影响，得出了Nu和f的多参数关联式；TANG等^[9]实验并模拟了四层缠绕管换热器壳侧空气流动和换热，实验结果与Messa的压降关联式偏差12.82%。REN等^[10-11]模拟了左右晃动或旋转搅动状态下，壳侧液化天然气流动和换热，模拟与实验结果的压降偏差为16.81%。FERNG等^[12]研究了一根缠绕管的换热器，分析了迪恩Dean数和管间距的影响，与实验结果吻合良好。王斯民等^[13]模拟分析了关键几何参数的影响，使用多目标遗传算法优化了结构。WANG等^[14]模拟发现交错分布的垫条能够增加流体扰动、增强换热；垫条数量的影响不大。以上研究为异形管缠绕管换热器的研究奠定了基础。RAINIERI等^[15]实验研究了螺纹管和光滑管绕管曲率和螺纹结构的影响，结果表明低Dean数下，绕管曲率的影响占优；在高Dean数下，螺纹结构占优。NAPHON^[16]对绕管外带翅片和不带翅片的缠绕管换热器进行了实验研究。GUPTA等^[17]实验研究了外翅片管(管内有波纹)不同于光管的压降关联式。尹接喜等^[18]模拟了并管结构缠绕管换热器。马飞^[19]建立了椭圆管和三叶管的缠绕管换热器模型，结果表明三叶管的综合性能最好。田杨等^[20]模拟发现与圆形管缠绕管换热器相比，水滴形管能大幅度降低壳侧压降。

异形管缠绕管换热器的数值模拟研究集中在改变圆管截面大尺寸方面，例如：三叶管、椭圆管，对改变圆管截面微小尺寸进而强化换热的研究鲜有报道。本文基于流动传热强化机理，提出了基于开槽强化管的缠绕管换热器，使用Fluent 18.0数值分析了壳侧空气流动、换热和综合性能随凹槽数量和深度的变化规律，并与圆形光管缠绕管式换热器进行了对比。

2 计算模型及数值方法 2.1 几何模型的建立

本文建立了缠绕管式换热器的三维模型，其主要结构参数如图 1所示，芯桶的外部直径${D_{{\rm{cylin}}}}$34 mm；和壳体的内部直径${D_{{\rm{shell}}}}$66 mm；换热管高度$H$100 mm；同层管间距$l$2.5 mm；换热管的外径，${D_{\rm{t}}}$8 mm；缠绕角度$\alpha $9.5°。为了减少进出口回流的影响，进出口有50 mm延长段。缠绕管上有半圆形凹槽，与管的中心线平行。不同凹槽数量、凹槽深度的缠绕管截面图如图 2所示。

2.2 网格划分及数值方法

由于管外壁凹槽深度很小，将凹槽部分独立划出后进行网格划分。凹槽内的计算域采用结构化网格，其余部分采用非结构化网格划分，如图 3所示。对凹槽数量为2，凹槽深度为0.3 mm几何模型的网格进行了独立性验证，网格数量为16 056 225 ~ 35 456 358，结果如图 4所示。当网格数量大于24 133 963时，对流换热系数h的变化小于1.5 %，压降Δp的变化小于1.0 %。所以取网格数量29 512 101满足要求。不同组合下计算域的网格单元数为2 500万~3 600万。

本文的计算工质为空气，在换热器壳侧单相换热，物性不变。入口边界条件是速度入口，为3.0、3.5、4.0、4.5和5.0 m·s^-1，空气进口温度为300 K，出口边界条件为压力出口；换热管壁温恒定，为330 K，其余壁面均为无滑移绝热壁面；壁面函数为默认的标准壁面函数，压力和速度耦合采用SIMPLE算法，动量和能量方程均为二阶迎风格式，各项残差小于为1×10^-6时，认为收敛。选择RNG k-$\varepsilon $湍流模型来计算，其基本控制方程包括连续性方程、质量守恒方程和能量守恒方程^[21]，笛卡尔坐标系下，表达式如下：

连续性方程：

$\nabla \mathit{\boldsymbol{u}} = 0$

(1)

动量守恒方程

$\nabla (\rho u\mathit{\boldsymbol{u}}) = \frac{\partial }{{\partial x}}(\eta \frac{{\partial u}}{{\partial x}}) + \frac{\partial }{{\partial y}}(\eta \frac{{\partial u}}{{\partial y}}) + \frac{\partial }{{\partial z}}(\eta \frac{{\partial u}}{{\partial z}}) - \frac{{\partial p}}{{\partial x}}$

(2)

$\nabla (\rho v\mathit{\boldsymbol{u}}) = \frac{\partial }{{\partial x}}(\eta \frac{{\partial v}}{{\partial x}}) + \frac{\partial }{{\partial y}}(\eta \frac{{\partial v}}{{\partial y}}) + \frac{\partial }{{\partial z}}(\eta \frac{{\partial v}}{{\partial z}}) - \frac{{\partial p}}{{\partial y}}$

(3)

$\nabla (\rho w\mathit{\boldsymbol{u}}) = \frac{\partial }{{\partial x}}(\eta \frac{{\partial w}}{{\partial x}}) + \frac{\partial }{{\partial y}}(\eta \frac{{\partial w}}{{\partial y}}) + \frac{\partial }{{\partial z}}(\eta \frac{{\partial w}}{{\partial z}}) - \frac{{\partial p}}{{\partial {\rm{z}}}}$

(4)

能量守恒方程

$\nabla (\rho \mathit{\boldsymbol{u}}T) = \frac{\partial }{{\partial x}}\left( {\frac{\lambda }{{{c_p}}}\frac{{\partial T}}{{\partial x}}} \right){\rm{ + }}\frac{\partial }{{\partial y}}\left( {\frac{\lambda }{{{c_p}}}\frac{{\partial T}}{{\partial y}}} \right){\rm{ + }}\frac{\partial }{{\partial z}}\left( {\frac{\lambda }{{{c_p}}}\frac{{\partial T}}{{\partial z}}} \right){\rm{ + }}{S_T}$

(5)

湍动能k方程

$\frac{\partial }{{\partial t}}\left( {\rho k} \right) + \nabla \left( {\rho k\mathit{\boldsymbol{u}}} \right) = \frac{\partial }{{\partial x}}\left[ {\left( {{\alpha _k}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial k}}{{\partial x}}} \right] + \frac{\partial }{{\partial y}}\left[ {\left( {{\alpha _k}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial k}}{{\partial y}}} \right] + \frac{\partial }{{\partial z}}\left[ {\left( {{\alpha _k}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial k}}{{\partial z}}} \right] + {G_k} - \rho \varepsilon $

(6)

湍动能耗散$\varepsilon $方程

$\frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \nabla \left( {\rho \varepsilon \mathit{\boldsymbol{u}}} \right) = \frac{\partial }{{\partial x}}\left[ {\left( {{\alpha _{_\varepsilon }}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial \varepsilon }}{{\partial x}}} \right] + \frac{\partial }{{\partial y}}\left[ {\left( {{\alpha _{_\varepsilon }}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial \varepsilon }}{{\partial y}}} \right] + \frac{\partial }{{\partial z}}\left[ {\left( {{\alpha _{_\varepsilon }}{\eta _{{\rm{eff}}}}} \right)\frac{{\partial \varepsilon }}{{\partial z}}} \right] + C_{1\varepsilon }^*\frac{\varepsilon }{k}{G_k} - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k}$

(7)

式中：u为速度的矢量形式，m·s^-1；u、v、w分别为x、y、z方向的速度，m·s^-1；c_p为流体的定压比热容，J·kg^-1·K^-1；p为压力，Pa；T为温度，K；ρ为流体的密度，kg·m^-3；λ为流体的导热系数，W·m^-1·K^-1；S_T为能量方程的广义源项；G_k为平均速度梯度引起的湍动能产生项；η为流体的动力黏度，Pa·s；η_tur为流体的湍动黏度，Pa·s；η_eff为流体的有效黏度，Pa·s；${\eta _{{\rm{eff}}}}{\rm{ = }}\eta {\rm{ + }}{\eta _{{\rm{tur}}}}$；${\eta _{{\rm{tur}}}} = \rho {c_\eta }\frac{{{k^2}}}{\varepsilon }$；${C_\eta }{\rm{ = 0}}{\rm{.08}}{{\rm{4}}_{}}{\rm{5}}$；${\alpha _k}{\rm{ = }}{\alpha _\varepsilon }{\rm{ = 1}}{\rm{.39}}$；$C_{1\varepsilon }^* = {C_{1\varepsilon }} - \frac{{\eta \left( {1 - {\eta \mathord{\left/ {\vphantom {\eta {{\eta _0}}}} \right. } {{\eta _0}}}} \right)}}{{1 + \beta {\eta ^3}}}$；$\eta = {\left( {2{E_{ij}}{E_{ij}}} \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}\frac{k}{\varepsilon }$；${\eta _0} = 4.377$；$\beta = 0.012$；E_ij为流体的时均应变率，${E_{ij}} = \frac{1}{2}\left[ {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right]$；${C_{{\rm{1}}\varepsilon }}{\rm{ = 1}}{\rm{.42}}$；${C_{{\rm{2}}\varepsilon }}{\rm{ = 1}}{\rm{.68}}$。

基于计算流体动力学(computational fluid dynamics，CFD)原理^[21]，可得到流体域范围内的速度，温度和压力。h和Δp可使用CFD-Post 18.0软件按式(8)和(9)处理：

$h = \frac{{{c_p}\dot m({T_{{\rm{outlet}}}} - {T_{{\rm{inlet}}}})}}{A}$

(8)

式中：$\dot{m} $为水的质量流速，kg·s^-1，T_outlet、T_inlet分别为进出口平均温度，K，A为换热总面积，m²。

$\Delta p = {p_{{\rm{inlet}}}} - {p_{{\rm{outlet}}}}$

(9)

式中：p_inlet和p_outlet分别为进出口平均压力，Pa。

2.3 模型的验证

将数值计算得到的Nu和单位长度压降Δp_m与文献[7]中实验所得的数据进行对比，如图 5所示。工质为甲烷气体。努塞尔数Nu和单位长度压降的平均偏差分别为5.67 %和7.90 %。换热管壁外其余壁面绝热导致模拟的Nu偏高；气体物性不变的假设，忽略了气体被加热时密度减小，导致模拟得出的压降值偏高。

3 结果和讨论 3.1 凹槽数量对流动和换热的影响

图 6~8分别显示了换热管缠绕角度为9.5°，凹槽深度为0.3 mm时，不同的壳侧入口速度下，h、Δp和换热器综合性能指标(performance evaluation criteria，PEC=Nu·f^1/3)随凹槽数量的变化。如图 6和7所示，随着凹槽数量的增加，h和Δp都有所增加，且增加趋势逐渐变缓。以壳侧入口速度3 m·s^-1为例，当凹槽数量为2时，h比圆形光管增大了11.86 %，Δp增大了17.75 %；当凹槽数量为6时，与圆形光管相比h增大了19.46 %，Δp增大了29.86 %。其余速度条件下，增加效果相当。这是因为凹槽开在传热管的两侧，凹槽的切线方向与流动方向接近垂直，增加了流体扰动，能够破坏沿流动方向发展的边界层，使热量通过涡流扩散的方式得到有效的传递，从而增大了对流换热系数。同时，流体扰动增强，湍动能增大，流体的内摩擦力变大，分子之间的碰撞更剧烈，造成了更大的能量损失，因此压降增大。

采用PEC来比较圆形管与凹槽管形的缠绕管式器壳侧的流动和换热性能，综合考虑了压降和换热，解决了换热增强和压降损失增大之间的矛盾，合理给定了两者之间的比重，能够评价换热器的综合性能。如图 8所示，随凹槽数量的增加，PEC在有2个凹槽时，比圆形光管先是有明显提高，再增加凹槽数量，PEC保持基本不变，即换热器的综合性能与凹槽数量无关，2个凹槽数量就可达到与更多凹槽数量一样的增强效果。以速度为3 m·s^-1为例，凹槽数量为2的情况下，PEC相对圆形光管增加了3.82%。在不同速度工况下，2个凹槽的缠绕管换热器相对圆形光管PEC增加率平均为4.21 %。这是由于凹槽数量增多时，对流换热系数得到增强的同时，压降也增大了，综合考虑两方面因素，性能并没有得到进一步提高。

3.2 凹槽深度对流动和换热的影响

如图 9所示，显示了凹槽数量为2，缠绕角度为9.5°时，h随凹槽深度在不同壳侧入口速度下的变化。从图 9可以看出，h随凹槽深度的增加先增大后减小，峰值在凹槽深度为0.3 mm时取得。在凹槽深度为0.3 mm时，在不同速度工况下，h比圆形光管的平均增长了11.79 %。这是因为，随凹槽深度增大，一方面流体扰动增强，流体的湍动能增大，流体分子剧烈碰撞，热量通过涡旋扩散的方式快速传递，对流换热系数增大，另一方面，在凹槽过深时，会引起流体在凹槽内部的滞流，热量通过一部分导热的方式来传递，形成了热阻，传热恶化。在凹槽深度较浅时，流体扰动占主导，换热增强，在进一步加深凹槽时，流体滞流引起的传热恶化占比越来越大，减小了对流换热系数。

如图 10显示了凹槽数量为2，缠绕角度为9.5°时，Δp随凹槽深度在不同壳侧入口速度下的变化。由图可以看出，Δp随凹槽深度先增大后维持不变，因为在凹槽深度为0.2、0.3 mm和0.4 mm情况下，Δp的偏差为0.05 %~1.41%，认为相等。不同速度工况下，凹槽深度为0.2 mm时，Δp比圆形光管平均增长了16.61%。由于当槽深度过大时，流体在凹槽内滞流后，流体在流经凹槽时，并没有得到更剧烈的扰动，湍动能的提高微弱。

PEC随凹槽深度的变化如图 11所示，可以看出，PEC随凹槽深度的增大先增大后减少，存在峰值，在凹槽深度为0.3 mm时取得。凹槽深度为0.3 mm时，不同速度工况下，相对圆形光管PEC增加率平均为4.21 %。值得注意的是，凹槽深度为0.4 mm时的增强效果不及0.1 mm时的增强效果，可见凹槽深度并不是越大越好。随凹槽深度增大，流体滞流在凹槽内部，使得换热恶化，压降损失维持不变，因此，缠绕管式换热器的综合性能变差。

4 结论

本文建立了基于开槽强化管的的缠绕管式换热器的三维模型，在圆形光滑缠绕管两侧沿管中心线开了不同数量和不同深度的半圆形凹槽，采用数值模拟的方法研究该种换热器壳侧空气流动、换热和综合性能随凹槽数量和凹槽深度的变化规律，与圆形光管缠绕管式换热器进行了对比。得出了以下结论：

(1) 与圆形传热管相比，凹槽的存在切断了边界层，增加了扰动，增强换热，换热器的综合性能有所提升。

(2) 当流体入口速度一定时，凹槽数量的增多，一方面增强了换热，另一方面增加了阻力压降；从换热器的综合性能来看，PEC并不随凹槽数量的增大而增大，而是维持基本不变，但相比圆形光管增强了约4.21 %。

(3) 对于凹槽深度的影响，当流体的入口流速一定时，对流换热系数和压降都随凹槽深度的增大先增大后减小，PEC在凹槽深度0.3 mm时取得最优值，此值与圆形光管相比，提高了约4.21 %。