刘献飞, 王恒，王方, 李志强, 朱彩霞, 张浩飞，王仕元
中原工学院能源与环境学院

       【摘  要】单螺杆膨胀机螺旋槽道内存在多种流型，掌握槽道内两相流动特性膜厚分布对减小气体泄漏，保证膨胀机高效稳定运行至关重要。本文采用VOF两相流数值模型，对垂直放置螺旋槽道内两相流动特性及膜厚分布进行了计算流体动力学分析。根据前期实验获得的螺旋槽内两相环状流的实验数据验证了数值模型的准确性，揭示了不同液相流体入射角度螺旋槽道内液相工质演变过程，理论分析不同液相流体入射角度下螺旋槽道外侧的液膜厚度分布规律，发现液相流体入射角度为0.9倍的螺旋升角时，液膜厚度分布最均匀。

       【关键词】螺旋槽道; 两相流; VOF模型; 流动特征；液膜厚度

       【基金项目】国家自然科学基金项目（51906265）；河南省高校科技创新人才支持计划（20HASTIT019）

引    言

       单螺杆膨胀机具有在低温下效率高以及对于工质的状态不敏感的优点，高压气体、饱和蒸汽、汽液两相或热液都可作为膨胀工质，大大拓宽了低温热源利用的温度范围。此外，单螺杆膨胀机所适用的功率范围在1kW～1000kW 范围内，克服了传统蒸汽轮机和燃气轮机不能做小的缺陷。因此，单螺杆膨胀机在工业节能、可再生能源利用、建筑节能和热泵系统等众多领域具有非常广阔的应用前景^[1]。目前关于单螺杆膨胀机的研究开发得到了国内外学者的高度重视，但由于单螺旋膨胀机涉及热力学、传热学、流体力学、机械设计、材料学及控制等多学科，理论研究难度很大。从国内外研究现状看，目前仍然侧重于膨胀机的加工和性能研究^[2,3]，对膨胀机螺旋槽道泄露间隙的液膜流动及分布特性研究数据缺乏、理论不足，而这正是制约膨胀机效率提高的关键因素^[4]。

       螺旋管内两相流动特性相对较小。Murai等人^[5]说明了空气-水两相流通过螺旋管时，离心加速度对流动结构分布的影响。Zhu等人^[6]基于双传感器电导探针研究了离心力对相分布的影响。Zhang等人^[7] 研究了螺旋管内油水分离机理，指出分离效率与旋转半径、压降、相密度差和液滴直径有关。Colombo等人^[8]采用欧拉模型对螺旋管内气液两相混合流动进行了数值模拟，得到了离心力对流场和相分布的影响。Mota和Pagano^[9]研究了不同流动参数条件和不同管道几何结构下螺旋管内气相分布。Jayakumar等人^[10]研究了螺旋管内两相流动时，螺旋管直径和螺距对热工水力特性的影响。Vashith和Nigam^[11]研究了两相流通过螺旋管时的速度分布、平均摩擦系数、相分布和入口长度。

       我们前期对螺旋矩形通道内环状流的流型演变和无量纲液膜厚度进行了实验研究^[12]，并且基于均相流模型^[13-15]，研究了不同结构参数对两相压降和相分布的影响，给出了螺旋通道内气液两相流体受离心力作用后相分布和液膜厚度。考虑到单螺杆膨胀机润滑油由油泵驱动，单独流入膨胀机间隙，用于润滑和密封^[1]。而目前对于润滑油单独进入间隙后，在螺旋矩形通道内的相演化行为、界面现象和膜厚分布的研究还很缺乏，本文利用VOF模型，分别定义气相和液相流动参数条件，给出了膨胀机螺旋槽道内液相演变规律及液膜厚度分布，这对如何调控膨胀机复杂截面螺旋槽道内液相流动、相界面特征及液膜厚度，避免高压气体向低压区域的漏泄损失，保证螺杆膨胀机高效、稳定运行具有至关重要的理论指导意义。 

1. 数值方法

       1.1 螺旋槽道物理结构示意图

       图1给出了螺旋槽道及其主要几何参数的示意图，矩形截面的宽度和长度分别为a=18mm、b=30mm。螺旋直径（从通道中心到轴测量）用D=195mm表示，两个相邻匝间的距离称为节距H=306mm，图中显示半节距h。螺旋线圈一圈的投影与垂直于轴线的平面形成的角度称为螺旋角β[^16]。流体入射角度（η）是指定义的流动方向与法向边界之间的夹角。螺旋槽道壁面最靠近线圈轴的一侧称为内侧，反之称为外侧。转角（θ）螺旋通道任何位置的局部截面与入口截面平面之间的角度。

