0 引　言

两相流与结构的流固耦合现象广泛存在，如LNG-FPSO船舶的复杂输流管道、海上石油开采平台的输送管道、动脉中的血液流动^[1]等。在实际工程中，为了更改管路走向，水平U型管的应用十分广泛。流体在管道中流动时会产生脉动激励力作用于管道内壁，引发管道结构振动，可能对安全生产和生态环境造成威胁，因此有必要对管道中流体流动形态及其诱导的管道振动响应进行深入研究。

多位学者对于管内流体流动形态及流激振动已经有了一定的研究。刘行等^[2]探讨了含冰率、流速、初始溶液浓度等因素对冰浆流动形态、流变特性及流动压降的影响。高岳等^[3]研究了立管振动响应特性以及振动对管内气液两相流动的影响。王志伟等^[4]分析了通过90°弯管后气液两相流的速度变化、截面含气率、压力分布及流型发展规律。

在水平管道中主要存在分层流、波状流、环状流、段塞流4种流型，一般而言，段塞流比其他流型更容易发生结构振动，但目前针对段塞流流体激励和振动响应数值仿真的相关研究较少。段塞流是一种不稳定的动态流动，通常表现为连续的无气液体或连续的无液气体^[5]。Faiza Saidj等^[6]研究了段塞流沿管道的发展和演变，Wang等^[7]发现段塞速度和长度是影响管道振动冲击的重要因素。

本文基于计算流体力学和计算结构动力学基本原理，应用Ansys有限元分析软件，对水平U型管中多种气液两相流流型进行对比分析，明确了段塞流工况是诱发U型管振动最剧烈的流型工况，进而研究了段塞流工况下的流体动态激励特性及管道振动响应特性。

1 数值仿真方法及相关原理

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

$ {\boldsymbol{M}}\ddot {\boldsymbol{u}} + {\boldsymbol{C}}\dot {\boldsymbol{u}} + {\boldsymbol{Ku}} = {\boldsymbol{F}}。$

(1)

式中： $ {\boldsymbol{M}} $ ， $ {\boldsymbol{C}} $ ， $ {\boldsymbol{K}} $ 分别为管道质量矩阵、阻尼矩阵和刚度矩阵； $ \ddot {\boldsymbol{u}} $ ， $ \dot {\boldsymbol{u}} $ ， $ {\boldsymbol{u}} $ 分别为管道加速度向量，速度向量和位移向量； $ {\boldsymbol{F}} $ 为管道所受的外界载荷。可利用显式求解中的中心差分法实现对式（1）的求解，管道振动的速度向量和加速度向量的迭代方程为：

$ \begin{split} & {{{\dot {\boldsymbol{u}}}_t} = \frac{1}{{\Delta {t^2}}}\left( {{{\boldsymbol{u}}_{t - \Delta t}} - {{\boldsymbol{u}}_t} + {u_{t + \Delta t}}} \right)}，\\ & {{{\ddot {\boldsymbol{u}}}_t} = \frac{1}{{2\Delta t}}\left( { - {{\boldsymbol{u}}_{t - \Delta t}} + {{\boldsymbol{u}}_{t + \Delta t}}} \right)} ，\\ & {{{\boldsymbol{u}}_{ - \Delta t}} = {{\boldsymbol{u}}_0} - \Delta t{{\dot {\boldsymbol{u}}}_0} + \frac{{\Delta {t^2}}}{2}{{\ddot {\boldsymbol{u}}}_0}} 。\end{split} $

(2)

式中： $ {{\boldsymbol{u}}_0} $ ， $ {\dot {\boldsymbol{u}}_0} $ ， $ {\ddot {\boldsymbol{u}}_0} $ 分别为初始时刻管道振动位移，速度，加速度； $ \Delta t $ 为差分求解的迭代时间步长； $ t $ 为指定时刻。由此，可计算得到指定时刻 $ t $ 的等效载荷：

$ \bar {\boldsymbol{F}} = {\boldsymbol{F}} - \left( {{\boldsymbol{K}} - \frac{2}{{\Delta {t^2}}}{\boldsymbol{M}}} \right){{\boldsymbol{u}}_t} - \left( {\frac{1}{{\Delta {t^2}}}{\boldsymbol{M}} - \frac{1}{{2\Delta t}}{\boldsymbol{C}}} \right){{\boldsymbol{u}}_{t - \Delta t}}，$

(3)

故 $ t + \Delta t $ 时刻的结构振动位移向量可通过以下方程进行计算：

$ \left( {\frac{1}{{\Delta {t^2}}}{\boldsymbol{M }}+ \frac{1}{{2\Delta t}}{\boldsymbol{C}}} \right){{\boldsymbol{u}}_{t + \Delta t}} = \bar {\boldsymbol{F}} 。$

