¹北京工业大学环境与生命学部，北京

²传热强化与过程节能教育部重点实验室，北京

³传热与能源利用北京市重点实验室，北京

收稿日期：2021年3月1日；录用日期：2021年4月12日；发布日期：2021年4月19日

摘要

本文通过采用CFD数值模拟分析法研究单螺杆膨胀机基元容积膨胀过程油气两相流的非稳态流动特性。考虑到单螺杆空间三维曲面结构的复杂性以及要实现单螺杆膨胀机基元容积膨胀过程的实时变化，利用STAR-CCM+软件对三维计算模型进行网格划分并采用变形网格技术实现动态膨胀过程的模拟仿真。模拟结果得到单螺杆膨胀机膨胀过程工质压力和温度变化规律，获得基元容积内流场的压力、温度和油气两相分布云图，并深入分析润滑油对工质状态的作用规律。

关键词

单螺杆膨胀机，膨胀过程，非稳态，变形网格技术，STAR-CCM+

Simulation Study on the Two-Phase Flow in the Expansion Process of Basic Element Volume of Single Screw Expander

Lili Shen^1,2,3, Yuting Wu^1,2,3, Wei Wang^1,2,3, Biao Lei^1,2,3, Ruiping Zhi^1,2,3

¹Department of Environmental and Life, Beijing University of Technology, Beijing

²Beijing Key Laboratory of Heat Transfer and Energy Conversion, Beijing

³MOE Key Laboratory of Enhanced Heat Transfer and Energy Conservation, Beijing

Received: Mar. 1^st, 2021; accepted: Apr. 12^th, 2021; published: Apr. 19^th, 2021

ABSTRACT

CFD simulation method has been used to investigate internal unsteady oil-gas two-phase flow characteristics of basic element volume of single screw expander during the expansion process. Due to the difficulty of three-dimension spacial curved surface of single screw structure and the real-time change of basic element volume of single screw expander during expansion process, the three-dimension computation model is meshed by using the STAR-CCM+ software and its morphing technology has been applied to simulate the dynamic two-phase flow of expansion process. The simulation results show that variation in the pressure and temperature of working fluid during the expander expansion process has been given, the cloudy chart of pressure, temperature and oil-gas two-phase distribution in the basic element volume has been obtained, and the effect of lubrication oil on the state of working fluid has been analyzed.

Keywords:Single Screw Expander, Expansion Process, Unsteady, Morphing Grid Technology, STAR-CCM+

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

计算流体力学(computational fluid dynamics, CFD)是一门通过计算机数值方法求解流体控制方程组来预测流体流动、传热、化学反应等相关物理现象的学科。CFD分析方法与实验方法相互补充共同解决实际问题，具有能够可视化、获取重要的物理量参数、减少实验次数、降低实验成本等优点，已广泛应用于航空航天、船舶、汽车、能源、化工、生物医学等领域。其中，在工业领域普遍使用的容积型流体机械压缩机和膨胀机也通常采用CFD方法进行分析研究，如活塞式 [1] [2]、滚动转子式 [3] [4]、涡旋式 [5] [6] [7]、滑片式 [8] [9] [10]、双螺杆式 [11] [12] [13] [14] 等等。然而，针对单螺杆结构的压缩机和膨胀机的CFD模拟仿真研究，目前还处在探索阶段，近几年才有公开发表的文献。

SUMAN等 [15] 和ZIVIANI等 [16] 基于四种流体力学仿真软件：ANSYS FLUENT，STAR-CCM+，CFX和PUMPLINX及其相对应的四种模拟方法(动网格技术、嵌套技术、关键帧重构方法和网格变形方法)，以单相空气作为工质，采用k-ε湍流模型分析比较单螺杆膨胀机稳态和瞬态的特性。研究结果得出，四种方法均可以进行稳态模拟，只有STAR-CCM+软件的嵌套网格技术可以对单螺杆膨胀机的瞬态特性进行模拟。CASARI等 [17] 针对基于有机朗肯循环系统的单螺杆膨胀机开展了CFD瞬态模拟研究，分析单螺杆膨胀机的流动特性和机械性能。然而，以上模拟过程中均未考虑润滑油的影响，只是模拟单相气体流动。因此，需对单螺杆膨胀机内油气两相流动开展非稳态模拟研究。文献 [18] [19] [20] 基于单螺杆机械螺杆槽道的几何结构，开展两相流流型、压降及液膜厚度等方面的模拟和实验研究，但这些研究没有与单螺杆膨胀机动态工作过程进行耦合研究。

