螺旋管内单相流动周向非均匀传热现象的数值模拟

2020-08-03顾汉洋叶亚楠

上海交通大学学报 2020年7期

关键词：螺旋管壁温离心力

王瑞，肖瑶，顾汉洋，叶亚楠

(上海交通大学核科学与工程学院，上海 200240)

符号说明

a—螺旋管内径， mm

c—螺旋管螺旋直径，mm

Cp—比定压热容， kJ/(kg·℃)

F—合力，N

g—重力加速度，kg/(m·s2)

h—换热系数，W/(m2·℃)

k—导热系数， W/(m·℃)

Nu—努塞尔数

Nuθ—θ处的局部努塞尔数

q—热流密度，kW/m2

Re—雷诺数

T—局部内壁温度， ℃

Tb—平均流体温度，℃

Tw—平均内壁温度，℃

u—流速， m/s

α—螺旋管螺旋升角，(°)

β—离心力与重力合力与竖直方向的夹角，(°)

θ—截面周向角度，(°)

Θ—无量纲温度

ρ—密度， kg/m3

φ—重力加速度与离心力加速度之比

螺旋管由于其结构紧凑、易于制造、换热效率高等优点被广泛地应用于食品工业、核工业、废热回收、制冷、航空航天和许多其他工业场景[1-3].在核工业中，螺旋管是近年来广受关注的小型模块化反应堆蒸汽发生器的理想选型.在设计螺旋管式蒸汽发生器时，充分了解螺旋管的传热特性才能保证蒸汽发生器换热能力与堆芯功率匹配使反应堆安全运行.因此，对螺旋管内流动传热的研究具有重大意义.许多研究指出螺旋管内流体存在复杂的流型[4].螺旋管内流体在离心力的作用下在横截面上产生了二次流.二次流的形态随边界条件的变化而变化，使得流体的速度分布不均，强化了螺旋管内的传热[5-7].由于实验方法不能够详尽地描述螺旋管内的局部单相流动换热特性，所以许多相关研究采用了数值模拟方法.这些研究中最常使用的湍流模型是k-ε模型和Reynolds应力模型.马越等[8]使用Reynolds应力模型模拟了高温气冷堆中氦气横掠外壁面时螺旋管截面上的壁面温度分布，发现螺旋管壁面最高温度位于截面顶部与内侧之间.史建新等[9]使用重整化群(RNG)k-ε模型模拟了特定结构的螺旋管内单相流体的周向流动与传热分布，发现螺旋管外侧流速大、内壁温低、换热系数高.Jayakumar等[10]使用了Realizablek-ε模型清晰地描述了螺旋管分别沿轴向和径向的局部努塞尔数波动，同时利用计算结果研究了螺旋管结构参数对传热的影响，并验证了已有的传热经验关系式的有效性.

综上所述，对螺旋管截面周向的传热分布研究比较初步，如文献[8]发现了螺旋管最高壁温位于顶部及内侧之间，却没有指出其具体位置及变化规律；如文献[10]给出了螺旋管周向传热分布的特征却没有从受力角度去解释成因.本文对螺旋管的周向局部特性做了更进一步的探讨，系统地研究了螺旋管传热周向分布的特征，总结了管壁周向温度分布随加速度比、螺旋直径、螺旋升角、螺旋管水力学直径等参数的变化规律及成因.

1 数值模拟

1.1 数学物理模型

所模拟的螺旋管物理结构参数如表1所示.应用计算流体力学(CFD)分析软件CFX16.0进行模拟计算，湍流模型采用Reynolds应力模型.文献[8]指出，Reynolds应力模型更多地考虑了旋转流动的特性，更适合于计算旋转流动中的流体；而常用的k-ε模型则采用各向同性的湍动黏度计算湍流应力，在计算旋转流动方面偏差较大.同时，文献[8]还验证了Reynolds应力模型的计算结果与经验关系式最为吻合，故本文采用Reynolds应力模型进行计算.

在CFX计算中advection scheme和turbulence numerical选项均设置为high resolution，即动量方程对流项和湍流输运方程对流项的离散格式均采用比较精细的离散格式.当质量、动量、能量等相关量的均方根(RMS)残差为1×10-7时，判断计算收敛.边界条件设置如下：入口为恒定质量流速入口，并设定恒定入口温度；出口为定压出口；壁面为无滑移的光滑壁面，壁面上设置恒定均匀热流密度加热；设重力加速度大小为9.8 kg/(m·s2)，方向为竖直向下.以水为工作介质，水物性采用IAPWS IF97的变物性程序包进行计算，计算区域水保持单相.计算涉及到的各工况边界条件如表2所示.

表2 工况边界条件Tab.2 Boundary conditions

1.2 网格的划分

分别对不同结构的螺旋管使用ICEM CFD软件创建结构化六面体网格.为了保证计算精度，进行网格的独立性实验.在工况1的边界条件下由不同数量的网格计算出的某截面(截面Tb=141.2 ℃)周向最高内壁温，如图1所示.其中：N为体积单元网格数量；Tm为最高内壁温度.由图1可见，随着网格数量的上升，最高内壁温度逐渐增大.当体积单元网格数量达到7×106后最高内壁温度基本保持不变，故认为此时的网格符合独立性要求.经过独立性验证后的各螺旋管体积单元网格数量如表3所示，

