水平管内气液两相流诱导振动的数值研究

马晓旭， 田茂诚， 张冠敏， 冷学礼

(山东大学 能源与动力工程学院，济南　250061)



水平管内气液两相流诱导振动的数值研究

马晓旭， 田茂诚， 张冠敏， 冷学礼

(山东大学 能源与动力工程学院，济南250061)

利用多场耦合有限元软件ADINA，结合任意拉格朗日-欧拉(ALE)动网格方法，对水平弯管内气液两相流诱导振动问题进行了数值模拟研究；重点考察了管结构振动响应特性及流型、体积含气率β和分相折算速度等对激振力的影响。计算结果表明：管结构动态响应表现为低频带内随机振动和周期性振动的混合，响应主频随分相折算速度的增加而逐渐增大。脉动激振力均方根值FRMS在环状流时最大，长泡状流时最小；β在60%～100%范围内且固定两相平均流速时，FRMS随β的增大而减小；FRMS随分相折算速度的增加而增大。最后通过建立数学模型总结出了水平弯管内气液两相流诱导振动的机理，并提出了针对性的减振措施。

气液两相流；诱导振动；流固耦合；脉动激振力；任意拉格朗日-欧拉法

流体诱导振动问题广泛存在于化工管道、换热器和核电站组件等工业应用中。这些振动会导致剧烈的噪声及设备元件的损坏，严重影响设备的寿命和安全运行[1]，相反，利用流体诱导振动也能实现强化传热的目的[2-3]。因此为更好的预防和利用振动能量，很有必要对流体诱导振动机理及其影响因素进行深入研究。

目前针对管外流的研究较多，例如唐友刚等[4-5]对海洋油气开发立管的涡激振动机理的系统研究，PETTIGREW等[6-7]综述了气液两相流横掠换热管束诱导振动的机理及其各影响因素等。而对于气液两相流且流经管内时，因其各相密度和速度差异，导致各处流动形态复杂多变，在遇到弯头、异径管、三通等元件时，流体诱导的脉动激振力同样会诱发极大振动。同时管内气液两相流又普遍存在于各种换热器如核电站蒸汽发生器和化工输流管道中，所以管内气液两相流诱导振动问题的研究具有重要的现实意义和工程应用价值，但是针对两相流内流激振的研究，尤其是利用计算流体动力学(CFD)的方法还比较少。

YIH等[8]较早对T形管内气液两相流脉动激振力进行研究，发现频谱图的主频要比管系统固有频率低很多，因此认为两相流不稳定流动引起的激发作用在实际工程中不会造成太大的危害，没有引起足够的重视。而RIVERIN等[9]对内径为20.6 mm的U形和T形管进行了实验研究，认为管内气液两相流诱导振动不能忽略，并用韦伯数对作用力均方根值进行了关联，关联式在较宽的几何结构和水力条件下取得较好的一致性。TAY等[10]实验得出液体物理特性中的表面张力和黏度二项对管内气液弹状流诱导激振力的影响很小。

随着计算机技术的发展，数值模拟作为一种重要方法广泛应用于各研究中。而由于气液两相流和管结构间流固耦合作用的复杂性，很少见到针对管内气液两相流诱导振动(以下简称GLTP-FIV)问题数值模拟的报道。ADINA软件以其稳定的非线性求解方法、强大的多物理场耦合计算能力而被誉为有限元软件中不可多得的精品。因此本文利用ADINA数值软件并基于任意拉格朗日-欧拉法，对水平管内GLTP-FIV激振特性进行三维流固耦合数值计算；研究了水平管内GLTP-FIV的动态响应和脉动激振力特性，并对脉动激振力的影响因素进行了分析；系统探讨了水平管内GLTP-FIV激振机理及减振措施，以期对气液两相流运行下设备管道的防振设计及安全操作起到指导作用。

1　理论方程

1.1流体域控制方程

本文所研究气液两相流诱导管结构振动问题，属于典型的流固耦合过程。而因流固耦合边界会发生变形，所以流体模型需采用基于任意拉格朗日欧拉Arbitrary-Lagrangian-Eulerian(ALE)坐标的网格运动算法，求解的基本未知量不仅包括压力、速度，还包括位移。

