基于漂移流模型的管道一维气液两相流数值仿真与分析_阿衍学.pdf

机械工程师MECHANICAL ENGINEER2023 年第 8 期网址： 电邮：MECHANICAL ENGINEER基于漂移流模型的管道一维气液两相流数值仿真与分析阿衍学1，李建东2，岳少博2，周穆洁1（1.黑龙江省机械科学研究院，哈尔滨 150040；2.西安前沿动力软件开发有限责任公司，西安 710075）0引言在各种工程应用中都会遇到多相流，比如石油工业中油气藏、井筒和地面集输系统1；核动力装置热工水力系统2-3；油气分离装置等4-6。在本文中，我们针对一维气-液两相流模型，该模型描述了管道内的两相流，是由两个两相的质量守恒方程和一个混合相动量方程组成的一维漂移流模型，其中液相和气相的速度由漂移通量定律联系起来。本文针对所建立的气液两相流模型，使用稳定和精确的数值公式来有效地进行数值求解。在建立数学模型时，我们认为气相和液相之间没有相变，即所考虑的两种流体之间没有传质。从气-液两相流双流体模型的控制方程入手，并推导了系统的雅可比方程及其特征值，表明了它在工程应用的复杂性。在数值求解中，我们利用通量校正输运的思想，实现了通量的计算，其中用HLL型求解器，对两步Lax-Wendroff方法获得的通量进行平均，使得具有极限的非线性平均不违反Godunov定理。1数学模型1.1守恒形式的控制方程针对气液混合物的一维两相流问题，为了简化计算，假设液体和气体之间没有传质，并对管道和流体做如下假设：1）流动管道的截面积沿长度方向是可以变化的；2）两种流体的理想状态方程遵循等温流动条件；3）两种流体都是可压缩的（气液分别采用各自的状态方程）；4）两相的压力相等。基于以上假设，可推导出如下公式：m=（1）l+g；（1）um=（1）lul+gug（1）l+g；（2）uV=（1）ul+ug。（3）式中：为气相体积分数，无量纲量，取值范围为01；l为液相密度；g为气相密度；ul为液相速度；ug为气相速度；um为混合相的质量平均速度；uV为混合相的体积平均摘要：两相流模型受到双曲线性质的影响，因为很难获得雅可比矩阵的简洁解析表达式，从而限制了开发Roe型黎曼求解器。针对管道气液两相流流动现象，利用一维气液两相流漂移流模型，建立了管道一维气液两相流控制方程，采用基于时间的限制器，使用Lax-Wendroff方法求解方程组，并以此开发了数值求解器。以垂直管中油气两相流为模型算例，计算了不同时刻的管中沿管道方向的含气率分布，通过分析比较，文中采用的漂移流一维气液两相流模型和数值算法能够满足管道一维气液两相流流动的工程应用。关键词：一维两相流；Lax-Wendroff；漂移流模型中图分类号:TQ 022.4文献标志码:A文章编号：10022333（2023）08013204Numerical Simulation and Analysis of One-dimensional Gas-liquid Two-phase Flow in Pipelines Based onDrift Flow ModelA Yanxue1,LI Jiandong2,YUE Shaobo2,ZHOU Mujie1(1.HeilongjiangAcademyofMechanical Sciences,Harbin 150040,China;2.Xian Advanced Dynamics Incorporated,Xian 710075,China)Abstract:It is well known that the multi-phase flow model suffers from the clear definition of hyperbolic nature as it isdifficult to obtain a neat analytical expression for Jacobian matrix which limits us to formulate a Roe type Riemann solver.This paper focuses on the phenomenon of gas-liquid two-phase flow in pipelines,then establishes a one-dimensional gas-liquid two-phase flow control equation in pipelines based on drift flow model,and develops a numerical solver using Lax-Wendroff method and time-based limiters.Taking the oil-gas two-phase flow in a vertical pipe as an example,thedistribution of void fraction along the pipeline direction at different times is calculated.Through analysis and comparison,the drift flow one-dimensional gas-liquid two-phase flow model and numerical algorithm can meet the engineeringapplication of one-dimensional gas-liquid two-phase flow in pipelines.Keywords:one-dimensional two-phase flow;Lax-Wendroff;drift flow model132机械工程师MECHANICAL ENGINEER网址： 电邮：2023 年第 8 期MECHANICAL ENGINEER速度。根据气液两相的质量守恒，可分别得到两相的质量守恒方程：（1）lAt+（1）lulAx=0；（4）gAt+gugAx=M。（5）式中：M为管道外部渗入的气体流量，如果没有气体渗入，则为0；A为管道横截面积；x为沿管长方向的长度。结合式（1）、式（2）可得混合相的质量守恒方程为mAt+mumAx=M。（6）式（5）中的通量项的形式可通过变形，代入式（1）和式（2），改写为如下形式：gug=gug+gumgum=gum+gugum()=mum+gug（1）lul+gug（1）l+g()=mum+g mug-（1）lul-gugm()=mum+g（1）lugul()m=mum+g1gm()ugul()。（7）结合式（6）、式（7），得到质量平均速度的质量守恒方程：gAt+mum+m1gm()ugul()()Ax=M。（8）类似地，混合相的动量守恒方程，可写成下式：mumAt+（1）lu2l+gu2g+()p()Ax=hfA+pdAdxmgAcos。（9）式中：为管道与水平面的夹角；p为混合相的压力，等于各分相压力；hf为管道的沿程损失；g为重力加速度。同样，在式（9）中，对通量项中的各相密度用平均质量密度替换，可得到下式：mu2m（1）lu2l+gu2g()=（1）lul+gug()2m（1）lu2l+gu2g()=（1）2l2ul2+2g2ug2+2（1）lgulugm（1）lu2l+gu2g()m=2（1）lgulug（1）lgu2l+u2g()m=（1）lgulug()2m。（10）把式（10）代入式（9），可得mumAt+mu2mA+（1）lgulug()2mA+()Ax=hfA+pdAdxmgAcos。（11）1.2气液相的状态方程遵循1.1节中的4条假设，可有液相和气相的状态方程如下：pl=pref+a2llref()；（12）pg=ga2g。（13）式中：pref和ref分别为液相的参考压力与参考密度；al和ag分别为声音在液相和气相中的传播速度；pl和pg分别为液相和气相的压力，由1.1节中第4条假设可知，pl=pg=p。（14）1.3漂移流模型两相间速度差可由实验得到的公式计算，如下式：ugul=vs1C+mum（1C）gvs()（C1）mg()（1C）+gC（1）()（1C）。（15）其中，C和vs是关于和相速度的函数，相关表达式为ug=ulC（1）+vs1C。（16）将式（16）代入混合相方程可得：um=（1）lulm+gmulC（1）+vs1C()；（17）ul=mum（1C）gvs（1）l（1C）+gC（1）。（18）则式（15）可以变为ugul=vs1C+mum（1C）gvs()（C1）mg()（1C）+gC（1()）（1C）=vs1C+mum（1C）gvs()（C1）m（1C）g（1C）()（1C）。（19）通常，C和vs由实验确定的经验公式给出，如下式：C（）=1，00.015；1+4.1176（-0.015），0.0150.1；1.35，0.10.4；1+0.5833（1），0.41.0|。（20）vs（）=0，00.015；5.8823（-0.015），0.0150.1；0.5，0.10.4；0.8333（1），0.41.0|。（21）1.4体积分数方程根据两相的平均密度方程和状态方程，我们可以推出体积分数与密度、压力之间的关系表达式：（1）l=mg；（22）133机械工程师MECHANICAL ENGINEER2023 年第 8 期网址： 电邮：MECHANICAL ENGINEER（下转第138页）图1两相流瞬态计算流程图程序启动对管道划分一维网格读取介质参数、初始条件设置或更新边界条件空间离散（HLL）线性方程组时间离散（两步Lax-Wendroff 格式）迭代求解存储数据时间结束否是输出结果l=pprefa2l+ref。（23）根据式（22）和式（23）可以得到体积分数的表达式为=1mg()a2lppref+refa2l。（24）1.5两相流求解计算流程在两相流瞬态计算中，气液相和混合相的求解参数，如体积分数、混合相密度和压力，都是时间的函数，在编写求解程序时，既要在空间上对方程离散，也要在时间上对方程组离散，在本文程序中，对空间离散采用HLL格式，对时间离散采用两步Lax-Wendroff方法，程序流程图如图1所示，按照流程图编制一维两相流求解程序。2仿真结果分析为了验证本文中建立的漂移流模型的准确性、数值方法精度和程序的稳定性，以经典的气液两相流管道流动算例作为验证算例，通过本文的计算结果的分析，来验证本文所采用的方法的准确性和精度。本算例模拟气液两相在长100 m的垂直管中的流动，液相介质为石油、气相介质甲烷，对管道划分100个单元，CFL条件数设置为0.5，管道两端边界条件采用齐次纽曼边界条件，迭代计算121步。图2图5分别为初始时刻和计算终了时刻天然气体积分数、混合相密度、混合相质量平均密度和混合相压力沿管道长度方向的分布曲线图。在计算终了时刻，在管道中间呈现梯度较大的分布，在气液界面处没有出现振荡现象，同时界面处也没有被拉平得到过渡段。计算结果表明本文所采用离散格式具有较高的精度和稳定性。在x方向6070m界面处理论上应是垂直过渡，这与本文所采用的网格数量有关，可以通过增加网格数量来提高计算精度。在管道系统中，通常是由大量的单根支管所组成，单根支管网格数量如果过大，那么对于整个管道系统来讲，会极大地增大计算量，消耗大量计算时间。在满足工程应用 的 情 况下，划分适当的网格来保证精度。3结语与双流t0 st60 s0.500.300.350.400.450.25010020406080 x/mt0 st60 sm/（kg m-3）620460010020406080 x/m480500520540560580600440t0 st60 sum/（m s-1）30-5010020406080 x/m0510152025-10t0 st60 sP/kPa1150010020406080 x/m650700750800850900950600105011001000图2初始时刻和终了时刻气相体积分数沿管长的分布图3初始时刻和终了时刻混合相密度沿管长的分布图4初始时刻和终了时刻混合相质量平均速度沿管长的分布图5初始时刻和终了时刻混合相压力沿管长的分布134机械工程师MECHANICAL ENGINEER2023 年第 8 期网址： 电邮：MECHANICAL ENGINEER（上接第134页）体两相流模型相比，漂移流两相流模型的方程形式相对简单，这样对方程处理起来相对方便，便于编程实现。本文采用的HLL和两步Lax-Wendroff方法，能够在较少网格数量的情况下，获得较为满意的计算精度，使得本文的方法可以很好地应用于工程实践中，特别是对于复杂的两相流管网仿真计算。同时本文的计算方法还可以推广用于气-水-油三相流动。参考文献1徐朝阳.井筒多相流瞬态流动数值算法及响应特征研究D.成都：西南石油大学,2015.2MARTINELLI R C,PUTNAM J A,LOCKHART R W.Two-phase,Two-component Flow in the Viscous Region J.Transactions of American Institute of Chemical Engineers,1946,42:681-705.3LOCKHART R W,MARTINELLI R C.Proposed Correlation ofData for Isothermal Two-phase,Two-component Flow in PipesJ.Chemical Engineering Progress,1949,45(1):39-48.4陈之航,曹柏林,赵在三.气液双相流动和传热M.北京：机械工业出版社,1992.5劳力云,张宏建,陈洪,等.水平管道气液两相流动态特性模型的研究J.高校化学工程学报,1999,13(6):449-553.6阎昌琪.气液两相流M.哈尔滨：哈尔滨工程大学出版社,2010.（编辑马忠臣）作者简介：阿衍学（1979），男，本科，高级工程师，主要从事机械设计研发及项目管理工作。收稿日期：2023-03-065.09 kg，试验数据如表4所示。对比抛光后的检测数据，抛光量、圆度、平行度、直线度均有变大。相比于15.22 kg配重抛光后的数据，更加接近极限值，可见配重减小是不可行的。考虑加大配重质量，使抛光部分的压紧力减轻，试验抛光情况。配重在15.09 kg的基础上加重0.5 kg。实际配重的质量为15.58 kg。试验数据如表5所示。配重15.58 kg的抛光后的数据和配重15.09 kg的抛光数据相似，也在极限状态。可见配重加大也不可行。同时，在实际抛光过程中，由于配重过大，导致抛光部分质量轻。曲轴旋转后，抛光部分出现了上下跳动的情况，影响抛光质量。从数据对比看，配重质量为15.22 kg更加符合抛光设计要求。根据抛光量和圆度、平行度的变化情况，拟将配重设定在15.315.35 kg，因制造误差，实际配重为15.32kg，试验数据如表6所示。实际抛光过程中，抛光部分没有出现上下跳动的情况，工作状态稳定，且抛光后的数据符合抛光设计要求，也均在中下差，轴颈无振纹、纹理均匀，符合设计要求。通过计算，压在轴颈表面质量在0.6950.700 kg之间为最佳。4结语将UG软件的建模分析功能与实物验证相结合，可减少实物验证数量和周期，也充分证明模拟分析方法的准确趋向性，经过推算该型抛光装置的抛光部分与配置部分的质量比关系，该型抛光装置配重质量在15.32 kg时为最佳工作状态。该分析方法加速了抛光装置的实物应用，也为其他类型曲轴抛光装置的配重计算提供了参考。参考文献1范淑华，项林川.大学物理（上册）M.武汉：华中科技大学出版社,2019:59-68.2机械设计手册 编委会.机械设计手册（上册）M.北京：化学工业出版社,1987:656-694.3何耿煌，李凌祥，程程.UG NX 10.0从入门到精通M.北京：中国铁道出版社,2016:263-275.4张玉，刘平.几何量公差与测量技术M.沈阳：东北大学出版社,2014:153-167.5王德泉.砂轮特性与磨削加工M.北京：中国标准出版社,2001:98-106.（编辑马忠臣）作者介绍：姚旭（1984），男，本科，工程师，主要从事高速柴油机曲轴制造技术研究工作。收稿日期：2021-12-15表4配重15.09 kg的抛光数据抛光量/mm 圆度/mm 平行度/mm 粗糙度Rz/m 直线度/mm表面质量评价0.0150.010.011.00.005 7表面无振纹、光泽均匀表5配重15.58 kg的抛光数据抛光量/mm 圆度/mm 平行度/mm 粗糙度Rz/m 直线度/mm表面质量评价0.0140.010.01Rz1.350.007 2表面无振纹、光泽均匀表6配重15.32 kg的抛光数据抛光量/mm圆度/mm 平行度/mm粗糙度Rz/m 直线度/mm表面质量评价0.0100.0050.005Rz1.280.003 1表面无振纹、光泽均匀138