图1 螺旋管1的网格独立性实验Fig.1 The grid independence experiment of Helical Pipe 1

表3 各螺旋管的网格数量Tab.3 Grid numbers of different helical pipes

优化后的螺旋管1网格如图2所示.螺旋管1横截面上总节点数为 1 464，第一层网格质心到壁面的无量纲距离Y+接近30，符合Reynolds应力模型的要求.

图2 螺旋管1的网格结构Fig.2 Grid structure of Helical Pipe 1

2 计算结果的实验验证

在计算换热系数时需要用到截面平均流体温度.由于流体的物性与温度有关，采用文献[10]的方式计算截面平均流体温度：

(1)

式中：A为平均流体轴向截面.

截面换热系数如下式：

式中：Tθ为θ处的局部内壁温；θ1，θ2为任意两个截面的轴向角度位置.

图3 工况1管道沿程计算和实验结果的对比Fig.3 Comparison of calculated and experimental results along the flow direction for boundary Condition 1

图4 螺旋管截面上T和T′的对比Fig.4 Comparison of T and T′ on helical pipe cross section

图5 换热系数计算值和实验值比较Fig.5 Comparison of calculated and experimental heat transfer coefficient

3 结果分析与讨论

截面上周向努塞尔数的计算式为

(4)

截面努塞尔数计算式为

(5)

采用归一化努塞尔数Nuθ/Nu讨论局部换热情况.螺旋管截面的周向传热分布受重力和离心力的共同影响，故提出以下无量纲参数：

(6)

即φ可用于表征重力与离心力合力的方向.此外将壁温也进行了无量纲化，无量纲的温度形式为

(7)

以下将就Θ与Nuθ/Nu的周向分布随φ和螺旋管几何结构参数的变化情况分别进行讨论.

3.1 φ对截面周向传热分布的影响

首先研究了φ变化较小的情况下周向传热分布的变化情况，选取工况1管道沿程的4个截面进行比较.工况1螺旋管截面传热分布情况如图6所示.由图6可知，由于这4个截面沿流动方向的流体温度递增、密度递减、黏度递减，质量流速恒定时u和Re均递增，φ由12.66递减至11.79.此时，Θ和h在周向上的变化曲线形状基本不变：螺旋管内壁温从θ=0°开始缓慢下降，在90°左右下降到最低，之后又开始缓慢上升，θ=225°之后温度开始急剧上升，在θ=285°～300°之间温度达到最高.总体而言，低温区温度变化比较平缓，高温区温度变化剧烈；与Θ变化趋势相反，h在高换热区域，即θ=45°～135°区域变化较为平缓，在θ=80°～100°之间达到最大，换热性能最强，之后又开始缓慢下降，θ=225°之后h急剧下降，在θ=290°～300°之间达到最低.随着Re的上升，离心力略有增大，重力和离心力的合力方向发生改变，使得最高壁温位置由θ=295°向θ=285°方向作小角度的移动，整体而言变化十分微小.因此可以认为，φ变化较小时截面的周向传热分布特性基本不变.此外，Re越大，尽管h的最低值基本不变，但h的最高值及平均值均有所上升，即Re越大，整体换热性能越强.

图6 φ微小变化时的无量纲温度分布和换热系数分布Fig.6 Dimensionless temperature distribution and heat transfer coefficient distribution with small changes in φ

φ变化较大的情况下，周向传热分布的变化情况如图7所示.分别选取工况1、2、3、4的某横截面进行研究.这4个工况的物理模型均为螺旋管1，但因流速差异较大，其φ值有明显的差异.不同φ对应的最高壁温点位置等分布特征如表4所示.综合表4和图7可以看出，无重力作用时，最高壁温点在θ=270° 左右即最内侧，最低壁温点在θ=90° 左右即最外侧.此时，螺旋管内流体在截面上主要受到离心力和压强梯度力的作用.离心力指向管外侧，压强梯度力指向管内侧.离心力的大小与流体轴向速度的平方呈正比，在螺旋管截面中心区的流体流速较高，且离心力大于压强梯度力，管内侧壁面处的高温流体经过管道中心区向管外侧流动，并与主流的低温流体混合后降温，而贴近管壁处的流体由于黏性力的影响，轴向速度变小，离心力小于压强梯度力，到达外侧壁面的低温流体在压强梯度力驱动下分两支沿两侧管壁流回管内侧并由两侧管壁加热后升温，如此循环往复形成了截面上下两个Dean涡，即图8中的二次流，也造成了θ=90°位置流体温度及壁温低，而θ=270°位置流体温度及壁温高的现象[11-13].而当重力相对作用逐渐增强时，由图7及表4可知，最高壁温点由内侧开始向顶部移动.当φ=63.96时，最高壁温点甚至移动到了θ=338°处，同时高换热区域也逐渐从外侧向管底部移动，并且该区域逐渐变宽，温度变化更加缓和.同时，Nuθ/Nu变化幅度也逐渐变小，这是因为φ增大时，螺旋管内流速降低、离心力变小、二次流强度变弱，截面周向传热不均匀性程度降低.