流体运动采用黏性不可压缩的Navier-Stokes方程来控制。在ALE坐标系下，流体域的连续性方程和动量方程描述如下：

(1)

动量方程：

(2)

σij=-pδij+2μeij

(3)

(4)

式中:ui、ωi分别为流体和网格运动的速度分量，ALE描述与Eulerian描述的区别在于将速度ui用相对速度ui-ωi替代。若ωi为零，即为纯Eulerian系；若ωi等于流体质点速度，即为纯Lagrangian系。ρ为密度常数，p为压力,μ为运动黏性系数，fi为体积力矢量;σij为应力张量分量，eij为速度应变张量分量。

1.2固体域控制方程

流体激振力作用下的管结构动力学基本方程为：

(5)

1.3流固耦合原理

应用在流固耦合界面的基本条件[11]如下：

位移协调方程：

(6)

力平衡方程：

(7)

(8)

式中：hd为结构节点的位移。流体均布力积分为集中力施加到结构节点上。

本文选取迭代法求解流固耦合，流体方程和结构方程按顺序相互迭代求解，然后各自将每一步得到的结果提供给另一部分使用，直到耦合系统的解达到收敛，迭代停止。一个时间步内的基本迭代过程为[12]：①首先形成流体域的离散方程各项，因为流固耦合交界面处的流体速度等于对应的固体速度，然后采用ALE算法求解流体控制方程得到流体域物理量，同时更新流体和固体网格边界及网格速度。②计算流固耦合界面处固体域节点上的流体载荷。③形成固体域的离散方程各项，并在流固耦合交界面处施加由第②步得到的流体载荷，以此求解固体控制方程得出固体域物理量，同时更新固体和流体网格边界及网格速度。④循环第①～③步，直至流体域、固体域的计算物理量收敛。

2　模型与方法

2.1物理模型与工况设置

物理几何模型如图1所示，为内径25.4 mm，壁厚2.3 mm的水平90°弯管。固定管两端AA′和DD′，空气和水的混合流体以恒定的平均速度从AA′端流入，然后从DD′端自由流出。AB段长350 mm，CD段长255 mm，弯头曲率半径OO′为50 mm，M点为脉动激振力及响应输出测点。

管材及流体物理特性如表1所示。

图1　几何模型及控制体积Fig.1 Geometric model and control volume

物理特性值ρ水/(kg·m-3)998.2ρ空气/(kg·m-3)1.2μ水/(kg·m-1·s-1)1.004×10-3μ空气/(kg·m-1·s-1)1.810×10-5K水/MPa2150K空气/MPa0.1416ρ管/(kg·m-3)1400E管/MPa3000泊松比v0.35σ/(N·m-1)0.0725

流体参数定义如下：

其中β为体积含气率；um为两相平均流速；usg和usl为分相(气相、液相)折算速度；A为总流通面积；Qg、Q1分别指空气和水的体积流量；D为管内径；ρL为液体密度；σ表示液体的表面张力系数；We为韦伯数。

将数值模拟工况绘制于BAKER[13]的水平管气液两相流流型判别图中，如图2所示。固定usl=0.3 m/s、0.5 m/s、0.8 m/s中某一液相折算速度时，气相折算速度选取从0.2 m/s到25 m/s中的10个工况。模拟工况覆盖长泡状流、弹状流和环状流等3种流型。

图2　模拟工况设置Fig.2 Simulation conditions

2.2数值方法及模型验证

固体模型在ADINA-Structures模块中建立，采用ADINA-native建模方式，管两端施加简支约束条件，弯管单元采用8节点3D solid单元，网格密度在线、面、体建模时依次设置。

流体模型在ADINA-CFD模块中建立，也采用native建模方式，流体单元采用的是8节点3D Fluid单元，并依据两相流的工况设置管内流体参数。

为获得本文数值解受网格疏密程度的影响，首先以液相折算速度usl为0.5 m/s、气相折算速度usg为1.5 m/s弹状流工况下管内气液两相流诱导振动问题为例，对数值模型的网格独立性进行验证，如表2所示。表中给出了激振力均方根值(FRMS)的对比结果。由表2可知，四种网格下的计算结果相差不大，网格C和网格D的数值计算结果偏差更是接近1%，表明当前网格密度对数值计算结果的影响很小，已达到网格独立性；综合考虑网格数量与计算精度，择优选取网格C的划分方法。此时固体与流体计算域的有限元单元模型如图3所示。其他流型均采用类似方法进行了网格独立性验证。