图1 螺旋槽道及其主要几何参数的示意图

       1.2 数值方法

       采用基于流体体积法的三维CFD模型求解控制方程组，采用几何重建方法对气液界面进行重建。该方案使用分段线性方法表示流体之间的界面，使得两相流的计算在气体和液体之间具有清晰的界面^[12]。采用有限体积法求解控制守恒和边界条件，进行了数值模拟。前期我们已经使用了不同的湍流模型来预测螺旋管内流动的流体力学和热特性，为了更精确地预测螺旋管内两相流的流体力学特性[14]和提高了螺旋管内流动的精度，这里选择了重整化群k-ε湍流模型快速应变流^[17]。

       连续性方程:

           (1)

       动量守恒方程: 

           (2)

       式中，P、g、F、μ分别为流场中的压力、重力加速度、作用于系统的物体力和流动流体的粘度。

       上面给出的方程依赖于所有相的体积分数，由性质ρ和μ隐式表示。在两相系统中，每个单元的密度和粘度可以估计为：

           (3)
           (4)

       通过求解体积分数的标量输运方程，可以追踪两相之间的界面：

           (5)

       对于液膜流动，表面张力起着重要的作用。FLUENT中的表面张力模型是Brackbill等人提出的连续表面力（CSF）模型 ^[18]。在该模型中，表面张力对VOF计算的贡献是动量方程的源项，由表面张力产生的体积力可计算为^[19]：

           (6)

       式中，σ是表面张力系数，k是根据单位法向散度定义的自由表面曲率：

           (7)

       1.3 数值计算

       本文采用的求解程序是基于商用CFD软件FLUENT，网格是使用商业软件GAMBIT创建的。图2描绘了用于三维计算域的网格（六面体网格和楔形网格）。在FLUENT求解过程中使用非稳态求解器，整个过程物理性质保持一致的，操作环境考虑重力作用。边界条件是入口分别为气相和液相速度入口，出口采用压力出口边界，出口方向变量的扩散通量设为零。在槽道壁面设置了一个稳定、无滑移、无渗透的边界，并考虑了强化壁面作用处理。压力插值采用压力交错选项（PRESTO）格式，压力-速度耦合采用压力隐式算子分裂（PISO）格式，动量方程采用二阶迎风差分格式，体积分数方程采用Geo重构格式，动量、湍流动能和湍流耗散率采用二阶迎风格式。基于两相流的动态特性，进行了时间步长为0.0001s的非稳态数值模拟，连续性、速度、k和epsilon的收敛准则为1.0e-5。

       为了更好的获得螺旋槽道内的两相流特征，首先对各网点精度进了独立性验证。数值模拟的网格测试考虑了75104、152478、237162和291403个网格的计算节点，网格间距分别为0.0015、0.0012、0.001和0.0095。使用不同网格的预测值如表1所示，结果表明，0.001的区间节点比例可以为研究矩形螺旋槽内的液相占有率和液膜厚度提供可接受的统计精度。

（a）    （b）
图 2 螺旋槽道的网格特征（a）计算域（b）入口截面
表1 不同网格数下的预测结果比较

       1.4 数据处理

       图3a示出了从商业软件FLUENT 6.3的输出获得的横截面含液率分布的原始图像，为了更好地测量液膜厚度，采用数字图像处理算法将液膜轮廓与其他流动特征分离开来。通过设置含液率阈值，得到了两相流结构中液相的等表面图，在本文研究中液相的阈值固定为0.9^[20]。图3b示出了测量螺旋槽道外侧的液膜厚度的示意图，为了更好地从图片信息中获取液膜厚度，本文采用了50倍的放大因子进行测量。

（a）                                           （b）
图 3（a） 截面局部含气率分布；（b）从气相分离液膜轮廓后液膜厚度的测量

2. 结果与讨论

       2.1 模型验证

       VOF方法已广泛应用于界面位置的跟踪和界面现象的求解，并进行了大量的工作^[8,11,17]。对于螺旋矩形通道内的两相流，关于液膜厚度的研究成果并不多见。数值模型的验证是根据我们以前的研究^[12]中环形流的实验数据进行的。由于离心力的影响，在螺旋矩形通道的外侧突出了环形两相流的液膜厚度。为了验证数值解的准确性，根据实验研究，进行了该数值模拟的结构参数和两相质量流量的计算。图4给出了实验和数值模拟所得液膜厚度的比较，图中还显示了实验数据的不确定度。实验数据与模拟数据的偏差在3.79%以内，因此，本研究建议计算模型有足够的精确度。