本文通过采用STAR-CCM+软件中的变形网格技术，对单螺杆膨胀机基元容积膨胀过程非稳态两相流动进行CFD模拟研究。获得膨胀过程中工质的压力和温度变化规律、润滑油的分布特征及润滑油对工质状态参数的作用规律等，为今后单螺杆膨胀机整机两相流动数值模拟提供基础指导。

2. 非稳态模型的建立

2.1. 三维空间几何模型

单螺杆膨胀机的主要由螺杆转子、星轮、机壳、主轴以及轴承、轴封等零件构成。螺杆螺槽、星轮齿和机壳内壁面构成封闭的工作容积，该工作容积也称为基元容积，工质流体在基元容积内膨胀以此推动螺杆和星轮运动，向外输出轴功。图1给出了单螺杆膨胀机从进气结束基元容积开始膨胀到膨胀结束开始排气之前的整个膨胀过程。

图1. 单螺杆膨胀机膨胀过程示意图

考虑到单螺杆膨胀机整机三维空间结构的复杂性，本文以基元容积为研究对象，建立封闭膨胀过程的三维几何模型，通过采用CFD的动态变形网格技术对单螺杆膨胀机膨胀过程中非稳态油气两相流动进行数值模拟，进而直观地分析单螺杆膨胀机膨胀过程中润滑油的分布特性及其对工质气体的作用规律。图2从两个不同角度展示了基元容积的三维几何模型图，基元容积由顶面、螺槽左侧面、螺槽右侧面、螺槽底面和星轮齿面五个面构成的曲面结构。

图2. 三维几何模型

由于单螺杆膨胀机实际运行的复杂性，在建立单螺杆膨胀机非稳态模型时，做了一些简化和假设：

1) 单螺杆膨胀机内的工质为空气，满足理想气体状态方程；

2) 润滑油为不可压缩流体，材料属性按照壳牌S4R46型号参数设置；

3) 忽略泄漏和摩擦，认为基元容积内流体流动为闭口系绝热流动。

2.2. 网格划分

STAR-CCM+搭载了CD-adapco独创的最新网格生成技术，可以导入复杂形状数、进行表面网格重构、自动生成体网格等网格生成所需的一系列工作。本文采用STAR-CCM+软件的切割体网格划分变形流体区域。整体操作简单易行，并不需要太多的人工参与。只需要定义每个面上的网格基本大小和尺寸，就可以在面上自动生成网格，生成面网格后，通过选择体网格的类型和尺寸，便可自动光滑地生成体网格。考虑到螺槽通道的复杂曲面结构，通过面控制对局部强弯曲、细窄区域进行局部网格细化，使其更好捕捉流动形态。基元容积开始膨胀和膨胀终了时流体区域的网格，如图3所示。

图3. 三维几何模型的网格

表1给出了网格无性验证结果，具体是通过判断膨胀终了时基元容积内平均压力和气体温度与网格数的关系。可以看出，随着网格数的增加，基元容积内压力呈现下降趋势且变化逐渐趋于稳定，而基元容积内气体温度对网格敏感度非常小。当网格数达到10万之后，压力恒定在0.152 MPa左右，因此，为了节约求解时间和计算资源，同时保证计算结果的准确性，本文采用的网格数为10.15万。

表1. 标网格无关性验证结果

2.3. 欧拉多相流分离流模型

欧拉法着眼于空间的点，基本思想是考察空间一个点上的物理量及其变化。欧拉多相流分离流模型是对每一相都建立连续性方程、动量守恒方程，能量守恒方程，通过压力和相间交换系数的耦合来计算求解。空间各点都有多相流体各自不同的速度、温度和密度，各相间有滑移。这些流体存在同一空间并相互渗透，但各有不同的体积分数，各相的体积分数是时间和空间的连续函数，各相体积分数之和等于1。

${\displaystyle {\sum }_i{\alpha }_i=1}$ (1)