表2　usl=0.5 m/s且usg=1.5 m/s时网格独立性验证

图3　计算域有限元单元模型Fig.3 Finite element model of computational domain

流体的湍流模型选用标准k-ε模型，压力选用STANDARD格式，压力速度耦合方式采用压力隐式计算法(PISO)格式。数值求解时采用VOF模型对气液自由界面流动进行追踪，表面张力通过连续表面作用力(CSF)模型进行计算。

边界条件：入口面采用速度入口(inlet)边界，出口面采用自由流出口边界，与流体接触的管壁内侧设置为流固耦合(FSI)边界条件。

图4(a)给出了激振力系数CF-RMS随平均流速um的变化及其与CARGNELUTTI[14]中理论模型的比较结果。

(9)

FRMS为脉动激振力的均方根值，通过下式计算：

(10)

式中，Fi为激振力瞬时值，N是激振力瞬时值的数量。

CARGNELUTTI[14]的理论模型为：

(11)

式中：ρm为两相平均密度。

图4　数值计算结果验证Fig.4 The verification of numerical results

图4(b)对比了弹状流与环状流激振力主频fpeak的无量纲斯特劳哈尔数St的模拟与AZZOPARDI[15]中的实验结果。

St=fpeakD/[um(1-β)]

(12)

式中，fpeak为脉动激振力的主频，D为管内径、β为体积含气率。

从图4(a)～(b)中可以看出，本文数值计算结果与已有理论模型与实验结果吻合良好，表明本文计算模型具有较可靠的数值精度。

3　结果与分析

3.1动态响应特性

本文采用响应加速度的功率谱密度(以下简称PSD)来表征管道振动的强弱。图5(a)～(c)是液相折算速度为0.5 m/s时三种不同工况下，弯头y方向加速度PSD图。从图中可以看出，两相流管道振动加速度谱的成分复杂，振动谱线多分布在0～30 Hz，属于低频振动，峰值频率多集中在0～5 Hz范围内。图5(a)～(c)依次为长泡状流、弹状流和环状流。长泡状流下的PSD图没有明显的尖峰，动态响应表现为较窄频带下的随机振动；当流体流动转至弹状流和环状流时，功率谱图有一个较高的尖峰，其它频率处也存在一些相对较低的峰值，此时周期性振动成分占据主要作用。因此水平管内两相流诱导振动响应表现为窄频带内的随机振动和周期性振动的混合。

图6描述了弯头x方向振动加速度峰值频率(以下简称主频)随分相折算速度的变化规律。由图可见，振动响应主频低于6 Hz，处于低频振动范围，相比其它流型长泡状流的最低。固定液相流速，响应主频随气相折算速度的增加而逐渐增大；相同气相流速下，响应主频随液相流速的提高也增大；说明分相折算速度对响应主频影响较大。

图5　y方向加速度PSD图Fig.5 PSD of y-direction acceleration

图6　x方向加速度主频随分相折算速度的变化Fig.6 Predominant frequencies of x-acceleration versus phase superficial velocity

3.2脉动激振力

气液两相流诱发的激振力是管系统振动的外部激励，因此下面着重分析激振力的变化规律及影响因素。

3.2.1流型的影响作用

从图7中可明显看到，长泡状流时的激振力较小，弹状流次之，环状流时激振力比长泡状流大一个数量级。图9中分相折算速度的影响分析进一步证实了流型的影响规律。分析原因可知，韦伯数较小时水平管内两相流流速较小，多为长泡状流型，气液界面间的表面张力较大而惯性力较小，从而导致两相流体冲击力较弱且动量变化较小。