图 4 实验和数值模拟的平均液膜厚度对比

       2.2 界面特征的演变规律

       单螺杆膨胀机槽道外侧的液相分布对防止气体泄漏具有重要作用，图5给出了螺旋槽道不同横截面处的含液率分布及演变规律。值得注意的是，当η=0时，由于螺旋槽道的节距产的扭曲和流体流动惯性的影响，流动会发生扭曲，在y/a=0位置时可以看到较大的持液率。随着流动的进一步发展，液体被离心力作用移向外侧。同时，在初始发展阶段，液体被节距产生的扭转力作用开始向y/a=1方向移动。在θ=60o角处，液体分布在水平中心线上接近对称，然后，由于二次流运动的惯性作用，液相的继续向y/a=1方向移动使得液相分布不对称^[13]，这种现象主要是由于螺旋通道入口的扭转效应引起的。随着流动的发展，扭转效应弱于曲率效应，对液体脉动分布的影响减弱。当η=0.9β时，根据数值模拟结果可以发现，扭转效应对螺旋槽液相分布的影响基本可以消除，主要受到离心力的作用，螺旋槽的外侧保持较高的持液率。

图5 螺旋槽道不同横截面处的含液率分布及演变规律（a）η=0；（b）η=0.9β

       2.3 液相入射角度对液膜厚度分布的影响

       从不同η下的等高线变化可以看出η对液相分布具有明显的影响，为了进一步揭示η对液相分布的影响，图6给出了不同η下的液膜厚度变化规律。从图6a可以看出，当η=0.80β时，液膜厚度在y/a=0位置明显较高，因为，这是由于η较小，流动仍然受到扭转力的明显影响。随着流动的不断发展，在y/a=0位置的液膜厚度先增大后减小，这是螺旋槽道节距产的扭曲和流体流动中的惯性，使得液体更容易在这一区域积聚。随着扭转力作用将流体从y/a=0的位置移到y/a=1的位置，使得y/a=1位置的液膜厚度逐渐增大。当η=0.85β时，由图6b可以看出，液膜厚度在y方向的不对称性显著减小。当η=0.90β时，图6c表示液膜厚度分布基本与均匀一致。当η继续增加到0.95β时，如图6d所示，在y/a=1区域出现较厚的膜厚分布。尽管η的值小于β，但由于螺旋槽道产生的较高的二次流导致液体向y/a=1移动^[21]。当η=β时，如图6e表示，y/a=1位置的液膜厚度明显高于y/a=0位置。对比不同θ下测量数据之间的薄膜厚度值可以发现，由于离心力始终在薄膜表面的法向上，随着流动的不断发展，外侧壁面薄膜平均厚度逐渐增加^[22]。

(a) η=0.80β                                            (b) η=0.85β

(c) η=0.90β

 (d) η=0.95β                                                      (e) η=β
图 6 不同进口角下液体厚度分布的变化

3. 结论

       本文借助VOF多相流模型，数值分析了两相流体的螺旋槽道内的流动特性，揭示了气液两相界面特征演变规律及液膜厚度分布规律，理论分析了液相入射角度对液膜厚度的影响规律，主要结论如下：

       1. VOF多相流模型能够很好地预测螺旋槽道内液相分布及液膜厚度变化规律，数值模拟结果和实验测量结果具有很好的一致性。

       2. 当液相入射角度为η=0时，螺旋槽道内液相分布呈现明显的波动性，当η=0.90β时，扭转效应对螺旋槽液相分布的影响基本可以忽略不计。 

       3. 不同入射角度下的液膜厚度对比可以发现，当η=0.9β时，螺旋槽道外侧壁面可以获得比较均匀的液膜厚度分布，当η=0.9β时，y/a=0位置区域的液膜厚度较高，反之y/a=1位置区域的液膜厚度较高。