连续性方程为：

$\frac{\partial }{\partial t}{\displaystyle \underset{V}{\int }{\alpha }_i}{\rho }_i\text{d}V+{\displaystyle {\oint }_A{\alpha }_i{\rho }_i\left(v_i-v_g\right)}\cdot \text{d}a={\displaystyle \underset{V}{\int }{S}_i^{\alpha }\text{d}V}$ (2)

式中， ${\rho }_i$ 为i相流体的密度， ${\alpha }_i$ 为i相流体的体积分数， $v_i$ 为i相流体的的速度， $v_g$ 为网格速度， $S_i^{\alpha }$ 为i相流体的质量源项。

动量守恒方程为：

$\frac{\partial }{\partial t}{\displaystyle \underset{V}{\int }{\alpha }_i}{\rho }_i{v}_i\text{d}V+{\displaystyle {\oint }_A{\alpha }_i{\rho }_i{v}_i\times \left(v_i-v_g\right)}\cdot \text{d}a={\displaystyle \underset{V}{\int }{\alpha }_i}\nabla p\text{d}V+{\displaystyle \underset{V}{\int }{\alpha }_i}{\rho }_{i}g\text{d}V+{\displaystyle {\oint }_A{\alpha }_i\left({\tau }_i-{\tau }_i{}^t\right)}\cdot \text{d}a+{\displaystyle \underset{V}{\int }{S}_i^v\text{d}V}+{\displaystyle \underset{V}{\int }{M}_i}\text{d}V$ (3)

式中， $p$ 为压力，且认为对于各相均相等， $g$ 为重力向量， ${\tau }_i$ 和 ${\tau }_i{}^t$ 分别分子应力和湍流应力， $S_i^v$ 为i相流体的动量源项， $M_i$ 是单位体积的相间动量转移量。相间动量传递表示各相相互施加的所有力之和，满足以下等式：

${\displaystyle {\sum }_i{M}_i=0}$ (4)

能量守恒方程为：

$\begin{array}{l}\frac{\partial }{\partial t}\left({\alpha }_i{\rho }_i{E}_i\right)+\nabla \cdot \left[{\alpha }_i{\rho }_i{H}_i\left(v_i-v_g\right)\right]+\nabla \cdot {\alpha }_i{v}_{g}p=\nabla \cdot \left({\alpha }_i{k}_{eff,i}\Delta T_i\right)+\nabla \cdot \left({{\rm T}}_i\cdot v_i\right)\\ +f_i\cdot v_i+S_{u,i}+{\displaystyle {\sum }_{j\ne i}{Q}_{ij}}+{\displaystyle {\sum }_{\left(ij\right)}{Q}_i^{\left(ij\right)}}\end{array}$ (5)

式中， $E_i$ 为总能量， $H_i$ 为总焓值， ${{\rm T}}_i$ 为粘性应力张量， $T_i$ 为温度， $f_i$ 为体积力， $Q_{ij}$ 为从j相到i相的相间换热率， $Q_i^{\left(ij\right)}$ 为从(ij)相界面到i相的换热率， $S_{u,i}$ 为i相流体的能量源项， $k_{eff,i}$ 为有效热导率，则有：

$k_{eff,i}=k_i+\frac{{\mu }_{t,i}{C}_{p,i}}{{\sigma }_{t,i}}$ (6)

式中， $k_i$ 为热导率， ${\mu }_{t,i}$ 为湍流粘度， $C_{p,i}$ 为比热， ${\sigma }_{t,i}$ 为湍流热扩散普朗特数。

上述控制方程组变量多且不封闭，针对单螺杆膨胀机内非稳态流场模拟也采用工程上广泛应用的湍流数值模拟方法——雷诺平均模型。根据对雷诺应力做出的假设和处理方式的不同，湍流模型又可分为雷诺应力模型和湍动粘度模型。其中，湍动粘度模型中标准k-ε模型可解决实际绝大部分工程问题，但用于强旋流、绕弯曲壁面流动时，会产生一定的失真。因此，本文采用改进后的Realizable k-ε模型，并采用分离求解法进行迭代计算。

2.4. 初始和边界条件