而由于弹状流型下液塞对管道截面的阻塞作用，气相的动能不能随气相流速的增大而成比例的增加。此能量一部分积聚在气弹内转化为内能，气弹内部压力升高；另一部分通过气弹前部与液塞尾部的动量交换转化成液塞的动能，带动液塞以接近气相的速度加速前进。当气弹内的能量累积到使其惯性力大于液相重力时，液相在靠近管顶的波峰处被击穿，气相在管中央成为连续相；而管顶处残留的液相在管壁的附着力与表面张力共同作用下贴附于壁面四周形成环状流液膜；流动即转变为环状流。所以，环状流时的流速更大，对壁面的冲击力更强；其次由于液相密度更大，相近流速下环绕管壁的液膜层比气弹传输更多的动量值，同时由于环状流气液界面的扰动特性，单位时间内弯管处液膜层的动量变化也远大于气弹的变化。

图7　流型对FRMS的影响Fig.7 Effect of flow regime on FRMS

3.2.2体积含气率β的影响

图8给出了有限的弹状流工况下，脉动激振力FRMS与体积含气率β间的关系。由图可知，固定气液两相平均流速且体积含气率β在60%～100%范围内时，脉动激振力随体积含气率的增加而降低，且平均流速越大时激振力降低速率越快。同一体积含气率β下，脉动激振力随平均流速的增加而增大。分析其原因，主要是因为当体积含气率增大时，单位时间流过弯头截面的气相份额增加，相应的液相份额减小，从而降低了两相流传输动量值，也减弱了液弹对管壁的冲击作用。换个角度说明弯头液相比重的增大，增加了两相流动量波动的程度及对管壁的碰撞冲击效应。而当增加气液平均流速时，两相流质量流量相应增大，对管壁结构的脉动冲击变大；流速的增大也会加大管内流体的湍流扰动强度和能量，相应的增强了作用在管壁结构上的激振力。

图8　FRMS值随β的变化规律Fig.8 FRMS versus volume void fraction

3.2.3分相折算速度us的影响

分相折算速度us对脉动激振力的影响见图9。图中可见，相同液相折算速度usl下，脉动激振力随气相折算速度的增加而增大，且增大的速率越来越快。固定气相折算速度usg，液相折算速度越大的脉动激振力也越大。长泡状流范围内激振力很小且波动较小，环状流下的激振力最大，弹状流次之。因为BLEVINS[16]指出流体激振力与流体密度和流速间通常存在F∝ρu2关系，由此式可见，激振力很大程度上与流速平方成正比。同时也是由环状流的流动特性所决定，因为环状流下围绕管壁的是一层环形液膜，而本研究中液膜密度约是气相的1 000 倍，所以相同流速下液膜传输动量更大，其对壁面的碰撞冲击效应也越明显。

图9　FRMS随分相折算速度us的变化规律Fig.9 FRMS versus phase superficial velocity

3.3诱导振动原因及机理

以上讨论已经得出振动响应与脉动激振力的变化规律及其影响因素，下面结合数学公式的推导来总结水平管内气液两相流诱导振动的物理机理。

取包含管壁与轴向两个有效截面构成的流管作为研究控制体，如图1所示。

流体系统内动量总量的输运公式为：

(13)

控制体内流体所受外力之和为：

(14)

考虑流体与壁面间的流固耦合作用，因此控制体表面所受力由两部分组成，一是管结构壁面对流体的作用力Fpipe，二是控制体进出口端面压力，所以

(15)

由动量定理可知，流体系统动量总量的时间变化率等于作用在系统上的外力之和，即

(16)

将式(13)～(15)代入式(16)并整理得出：

(17)

由牛顿第三定律得出，流体对管壁的激振力为

(18)

对于气液两相流，考虑均相流假设时，式(18)转换为：

(19)

式中:α为截面含气率。

分析式(19)各组成部分，管内气液两相流诱导激振力主要由四部分组成。

(1) 是控制体动量变化产生的作用力；这也正与YIH等[8]中激振力与流体动量的时间变化率有关的结论相一致。

(2) 是控制面动量变化产生的作用力；对于本研究控制体管壁表面没有流体通过，故仅包含单位时间内控制体端面流体动量的净通量。

上述两项主要由气液两相流不稳定性如速度、截面含气率的波动和流向改变所决定。而当两相流在流动转向时，上述两项引起的作用力主要包括流体离心力和流体冲击力两部分。

(3) 是重力项；即控制体内流体的质量力变化，主要由气液两相密度、速度的差异性和波动性引起。