参考文献

       [1] Wang Wei, Wu Yuting, Ma Chongfang et al, Preliminary experimental study of single screw expander prototype, Appl Therm Eng 2011;31:3684–3688. 
       [2] Ziviani D, Van M Den Broek, Paepe M De, Geometry-Based Modeling of Single Screw Expander for Organic Rankine Cycle Systems in Low-Grade Heat Recovery, Energy Procedia 2014;61:100–103.
       [3] Li Guoqiang, Lei Biao, Wu Yuting et al, Influence of inlet pressure and rotational speed on the performance of high pressure single screw expander prototype, Energy 2018;147:279–285 
       [4] Shen Lili, Wang Wei, Wu Yuting et al, Theoretical and experimental analyses of the internal leakage in single-screw expanders, Int. J. Refrigeration. 86 (2018) 273–281. 
       [5] Murai Y, Yoshikawa S, Toda S, Structure of air-water two-phase flow in helically coiled tubes, Nucl Eng Des 2006;236:94–106. 
       [6] Zhu Hongye, Li Zhaoxu, Yang Xingtuan et al, Flow regime identification for upward two-phase flow in helically coiled tubes, Chem Eng J 2017;308:606–618.
       [7] Zhang Jun, Guo Jun, Gong Daotong et al, An investigation on oil-water separation mechanism inside helical pipes, J Hydrodyn Ser B 2006;18:343– 347.
       [8] Inzoli F, M Colombo, Cammi A et al, CFD study of an air–water flow inside helically coiled pipes, Prog Nucl Energy 2015;85:462–472. 
       [9] da Mota FRM, Pagano DJ, Simulation and experimental study of phase segregation in helical pipes: A new method for flow conditioning, Flow Meas Instrum 2014;35:99–108. 
       [10] Jayakumar JS, Mahajani SM, Mandal JC et al, Thermal hydraulic characteristics of air-water two-phase flows in helical pipes, Chem Eng Res Des 2009;88:501–512. 
       [11] Vashisth S, Nigam KDP, Prediction of flow profiles and interfacial phenomena for two-phase flow in coiled tubes, Chem Eng Process Process Intensif 2009;48:452–463. 
       [12] Liu Xianfei, Xia Guodong, Yang Guang, Experimental study on the characteristics of air–water two-phase flow in vertical helical rectangular channel, Int J Multiph Flow 2015;73:227–237.
       [13] Xia Guodong, Liu Xianfei, Zhai Yuling et al, Single-phase and two-phase flows through helical rectangular channels in single screw expander prototype , J Hydrodyn Ser B 2014;26:114–121. 
       [14] Liu Xianfei, Xia Guodong, Zhai Yuling, Numerical analysis of the two-phase pressure drop and liquid distribution in single-screw expander prototype, Sci. bull 2014;59:4388–4396.
       [15] Xia Guodong, Liu Xianfei, An investigation of two-phase flow pressure drop in helical rectangular channel, Int Commun Heat Mass Transf 2014;54:33–41. 
       [16] Lu Xing, Du Xuepeng, Zeng Min et al Shell-side thermal-hydraulic performances of multilayer spiral-wound heat exchangers under different wall thermal boundary conditions, Appl Therm Eng 2014; 70: 1216–1227. 
       [17] Liu Xianfei, Wang Fang, Li Zhiqiang et al, Parametric investigation of thermal-hydrodynamic performance in the innovative helically coiled heat exchangers in the heat pump system, Energ Buildings, 2020; 216:109961
       [18] Shine SR, Kumar SS, Suresh BN, Numerical study of wave disturbance in liquid cooling film, Propuls Power Res 2013;2:107–118.
       [19] Gu Fang, Liu Chunjiang, Yuan Xiguang et al, CFD simulation of liquid film flow on inclined plates, Chem Eng Technol 2004; 27:1099–1104. 
       [20] Yu Jiawen, Jiang Yiqiang, Cai Weihua et al. Forced convective condensation flow and heat transfer characteristics of hydrocarbon mixtures refrigerant in helically coiled tubes, Int J Heat Mass Transf 2018;124:646–654.
       [21] Zhao Houjian, Li Xiaowei, Wu Xinin, Numerical investigation of supercritical water turbulent flow and heat transfer characteristics in vertical helical tubes, J Supercrit Fluids 2017;127:48–61.
       [22] Kang Zhongtao, Li Xiandong, Mao Xiaobing, Acta Astronautica Experimental investigation on the surface wave characteristics of conical liquid film, Acta Astronaut 2018;149:15––24.

       备注：本文收录于《建筑环境与能源》2020年10月刊总第37期（第22届全国暖通空调制冷学术年会文集）。版权归论文作者所有，任何形式转载请联系作者。