本文模拟单螺杆膨胀机膨胀过程内油气两相的非稳态流动，初始条件按照进气状态参数设置，其中初始绝对压力为0.7 MPa，进气温度为293 K，润滑油初始温度为293 K，工质气体的体积分数为99.5%，润滑油的体积分数为0.5%。

计算域上壁面边界均为动态移动壁面边界条件，其中，顶面、螺槽左侧面、螺槽右侧面和螺槽底面按照螺杆转速设置，星轮齿面按照星轮的转速设置。具体是根据传动比11/6设置螺杆和星轮转子的转速，即：螺杆转子转速为3000 rpm，则星轮转速为1636.3636 rpm。

2.5. 求解流程

根据前面已设置好的物理模型、边界条件和初始条件等，便可进行数值求解计算。本文通过利用STAR-CCM+软件中的宏控制，通过采用JAVA语句编写程序，实现程序调用和语句操作。具体求解思路如流程图4所示，先根据内容积比设置，确定初始膨胀时间作为开始计算时间，获取该时刻对应基元容积的三维几何模型，之后便可对计算域进行网格划分求解一个时间步长，此时基元容积的动态边界会根据螺杆和星轮转速设置的大小进行啮合旋转运动，基元容积开始进行膨胀，随着膨胀容积的增加，需要实时替换基元容积三维几何模型，直到膨胀过程结束，输出计算结果。

图4. 数值计算流程图

3. 模型验证

为了验证模型的可靠性，针对117 mm标准CP型单螺杆膨胀机，基本参数如表2所示，通过将单相空气绝热膨胀过程的压力、温度变化的模拟计算值与理论计算值进行对比，如图5和图6所示。其中，虚线为理论计算值，实线为模拟计算值。从图中可以看出，在四种不同内容积比下(图中用τ表示)，模拟计算膨胀过程的压力和温度与理论计算值曲线变化基本重合，说明二者吻合较好。因此，可通过采用该模型对单螺杆膨胀机膨胀过程的油气两相流动进行模拟研究。

表2. 单螺杆膨胀机基本结构参数

图5. 压力模拟值与理论值的对比

图6. 温度模拟值与理论值的对比

4. 计算结果分析

本文以117 mm标准CP型单螺杆膨胀机为例，通过模拟喷油单螺杆膨胀机在螺杆转速3000 rpm，初始膨胀压力为0.7 MPa，温度为293 K，通过对比油气两相与单相流体的状态参数，分析润滑油对工质气体的作用规律，并给出膨胀过程中基元容积机内气体压力、温度及不同油气两相分布的情况。

4.1. 润滑油的影响

如图7所示，不同内容积比下的单螺杆膨胀机膨胀过程中基元容积内单相流和两相流的工质压力变化曲线。其中，两相流工质的压力采用实线表示，单相流工质的压力采用虚线表示，可以看出，在四种内容积比下，单相流和两相流工质的压力随螺杆转角的曲线几乎重合，说明考虑润滑油后对工质的压力影响非常小。

图8给出了不同内容积比下的膨胀过程中基元容积内单相流和两相流的工质温度变化关系，图中两相流工质的温度采用实线表示，单相流工质的温度采用虚线表示。可以看出，在相同的内容积比下，两相流中工质温度随着膨胀过程的进行会逐渐高于单相流工质的温度，主要是因为润滑油的比热容高于工质气体，随着膨胀过程的进行，气体在膨胀降温过程中会从润滑油中吸收一部分热量。结合图9两相流中润滑油和气体温度随着螺杆转角的变化曲线可知，润滑油在膨胀之初温度开始缓慢降低，且温度总高于同期工质气体温度，因此，两相流中油气相互作用使得工质温度会高于单相流的情况。此外，随着内容积比的增加，二者差值也逐渐增大，如内容积比为5.4的膨胀终了时单相流和两相流气体温差达到23 K，而内容积比为1.8的膨胀终了时二者气体温差不到7 K。这是由于内容积比越大，进气孔口越小，进气过程缩短，膨胀过程延长，润滑油对工质气体温度作用时间也越长，使得内容积比越大，单相流和两相流气体温度相差越大。

图7. 不同内容积比下单相流与两相流工质压力随螺杆转角的变化关系