(4)

流体域控制方程包含质量和动量守恒方程，如下式：

$ \nabla \cdot {\boldsymbol{V}} = 0，$

(5)

$ \rho \left( {\frac{{\partial {\boldsymbol{V}}}}{{\partial t}} + {\boldsymbol{V}}\nabla \cdot {\boldsymbol{V}}} \right) = - \nabla p + \rho {\boldsymbol{f}} + \mu {\nabla ^2}{\boldsymbol{V}} 。$

(6)

式中： $ {{{\boldsymbol{V}}}} $ ， $ {{{\boldsymbol{f}}}} $ 分别为速度矢量，体积力矢量； $ \rho $ ， $ p $ ， $ \mu $ 分别为流体密度，压力，动力粘度。

本文采用几何重构的VOF模型，其适用于求解具有较明显相界面的两相流，具有较高的精确度。在每个时间步内，根据流体域有限控制体内的空泡份额来捕捉相界面。湍流模型为标准 $ k - \varepsilon $ 模型， $ k $ 为湍动能， $ \varepsilon $ 为湍流耗散率。相间作用力模型为连续表面张力框架模型，同时设定入口、出口和壁面边界条件。引入VOF模型后，流体域的质量和动量控制方程可以表达为带相体积分数的方程：

$ \frac{\partial }{{\partial t}}\left( {{r_\alpha }\rho } \right) + \nabla \cdot \left( {{r_\alpha }\rho{\boldsymbol{ V}}} \right) = 0，$

(7)

$ \begin{split} \frac{\partial }{{\partial t}}\left( {\rho {\boldsymbol{V}}} \right) + \nabla \cdot \left( {\rho {\boldsymbol{V}} \otimes {\boldsymbol{V}}} \right) =& \nabla \cdot \left( {\mu \left( {\nabla {\boldsymbol{V}}} \right) + {{\left( {\nabla {\boldsymbol{V}}} \right)}^{\rm{T}}}} \right) -\\ & \nabla p + \rho f + \sigma k{\boldsymbol{n}}。\end{split} $

(8)

式中： $ {r_\alpha } $ ， $ \sigma $ ， $ k $ ， ${\boldsymbol{n}}$ 分别为相体积分数，相间表面张力系数，相界面局部曲率，相界面处的单位法向量。根据VOF方法基本原理，相体积分数满足以下关系：

$ \sum\limits_{\alpha = 1}^2 {{r_\alpha }} = 1，$

(9)

流固耦合界面的位移协调和应力平衡方程如下：

$ {d_{{f}}} = {d_{{s}}}，$

(10)

$ {n_{{f}}}{\tau _{{f}}} = {n_{{s}}}{\tau _{{s}}}。$

(11)

式中： $ {d_{\text{f}}} $ ， $ {d_{\text{s}}} $ 分别为流体位移和固体位移； $ {\tau _{\text{f}}} $ ， $ {\tau _{\text{s}}} $ 分别为流体应力和固体应力。

2 CFD数值模型及验证 2.1 模型及相关参数定义

本文研究对象为水平U型管，其在需要更改管路走向的工程管道系统中应用广泛。如图1所示，管道位于水平面XOY内，重力加速度沿Z轴负方向，其值为9.81 $ {\rm{m}}/{{\rm{s}}^2} $ 。 $ {L_{\text{1}}} = 5\;{\rm{m}} $ ， $ {L_{\text{2}}} = 6\;{\rm{m}} $ ，弯管半径 $ R = 1\;{\rm{m}} $ ，管道内径为240 mm。流体工质为水—空气混合流体，物理特性如表1所示，相关参数定义如下：

$ {J_{{w}}} = \frac{{{Q_{{w}}}}}{A} {J_{{v}}} = \frac{{{Q_{{v}}}}}{A} ，$

(12)

$ J = \frac{{{Q_{{w}}} + {Q_{{v}}}}}{A} = {J_{{w}}} + {J_{{v}}} ，$

(13)

$ \beta = \frac{{{Q_{{v}}}}}{{{Q_{{w}}} + {Q_{{v}}}}} = \frac{{{J_{{v}}}}}{J}。$

(14)

式中： $ {Q_{\text{w}}} $ ， $ {Q_{\text{v}}} $ 分别为水和空气的体积流量； $ A $ 为管道流通横截面积； $ {J_{\text{w}}} $ ， $ {J_{\text{v}}} $ 分别为液相、气相折算速度； $ J $ 为两相平均流速； $ \beta $ 为体积含气率。

2.2 网格无关性及数值仿真框架验证