(4) 是压力项；即进出口端面压差随时间的波动。

当管结构的弹性、阻尼、固有频率等参数固定不变时，管内气液两相流诱发的脉动激振力成为管结构振动的重要因素。因此，通过上述数学模型的分析并结合计算结果的讨论，总结水平弯管内气液两相流诱导振动的原因及机理如下：① 两相流动量波动激振；② 转向时液弹或液膜的碰撞冲击效应；③ 进出口压差波动激振。

3.4气液两相流管道减振措施

根据上述结果的分析和讨论，可采取以下几个措施以达到气液两相流管道减振的目的：

(1) 改变结构参数

上述分析得出振动响应是低频振动(0～30 Hz)，其气液两相流激励频率也将属于低频范围，所以应尽量增大管结构固有频率，以避免共振现象的发生。而管结构固有频率的增加可通过减小管道的跨度、提高支承刚度和增加管道壁厚等来实现。

(2) 调整操作参数

改变两相流速度、温度、压力等操作参数以避开激振力较大的弹状流和环状流，同时也减弱了液弹或液膜对管道转向处的碰撞冲击作用。

(3) 优化管路布局

为降低流体转向对管壁的冲击力，在满足静力分析要求的前提下应尽可能减少转弯的数量，并尽量减小转弯角度等。

4　结　论

本文利用有限元ADINA数值方法对水平弯管内气液两相流诱导振动问题开展了三维非稳态数值研究。在验证数值模型可靠的基础上，重点研究了诱导振动的响应特性及脉动激振力的影响因素。最后探讨了诱导振动的机理并提出了针对性的减振措施，主要结论如下：

(1) 利用ADINA流固耦合模型和VOF模型可以实现管内气液两相流诱导振动的模拟，为进一步研究流体诱导振动问题提供了新的方法和途径。

(2) 水平管内气液两相流诱导振动属于低频振动，动态响应表现为窄频带内随机振动和周期性振动的混合。响应主频随分相折算速度的增加而逐渐增大。

(3) 环状流下的激振力最大，长泡状流时最小。体积含气率β在60%～100%范围内且固定两相平均流速时，激振力均方根值FRMS随β的增加而减小。FRMS随分相折算速度的增加而增大。

(4) 水平弯管内气液两相流诱导振动的原因及机理为：两相流动量波动激振；转向时液弹或液膜的碰撞冲击效应；进出口压差波动激振。

(5) 气液两相流管道系统的减振可采取增加管道固有频率、避免液膜或液弹发生、减少流路转弯数量等措施来实现。

本文研究仍然存在不完善之处，研究结果也只是初步探讨，更为全面的了解水平管内气液两相流诱导振动的机理及其它影响因素还有待更深入的研究。

[1] 程林. 换热器内流体诱导振动[M]. 北京：科学出版社,2001.

[2] 程林, 田茂诚, 张冠敏，等. 流体诱导振动复合强化传热的理论分析[J]. 工程热物理学报, 2003, 22(3): 330-332.

CHENG Lin,TIAN Maocheng,ZHANG Guanmin,et al. Theo-retical analysis of complex heat transfer enhancement by flow induced vibration[J]. Journal of Engineering Thermophyscs,2003, 22(3): 330-332.

[3] 田茂诚, 林颐清, 程林，等. 汽-水换热器内流体诱导振动强化传热试验[J]. 化工学报, 2001, 52(3): 257-260.

TIAN Maocheng, LIN Yiqing, CHENG Lin,et al. Experimen-tal investigation of heat transfer enhancement by flow induced vibration in steam-water heat exchangers [J]. Journal of Chemical Industry and Engineering,2001,52(3): 257-260.

[4] 唐友刚, 樊娟娟, 张杰,等. 高雷诺数下圆柱顺流向和横向涡激振动分析[J]. 振动与冲击, 2013, 32(13): 88-92.

TANG Yougang, FAN Juanjuan, ZHANG Jie, et al. In line and transverse vortex-induced vibration analysis for a circular cylinder under high Reynolds number[J]. Journal of Vibration and Shock, 2013, 32(13): 88-92.

[5] 黄维平, 刘娟, 唐世振. 考虑流固耦合的大柔性圆柱体涡激振动非线性时域模型[J]. 振动与冲击,2012, 31(9): 140-143.