图8. 不同内容积比下单相流与两相流工质温度随螺杆转角的变化关系

图9. 不同内容积比下工质和润滑油温度随螺杆转角的变化关系

4.2. 压力变化

图10给出了内容积比为3.0时，基元容积在膨胀过程中不同螺杆转角下的两相流的压力云图。从螺杆转角为65˚开始膨胀直至150˚膨胀结束，压力从0.7 MPa降低到0.152 MPa。结合图7中压力随螺杆转角的变化曲线，可以看出，膨胀过程中。压力随着容积增大呈现先快速下降后缓慢下降的趋势，并且从压力云图中可以看出在给定时刻(螺杆转角)基元容积内压力整体分布均匀。

图10. 基元容积的压力云图

4.3. 温度变化

图11给出了内容积比为3.0时，不同螺杆转角对应基元容积内两相流的温度云图。从图中可以看出，给定时刻(螺杆转角)基元容积内温度等温线分布较多，说明局部区域温度分布不均匀。结合图8中工质温度随螺杆转角的变化曲线可以看出，基元容积内的温度随着膨胀过程的进行呈现先快速下降后缓慢下降的趋势。从螺杆转角为65˚开始膨胀直至150˚膨胀结束，温度从293 K降低到202 K。

4.4. 油气体积分数的影响

通过定量表示基元容积内油和气体的体积分数，可以直观地看出二者的分布情况。其中，体积分数越接近1，表明其体积分数越高，而体积分数越接近0，表明该项占比越小。文献 [21] 指出喷油压缩机内一般油气体积比在1%左右，因此，本文通过选用油气体积比为0.5%和1%两种情况进行数值模拟。

图11. 基元容积的温度云图

膨胀过程中基元容积内油气两相分布云图，如图12(a)和图12(b)所示。从图中可以看出，膨胀初期油气两相均匀分布，膨胀过程中由于螺杆和星轮啮合旋转导致的离心力使得润滑油向顶面、螺槽左侧面及星轮齿面聚集。此外，随着油气体积比的增加，相同螺杆转角下油相分布面积增多。

图12. 不同油气比下油气两相分布云图

为了进一步分析油气体积比对工质压力和温度的作用规律，图13给出了不同油气体积比下膨胀过程基元容积内工质压力和温度随着螺杆转角的变化曲线。从图中可以看出，不同油气体积比对工质压力影响非常小，而随着油气体积比增大，工质温度在膨胀过程中下降速度减缓，且在膨胀终了时温度会升高。说明喷入润滑油的量对于单螺杆膨胀机的工质气体的压力作用很小，对工质气体的温度影响较大，尤其是膨胀终了时气体的温度。

图13. 两种不同油气比下工质压力和温度随着螺杆转角的变化关系

5. 结论与展望

本文将通过采用STAR-CCM+软件中的变形网格技术实现单螺杆膨胀机膨胀过程中基元容积内油气两相的非稳态流动。主要结论如下：

1) 两相流中润滑油对工质压力影响甚微，对工质温度影响较大，且随着内容积比的增加，影响更加显著。

2) 基元容积两相流动的压力云图总体分布变化均匀，而温度云图由于润滑油的影响分布不均匀且局部存在温度差。

3) 基元容积润滑油在膨胀过程中会向顶面、螺槽左侧面以及星轮齿面聚集。

4) 油气体积比的增加会导致局部润滑油聚集面积增大，且会使膨胀过程中工质气体温度升高。

虽然本文研究了单螺杆膨胀机基元容积内油气两相非稳态流动，但由于做了一些较为理想化的假设，后期需要深入的工作还有很多，如考虑泄漏因素模拟单螺杆膨胀机整机实际运行工作过程中的两相非稳态流动。

基金项目

国家重点研发计划项目(2017YFB0903603)和国家自然科学基金项目(51706004)。

文章引用

沈丽丽,吴玉庭,王 伟,雷 标,智瑞平. 单螺杆膨胀机基元容积膨胀过程两相流模拟研究
Simulation Study on the Two-Phase Flow in the Expansion Process of Basic Element Volume of Single Screw Expander[J]. 机械工程与技术, 2021, 10(02): 188-199. https://doi.org/10.12677/MET.2021.102022

参考文献