在均质流假设下，两相流在某一截面的动量通量 $ M $ 可由空泡分数计算：

$ M = A{J^2}\left( {{\rho _{\text{v}}}{\alpha ^A}\left( t \right) + {\rho _{\text{w}}}\left( {1 - {\alpha ^A}\left( t \right)} \right)} \right) ，$

(15)

$ \Delta M = F\Delta t。$

(16)

式中：A为截面积；J为混合流体平均速度； $ {\rho _{\text{v}}} $ ， $ {\rho _{\text{w}}} $ 分别为气体，液体密度； $ {\alpha ^A}\left( t \right) $ 为t时刻该截面的空泡分数； $ F $ 为流体激励力。由上式可知，气液两相流在某一截面的动量通量同流体的空泡分数成正比，动量通量的变化同流体激励力成正比，故可通过评估空泡分数来分析流体激励力特性。以水平U型管段塞流工况（液速 ${J_{\text{w}}} = 1{\text{ }}{\rm{m}}/{\rm{s}}$ ，气速 ${J_{\text{v}}} = 2.44{\text{ }}{\rm{m}}/{\rm{s}}$ ）为例，取U型管的3个监测面（见图1）空泡分数时间序列均方根值的平均值，作为网格无关性验证的衡量指标，计算式如下：

$ \alpha _{RMS}^A = \sqrt {\frac{1}{3}{{\sum\limits_{i = 1}^3 {\left( {\alpha _i^A} \right)} }^2}} 。$

(17)

由表2可知，随着单元和节点数目的增多，相对误差越来越小，方案4与方案3的相对误差仅为0.1643%。对结果精度和计算代价进行权衡，选择方案3作为仿真网格，如图2所示。对其他工况下的管道气液两相流数值仿真，均进行了相应的网格无关性验证来确定网格密度。

为验证数值仿真框架的可靠性，参照Abdalellah^[8]等的段塞流实验进行对比验证，数值仿真的边界条件依据实验条件设置，图3为实验与CFD模拟在相同管道段中段塞流形态的时序对比。其中，云图1表示气相，2表示液相，虚线3表示液塞边界。从液塞数量上看，实验和数值仿真的观测结果一致；从液塞形态上看，实验和数值仿真的观测结果基本吻合，说明了CFD模型的有效性。

3 结果与讨论 3.1 气液两相流流体激励特性 3.1.1 典型流型对比分析

为了更直观地说明段塞流诱导管道振动的剧烈程度，对于水平U型管，设置不同的气液两相流入口边界条件，研究分层流、波状流、环状流的管内流动特性，并与段塞流的管内流动特性作对比。4种工况的选取依照Baker流型图，如表3和图4所示。

选取U型管跨中横截面，对比4种流型下空泡分数的时域特性，如图5所示。

4种流型工况的混合流平均流速差别较大，仿真总时长设置不一。为便于展示，波状流和环状流参考上方横轴，段塞流和分层流参考下方横轴。由图5可知，分层流空泡分数稳定在0.74附近；环状流在流型充分发展后，空泡分数基本稳定在0.78附近，并伴随小范围波动；波状流相较环状流呈现出较大幅度的波动。段塞流具有与其他3种流型截然不同的特性，其空泡分数在0与0.75之间呈锯齿状反复交替，这与段塞流在管道中液塞与气体交替前进的特征相符。

根据式（15）和式（16），截面空泡分数频域特征可反映激励力频域特征，4种流型下跨中横截面空泡分数频域曲线如图6所示。可以看出，段塞流诱导激励力幅值最大，分层流诱导激励力幅值最小。4种流型在水平U型管中诱发的流体激励力都属于低频激励力，主频均在5 Hz以下。其中，段塞流在0.71 Hz处达到空泡分数最大峰值。

3.1.2 段塞流特征参数分析

设定气相折算速度为 $ {J_{{v}}} = 2.44{\text{ }}{\rm{m/s}} $ ，液相折算速度分别为 $ {J_{ {w}}} = 0.8\;{\rm{m/s}} $ ， $ 1\;{\rm{m/s}} $ ， $ 1.15\;{\rm{m/s }}$ ， $ 1.3\;{\rm{m/s}} $ ，分析4种段塞流工况下水平U型管中段塞速度、频率和长度3种特征参数与液相折算速度的关系，同时与部分学者水平直管段塞流试验进行对比^[8-10]，如图7～图9所示。

段塞速度是液塞在段塞流单元中的平均速度，其对液塞生长、发展和消失的机理研究起重要作用。从数值仿真结果中，每间隔0.5 s获取一次段塞流形态云图，取中部水平直管段捕获的所有液塞的平均速度为本次仿真的段塞速度（见图7），在气相折算速度不变的情况下，水平U型管中段塞速度随液相折算速度的增大而线性增大。