图7 不同φ下横截面上的无量纲温度总体分布和Nuθ/Nu分布Fig.7 Dimensionless temperature distribution and Nuθ/Nu distribution on cross sections at different φ values

图8 无重力作用下的螺旋管二次流Fig.8 Secondary flow of helical pipes without gravity

图7对应曲线的详细温度分布及二次流情况如图9所示，其中Tbl为局部Tb.图9(a)中φ=63.96，即重力占主导作用，这时截面中心的流体从θ=340°方向向大约θ=150°方向流动，θ=150°一侧流速最大、温度最低，而θ=340°方向温度最高，整个温度分布呈现出明显的分层.从图9(a)～(d)离心力作用逐渐增强，截面中心二次流方向逐渐由倾斜变为水平，温度分布的分层情况也和二次流流动状况相符合，与图7的规律一致.此外由表4可知，当有重力作用时，壁温最高点与壁温最低点的角度差在θ=185°～220°之间，超过θ=180°.结合图9的二次流流动情况可以发现，重力对内外侧流体二次流的影响不同，图9(d)中仅有离心力作用，截面中心的二次流水平流动，因而最高壁温和最低壁温位置相差接近180°；图(b)和(c)中重力作用增强，重力的浮升力效应使得内侧的热流体上浮，该侧的中心水平二次流变成了由θ=300°～320° 位置斜向下的流动，而外侧流体的流动方向也由水平方向略微向下偏移，但由于内侧流体二次流流动方向在重力作用下偏移更大，导致了最高壁温位置与最低壁温位置的角度差超过180°，这反映出重力对内外侧二次流不同的影响.

表4 图7中不同φ下的周向传热特征

图9 不同φ下横截面上的温度分布及二次流情况Fig.9 Temperature distribution and secondary flow on cross section at different φ values

3.2 螺旋管结构参数对截面周向传热分布的影响

为探究a对截面温度分布的影响，在工况7和工况8的计算结果中分别选取了φ=22.25，20.22两个截面进行比较.由表2可以看出，工况7和工况8的物理模型除a不同外，其他参数均相同.不同水力学直径下横截面的温度及Nuθ/Nu的分布如图10所示.由图10可知，a对螺旋管周向传热分布没有明显的影响，相同φ下工况7与工况8的Nuθ/Nu分布曲线基本一致.

图10 不同水力学直径下横截面的温度及Nuθ/Nu分布Fig.10 Temperature and Nuθ/Nu distribution of cross sections at different hydraulic diameters

为探究c对截面温度分布的影响，在工况5和工况7的计算结果中分别选取了φ=20.23，19.27两个截面进行比较.不同螺旋直径下截面温度及Nuθ/Nu的分布如图11所示.由表2可以看出工况5和工况7的物理模型除c不同外，其他均相同.图11的结果表明，当φ一定时，c变大，周向壁温峰值和谷值位置不变，而Nuθ/Nu波动幅度明显变小.结合上文所述的二次流的形成机理对此现象解释如下：当φ不变时，重力和离心力的合力方向不变，二次流动方向不变，所以周向壁温峰值和谷值位置不变；当φ不变时，离心力大小不变，螺旋管截面中心区的二次流主要在离心力作用下产生，因此强度变化不大；当c增大即螺旋管曲率减小时，截面上压强梯度力变小，贴近壁面处的二次流由压强梯度力克服离心力产生，这部分二次流强度明显减弱，导致了二次流强度整体变弱，截面周向传热不均匀程度减小，Nuθ/Nu波动幅度减小.

图11 不同螺旋直径下截面温度及Nuθ/Nu分布Fig.11 Temperature and Nuθ/Nu distribution of cross sections at different hydraulic diameters

为探究α对截面温度分布的影响，在工况7和工况6的计算结果中分别选取了φ=23.53，17.5的两个截面进行了比较.从表2可以看出工况7和工况6的物理模型除α不同外其他均相同.不同螺旋升角下横截面的温度及Nuθ/Nu分布如图12所示.图12结果表明，α对螺旋管周向传热并没有明显的影响，相同φ下工况7与工况6的Nuθ/Nu分布曲线基本一致.螺旋管受力分析如图13所示，其中F为向心力与重力的合力.由图13可知，重力与横截面1-1的夹角为α，故重力加速度在横截面1-1上的分量为gcosα，对于螺旋线运动，离心力方向为水平背离螺旋线中心轴方向，又知离心力方向与螺旋管轴线切线和重力方向垂直，可以得出离心力与重力在横截面1-1上的分量垂直，故离心力方向如图12中的横截面1-1所示，其加速度大小为2u2/c.因此可以得出F与离心力的夹角β满足：