HUANG Weiping,LIU Juan,TANG Shizhen. Nonlinear model of vortex induced vibration of flexible cylinder in consideration of fluid-structure interaction[J]. Journal of Vibration and Shock, 2012, 31(9): 140-143.

[6] PETTIGREW M J,TAYLOR C E.Two-Phase flow-induced vibration:an overview[J].ASME Journal of Pressure Vessel Technology, 1994, 116(3):233-253.

[7] PETTIGREW M J, TAYLOR C E. Vibration analysis of shell-and-tube heat exchangers: an overview. Part 2: vibration response, fretting wear,guidelines[J].Journal of Fluids and Strutures,2003,18:485-500.

[8] YIH T S, GRIFFITH P. Unsteady momentum flux in two-phase flow and the vibration of nuclear system components[C]// Proceedings of the International Conference on Flow-Induced Vibrations in Reactor System Components, Argonne Ⅲ, 1970, Report ANL-7685: 91-111.

[9] RIVERIN J L, DE LANGRE E, PETTIGREW M J. Fluctuating forces caused by internal two-phase flow on bends and tees[J]. Journal of Sound and Vibration,2006,298:1088-1098.

[10] TAY B L, THORPE R B. Effects of liquid physical properties on the forces acting on a pipe bend in gas-liquid slug flow[J]. Chemical Engineering Research Design,2004,82:344-356.

[11] ZHANG H, BATHE K J. Direct and iterative computing of fluid flows fully coupled with structures[J].Computatiotional Fluid and Solid Mechanics,2001,1440-1443.

[12] ZHANG Q, HISADA T. Analysis of fluid-structure interaction problems with structural buckling and large domain changes by ALE finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2001,190: 6341-6357.

[13] GHIAASIAAN S M. Two-phase flow, boiling, and condensation[M].Cambridge:Cambridge University Press, 2008.

[14] CARGNELUTTI M F, BELFROID S P C, SCHIFERLI W. Two-phase flow-induced forces on bends in small scale tubes[J]. ASME Journal of Pressure Vessel Technology, 2010,132(4): 1-7.

[15] AZZOPARDI B J, BAKER G. Characteristics of periodic structures in gas/liquid two-phase flow [C]//Proceedings of UK/Japan Two Phase Flow Meeting.Guildford,UK,2003.

[16] BLEVINS R D. Flow-induced vibration[M]. Malabar,FL: Krieger Publishing Company,2006.

Numerical investigation on gas-liquid two-phase flow-induced vibration in a horizontal tube

MA Xiaoxu, TIAN Maocheng, ZHANG Guanmin, LENG Xueli

(School of Energy and Power Engineering, Shandong University, Jinan 250061, China)

Based on the arbitrary Lagrangian-Eulerian(ALE) dynamic mesh technique, gas-liquid two-phase flow induced vibration(GLTP-FIV) in horizontal elbows was investigated numerically with the multi-fields coupling FE software ADINA. The vibration response characters and the effect of flow regime, volume void fraction(VVF) and phase superficial velocity(PSV) on excitation forces were discussed. The results show that the dynamic response corresponds to a combination of low-frequency band random vibration and periodic vibration. The predominant frequency increases with the rise of PSV. The root mean square of excitation forces(F-RMS) increases gradually when the flow regime is in elongated bubbly, slug and annular flow. For a given two-phase mix averaged velocity, F-RMS decreases with VVF when it is in 60%-100%. And F-RMS also increases with PSV. Finally, the mechanisms of GLTP-FIV in elbows were summarized by establishing the mathematical model. And some vibration reduction measures were presented.

gas-liquid two-phase flow; flow-induced vibration(FIV); fluid-structure interaction(FSI); pulsating excitation forces; arbitrary lagrangian-eulerian(ALE) method

山东科技发展计划项目(2012GGX10421)

2015-03-25修改稿收到日期：2015-07-21

马晓旭 男，博士生，1984年11月生

田茂诚 男，博士，教授，博士生导师，1965年10月生E-mail:tianmc65@sdu.edu.cn

TB123

A

10.13465/j.cnki.jvs.2016.16.033