段塞流在管道中的流动是复杂的随机过程，并不具有规律的周期性，段塞频率只是用周期性来近似模拟段塞流流动。对中部水平管段的段塞进行追踪，在管道相同位置出现相同段塞时认为是段塞的一个周期。由图8可以看出，在气相折算速度不变的情况下，段塞频率随液相折算速度的增大而增大。此外，在气相和液相折算速度分别为2.44 m/s和1 m/s时，段塞频率为0.7 Hz，这与相同工况下获取的段塞流激励力主频0.71 Hz相吻合，表明段塞频率同流体诱导激励力主频是相对应的关系。

段塞长度和段塞频率是强相关参数，段塞频率增大，段塞长度减小。从数值仿真结果中，每间隔0.5 s获取一次管道段塞流形态云图，取中部水平直管段捕获的所有液塞的平均长度为本次仿真的段塞长度。由图9可以看出，在气相折算速度不变的情况下，随着液相折算速度的增大，段塞长度减小。

因管道形态、内径、壁厚等参数不一致，本文与其他学者的结果数值并不一致，但在特征参数随液相折算速度的变化趋势上是统一的，这揭示了段塞流包括段塞速度、频率及长度在内的特征参数变化规律在水平管道中具有一致性。

3.2 段塞流工况下管道振动特性分析 3.2.1 固体域模型及管道系统模态计算

管道采用钢材料，壁厚为20 mm，弹性模量为200 GPa，泊松比为0.3，密度为7850 kg/m³，生成网格如图10所示。管道两端设定固定约束，将管道内壁设为流固耦合面进行数据交互。

对管道进行干湿模态分析^[11]，前6阶结果对比如表4所示。由于水的密度和声速特性远大于空气，在湿模态分析时采用极端工况，即均为单相液态水，故气液两相流管道系统固有频率介于表中所列干湿模态固有频率之间。已知段塞流工况下空泡分数主频（流体激励力主频）为0.71 Hz，与管道固有频率的错开率大于78%，故流体与管道不会产生共振。

3.2.2 管道振动响应分析

为探究水平U型管在段塞流作用下的振动响应特性，设定液速 $ {J_{{w}}} = 1{\text{ }}{\rm{m/s}} $ ，气速 $ {J_{{v}}} = 2.44{\text{ }}{\rm{m/s}} $ ，在管道两弯头设置2个对称的监测点（见图1），提取其位移时域曲线，并进行频谱分析，如图11所示。

可知：

1）监测点1和监测点2的时域、频域曲线基本重合，两监测点在几何上对称，在振动响应上也对称。

2）X和Y方向峰值振动频率为6 Hz，接近管道二阶固有频率，Z方向峰值振动频率为3.8 Hz，接近管道一阶固有频率。这与管道固有振型相对应，根据模态分析结果，管道一阶振型主要沿Z方向，二阶振型主要沿X和Y方向。

3）两端固定约束的水平U型管在段塞流工况下响应幅值呈现出X方向最大、Z方向次之、Y方向最小的特性。气液两相流以较大速度流入管道，沿着XOY平面冲击原有流场使其逐渐趋于段塞流分布，且管道在X方向上跨距较Y方向更大，故X方向的响应幅值也较大；而管道一阶振型主要沿Z方向，其振动响应幅值介于X和Y方向之间。

4 结　语

本文基于计算流体力学和计算结构动力学方法，对气液两相流诱导水平U型管振动进行了数值仿真研究。在验证数值仿真框架可靠性的基础上，对比分析气液两相流典型流型，重点研究了振动最剧烈的段塞流的流体特征参数和其诱导的管道振动响应，主要结论如下：

1）基于Ansys几何重构VOF模型、 $ k - \varepsilon $ 湍流模型可以准确预报管道内气液两相流流型分布，结合流固耦合模块可实现对管道振动的数值仿真。

2）4种典型流型在水平U型管中诱发的流体激励力都属于低频激励力，主频均在5 Hz以下。段塞流诱导的激励力幅值最大，波状流和环形流次之，分层流最小。

3）对于水平U型管，在气相折算速度不变的情况下，随着液相折算速度的增大，段塞速度和段塞频率增大，段塞长度减小。此外，段塞流诱导激励力主频同段塞频率相互对应，二者从不同层面反映了段塞流诱导管道振动的原因。

4）U型管两弯头结构对称，振动响应也基本对称。X和Y方向峰值振动频率接近管道二阶固有频率，Z方向峰值振动频率接近管道一阶固有频率，与管道固有振型相对应。