基于数字孪生技术的河流穿越涡激振动评估

马剑林，连江桥，吴志锋，关佳智，熊 健，曲油伊

1.国家管网集团西南管道有限责任公司，四川成都 610000

2.中国石油天然气管道工程有限公司，河北廊坊 065000

据了解，河床侵蚀会使管道暴露，进而造成管道悬跨。这些暴露的无支撑管道会受到振荡力的作用，导致管道产生循环位移，通常称为涡激振动（VIV），这些涡激振动会导致管道疲劳损伤或结构过载。为此，管道行业开发了相关的分析工具来评估管道的VIV，由于这些工具是借助于流体动力学分析和实验研究而开发的，仅考虑了普通管道的情况，将管道作为一个理想的圆柱体考虑，但是许多老式管道在穿越河流时会加装配重块等附加结构，这类情况使得传统的分析工具难以适用。

数字孪生（Digital Twin）以数字化的方式建立物理实体的多维、多时空尺度、多学科、多物理量的动态虚拟模型[1]，以此来仿真和刻画物理实体在真实环境中的属性、行为、规则等。数字孪生技术作为实现数字化转型、解决物理世界与信息世界之间交互共融、践行智能制造目标的关键技术，利用其对管道本体的特征、行为、形成过程和性能等进行描述和建模，建立与现实管道完全对应和一致的虚拟模型，可实时模拟管道自身在现实环境中的行为和性能。

为了评估发生洪水时管道的完整性以及在高流速下河床侵蚀可能造成的管道悬跨问题，研究基于数字孪生体，采用新型的双向流固耦合作用（FSI）技术来评估现有经验模型未考虑的普通管道（仅圆柱体）情况。本文概述了已开发并成功应用于评估各种河流穿越结构完整性的技术，该分析技术已被用于考虑对一系列场景具有敏感性的管道响应，其中也考虑了其他因素的影响，如平均河流流速、穿越角度、河流宽度和未支撑的管道长度、管道支撑条件以及配重块的几何形状与尺寸和位置等。

1 VIV和新型评估程序开发的背景

当流体流过物体时，流体会产生涡流，因此改变了物体表面压力和剪应力的分布，使物体受到随时间变化的阻力和升力。物体表面压力和剪应力分布变化的频率由物体尾流的频率决定。如果让物体对这些水动力作出反应（位移），它就会振动。物体上的升力会引起横流振动，同样地，阻力会引起直流振动，这些振动被称为涡激振动（VIV）。如果振动涡旋力的频率与结构的固有频率相匹配，结构就会发生共振，结构的完整性就会受到不稳定振动的影响；当振力频率与结构固有频率不匹配时，结构的振动响应也能保持稳定。但是随着时间的推移，这种振动仍然会对结构造成疲劳损伤。

由于水流对河床的侵蚀，河流穿越位置的管道可能会产生暴露和悬跨，如图1 和图2 所示。河流穿越处的管道可以分为三部分，分别是水中无支撑管道、水中有支撑管道和土体中的管道，它们的特性会影响管道的振动响应。水流过暴露的管道会产生振荡的升力和阻力，升力会导致管道沿垂直方向的循环位移，称为横流振动；而阻力会导致管道沿水平方向的循环位移，称为直流振动。

图1 管道穿越河流示意

图2 河流穿越处无支撑的管道

大多数现有文献和设计标准[2]都是基于数值模拟和实验推导的数据，或者是已经建立了的半经验公式来预测管道VIV诱导振动，但是这些预测是基于理想圆柱体上的二维流动，只考虑了圆柱横截面的横向流动并允许刚体运动。在实际情况下不能达到理想的圆柱体情况，河流穿越处的管道大多数都有附加结构，而且河流流动是三维的，河流流动可能不垂直于管道。同时，管道的振动不遵循刚体运动规律而是可变形的。为了考虑这些与现有管道VIV分析模型的不一致影响，研究开发了一种更先进的分析技术，可以处理和评估每种情况下这些因素的影响。

2 流固耦合模拟

1）双向流固耦合（FSI）分析方法。为了探究考虑管道支撑、河流流向等因素的影响，而提出了该种分析方法。FSI 分析采用计算流体动力学分析（CFD）来评估管道载荷，采用有限元分析（FEA）来估计管道响应。如图3所示。该过程首先完成有限元分析，确定由于重力和浮力载荷引起的管道位移，然后进行CFD分析以确定管道上的力，再进行有限元分析，估计新的管道位置，CFD和FEA分析在时域内完成耦合，以更新管道位移引起的加载，并按照小时间步长重新评估管道的新响应。

图3 双向流固耦合作用过程

2）控制方程。通过管道的流动被视为非定常流、三维流、湍流、等温流和不可压缩流。计算流体力学中的控制方程取决于湍流模型的选择。对于所有的研究案例，发现在湍流区流动雷诺数（Re）下降。雷诺数计算公式为：

式中：ρ 为水的密度，kg/m³；V 为水的流速，m/s；D为管径，mm；μ为水的动态黏度，Pa·s。

由于管道上存在配重块等附加结构，在管道尾流中会出现较大的流动分离和产生多个频率。在这种情况下，使用传统的湍流模型如非定常雷诺平均湍流模型（U-RANS）是不合理的，主要是因为它无法解决具有附加结构的管道尾迹中的大多数频率。解析尾迹频率对于VIV 评估至关重要，考虑到计算时间和成本，研究得到了一种分离涡模拟方法（DES）的混合湍流模型。DES 湍流模型在靠近水的计算单元中使用定常雷诺平均湍流模型（RANS 模型）对结构或域边界界面（即壁面）进行流体边界层求解，在远离壁面的计算单元如管道的尾流中使用大涡模拟（Large Eddy Simulation，LES）等高保真湍流模型[3]。DES 结合了URANS 和LES 湍流模型的优势，根据计算网格的大小和模拟的时间步长，捕获了管道尾流中的大部分频率。本研究采用k-ω 基于剪切应力传输（Shear Stress Transport，SST）的DES 模型，其中k为湍流动能，ω 为湍流频率。在这里，DES 模型的质量守恒和动量守恒方程记为：

LES 湍流模型采用Smagorinsky-Lily 亚网格尺度（SGS）模型[4]。有限元法的控制方程为：

式中：mij为结构质量，kg；为结构加速度，m/s2；Kij为刚度矩阵；uj为结构位移，m；Fb为外力，N。

分别从FEA 分析和CFD 分析中计算出的位移和力，使用外部耦合求解器进行交换。采用一对一映射或插值技术，将有限元分析中计算的位移插值到CFD 域的流体结构界面。然后利用位移来更新流体域，并利用动态网格技术对界面进行相应的位移。类似地，在有限元分析域的流体结构界面上，通过CFD 求解器计算得到的力被映射或插值。为了获得最大程度的一对一映射，减小插值误差，两种求解器的流体结构界面元素保持相似并共享。

3 普通管CFD和FSI分析的数值设置

作为对CFD/FSI 模型的验证，将FSI 分析计算的结果与现有半经验模型的结果进行比较，验证考虑了普通管道上的VIV 问题。将FSI 结果与CFD结果（不耦合结构求解器）以及标准模型（如DNV-RP-F105）计算的结果进行比较。

3.1 问题描述

对于直径为508 mm、河流流速为4 m/s、与河流夹角为90°的管道，对暴露无支撑河流穿越的普通管道进行CFD 分析的样本如图4所示。水流角度定义为管道轴线与水流方向之间的角度，见表1中案例一状况。

图4 CFD和FSI分析用普通管/m

表1 研究案例及其属性

3.2 计算网格和边界条件

在普通管道周围的流体包络线中生成四面体元素，在FSI 计算过程中，四面体单元具有适应大位移的灵活性。元素聚集在管道外表面以捕获边界层。第一个元素的大小为y+＜5，以满足“壁面定律”条件。在尾流区域生成更细的元素来捕获柱体的尾流。在FSI 计算过程中，流体网格动态更新以适应管道的位移。

同样，在有限元模型中，管道以与四面体单元进行网格划分，在流体结构界面处存在一对一映射，如图5所示。

图5 有限元网格

在CFD 域的入口施加4 m/s 的河流流速，出口施加大气压力。在对称边界条件下应用域的侧壁；在自由滑移边界条件下，应用域的上下壁。假定管道外表面光滑，流体速度为零（无滑移状态）。普通管道的管道尾流中可能不会出现较大的流动分离，流动可以用更简单的RANS 湍流模型来解决，因此可以构建k-ω SST 模型。在FSI 计算中，将管道外表面与有限元分析求解器耦合，并在有限元分析模型中采用固定-固定端条件。

4 渡河管CFD和FSI分析的数值设置

采用双向FSI 分析技术分析河流穿越场景，与标准化分析技术进行比较，并考虑现有VIV评估模型中未包括的因素，如管道上存在配重块。本文讨论以下两种渡河带有配重块的情景。

4.1 问题描述

如表1 所示，本文研究了三例河流穿越案例。第一个案例是裸露的普通管道（无附加结构），用于将FSI 分析与传统的VIV 评估进行比较。案例二和案例三被用于不同条件下，使用FSI 评估技术来研究具有配重块的管道跨越水流的场景。

4.2 研究案例及其属性

研究案例及其属性见表1。

4.3 计算网格和边界条件

在无支撑的管道跨度周围创建流体包络线，并在流体域中生成四面体元素，如图6所示。

图6 用于CFD和FSI分析的CFD网格

同样，在FSI 接口上建立了一一对应的有限元分析模型。在河道中生成实体四面体单元，在钢管中采用三角形壳单元，在河床上的管段和埋在土壤中的管段采用线单元。在线单元的轴向、横向和横向节点上附加弹簧单元，模拟管土相互作用。有限元分析模型如图7所示。

图7 有限元网格

河流流速如表1 所示，是在流体域入口处定义的。在CFD 域的出口定义大气压力。在流体域的顶部和底部分别施加自由滑移条件。域侧壁采用对称边界条件，无支管跨外表面采用光滑条件，这段管道被定义为流体结构界面。在有限元分析模型中，将管端定义为固定-固定端。模型中定义的运行参数如表1 所示。计算了重力引起的虚拟加速度，并考虑了浮力的影响。

所使用的土壤弹簧的特性如图8所示。

图8 土壤特性

5 结果与讨论

5.1 普通管无配重(案例一)

对于普通管道（无附加结构）情况，在进行双向流固耦合（FSI）分析之前，对有限元分析模型进行了模态分析，计算了振型和相关的固有频率。模态分析结果表明，该结构在直流模式和横流模式下的固有频率均为7 Hz。

在将CFD 求解器与FE 分析求解器耦合进行双向FSI 分析之前，仅进行CFD 分析，将CFD 分析结果与FSI 分析计算结果进行比较。进行CFD 分析的另一个原因是初始化FSI 计算，这有助于减少计算时间，更快地达到所需的解。图9 显示了圆柱体尾流中产生的涡流。这些反向旋转的双涡被广泛地称为卡门涡，并以一定的频率脱落，通常以一种称为斯特劳哈尔数（St）的无维形式表示，其形式为：

图9 涡流脱落

式中：fs为脱落频率，Hz；D 为直径，m；U 为流速，m/s。

图10 中给出的图形显示了由于这些涡的脱落而对圆柱体的升力和阻力产生的波动作用。从这些力的频谱分析，脱落频率为2.55 Hz，如图11所示。这种情况下的雷诺数是3E6，使用式（5）计算得到St的值为0.32，对于这个雷诺数，公开文献中也记述过St的值数为0.32。

图10 普通管上的波动升力和阻力（CFD分析）

图11 升力谱分析（CFD分析）

从CFD 分析结果可以明显看出，旋涡脱落频率2.55 Hz 与结构固有频率7 Hz 并不接近。从这个观察来看，就VIV 而言，可以注意到振动是稳定的，因为结构不会发生共振。然而，使用CFD 分析并不能确定地进行评估。通过双向FSI 分析，观察允许结构响应升力和阻力波动的效果，分析表明在FSI计算过程中观察到涡旋结构没有变化。

在管道中心的节点上跟踪管道位移。管道位移时程如图12 所示。发现振动是稳定的，因为振动的振幅没有指数增加。DNV RP-F105 使用的参数，减速速度V*定义为：

图12 从FSI分析中捕获的振幅

式中：fn为结构固有频率，Hz；D 为管径，m；U为流速，m/s。

然后使用半经验方程或图形方法计算出降低速度，以计算横流和直流方向上的振动幅度。将DNV RP-F105 模型预测的振幅与FSI 结果进行比较，如图12 所示。结果表明，DNV RP-F105 在直流方向上估计振动幅值为8 mm，在横流方向上没有预测振动。然而，双向FSI 分析表明，在流入方向上，振动幅度大于8 mm，而且观察到横流方向也有一定的振动，如图12 所示。基于这一结果，对于这一特定场景，DNV RP F105模型对直流和横流振幅的估计都较低。

此外，与仅CFD 分析中仅识别出一个脱落频率不同，双向FSI 分析的升力谱分析显示存在多峰，其中两个频率占主导，如图13所示。这表明，CFD 分析和双向FSI 分析的结果可能不同，可能导致不同的VIV 评估结果，在文献[5]中也有类似的结论。造成这种情况的原因可能是，当管道在水中位移时，附加质量等影响在CFD 分析中没有体现出来，但被认为是VIV评估的重要因素。

图13 升力谱分析（FSI分析）

为了保证所采用的FSI 技术能够捕捉到不稳定的振动场景，将普通管道的固有频率降低到接近旋涡脱落频率（2.55 Hz），本文通过创建一个细长的管道，在管道悬跨段两侧增加10 m 来实现，增加的管段不参与CFD 求解，只包含在有限元分析模型中。管道的几何形状如图14 所示，添加的管道段采用线元建模，以减小有限元模型尺寸。

图14 不稳定振动场景的有限元模型

为了保证计算过程中网格在受到不稳定振动引起的较大位移后仍能保持计算质量，在计算过程中对网格进行了连续监测。图15 显示了计算过程中的变形网格，并显示了计算过程中捕获的涡流。

图15 分析过程中变形的网格和捕获的涡流

FSI分析能够预测不稳定振动，结果如图16 所示。在管道悬跨段两侧各增加10 m 的管道，使结构的固有频率和脱落频率更接近，从而促进了不受控制的振动，管道的垂直运动不受控制。这种不受控制的横流振动影响了直流振动。如图16 中横流位移所示，不受控制的振动由管道位移的指数增长表示。这个结果证实双向FSI 分析技术能够捕捉到这样的场景。

图16 不稳定振动场景下的直流振动和横流振动振幅

5.2 河流交叉应用程序（案例二和案例三）

在进行FSI 计算之前，先进行模态分析，然后进行CFD 分析。进行模态分析的目的是计算模态振型及相关频率。进行CFD 分析主要是为了初始化FSI 计算和估计尾流频率。对于案例二，升力在0.5 Hz 和1.25 Hz 频率波动，具有显著的谱密度。阻力在0.1 Hz 和0.4 Hz 频率波动。同样，对于案例三，升力在主频为0.96 Hz 时波动，阻力在主频为0.39 Hz 和1.71 Hz 时波动，在案例三中的力波动频率高于案例二。但是，仅通过比较固有频率和脱落频率，无法对VIV 进行评估，因为CFD 分析本身无法捕捉到附加质量等影响，而附加质量在VIV中起着重要作用。因此，双向FSI 分析变得非常重要。进行双向FSI 分析以估计两种情况下的VIV 振幅。利用q 准则定义的等面图来识别尾流中的涡，q准则是速度梯度张量的第二不变量，其形式为:

在尾流中发现了不同的湍流相干结构，如“马蹄形”涡。由于河流权重的存在导致分离，并且由于DES 模型能够捕获此类湍流结构，因此在尾流中识别出了这些涡，并捕获了相关频率。

图17 显示了振动的直流振幅和横流振幅。图中对于横流位移来说，曲线负值表示结构自重引起的位移。在案例二中，振动幅度保持不变，连续几次循环管道位移为0.17 m，表明振动情景稳定。在这种情况下，横流振动不太显著，且随时间衰减。相反，在案例三中，直流振幅随时间呈指数增长，表明振动不稳定。横流振幅随时间逐渐衰减。从CFD 分析中的力谱分析可以看出，案例二中的力波动频率高于案例三，这可能是导致案例三振动不稳定的原因。由于河流的重量和急剧的水流角度的存在导致了流动的不对称，另外“马蹄形”的漩涡等不稳定性也可能导致案例三中的不稳定振动场景。案例二中的管道由于受到具有较大位移的稳定振动，可以完成疲劳寿命评估，以了解稳定振动造成的威胁。第三种情况的振动是不稳定的，这种情况被认为是最终失效的情况。对于案例二，对有限元分析结果进行后处理，评估应力场，并在观察到最大等效应力的节点上计算轴向、径向和环向应力的范围。计算得出主应力范围为209.5 MPa。出于保守性考虑，假设环焊缝位于最高应力位置，并使用平均曲线法完成疲劳寿命评估，根据振动频率，连续高流量条件下的结构疲劳寿命估计为4 d。

图17 用于FSI分析的直流和横流振动幅值

5.3 灵敏性分析

为了了解温度差和管道在河床上的长度等参数对结构固有频率的影响，研究使用模态分析进行+敏感性研究。本研究考虑了表1 中描述的案例三的管道几何形状。基本参数情况如表2所示，其中不考虑温差，在基本参数情况下悬跨两侧河床上的自由悬浮管道长度为15 m。

表2 基本参数

表3 为模型中考虑正、负温差时结构的固有频率。当模型中考虑负温差时，产生热拉应力，直流式和横流式频率分别从2.29 Hz 和2.36 Hz 增加到2.42 Hz 和2.62 Hz。同样，当考虑正温差时，直流频率由2.29 Hz减小到2.23 Hz，横流频率由2.36 Hz减小到2.27 Hz。从本研究中发现，考虑正温差进行VIV 评估是保守的，因为其降低了结构的固有频率。

表3 温差对结构固有频率的影响

最后，观察不同管长停在河床上效果的敏感性分析，结果如表4所示。横流模式下，管道长度对固有频率的影响较小，管道长度的增加会降低横流模式下的固有频率。在河床上停留的管道长度对直列模态固有频率没有影响。研究发现，就结构固有频率而言，在河床上停留的管道长度对管道过河VIV评价的影响不显著。然而，完全忽略它（假设河床上没有管道）可能会增加结构的固有频率。

表4 位于河床上的管道长度对结构固有频率的影响

6 结论

介绍了一种先进的基于数字孪生双向流固耦合（FSI）技术，表明该技术能够考虑其他传统分析方法（如带有配重块的河流穿越）难以解决的一系列河流穿越场景，并为VIV提供更好的管道响应预测，以支持完整性评估。证明数字孪生技术具有巨大的应用潜力，可应用于模拟、分析、优化、监控与预测管道在河流穿越环境中的系列行为与过程。通过模拟管道在河流穿越环境中的运行情况，掌握要素行为、状态等运行参数，监控管道运行等情况，同时通过预测模型进行质量管理，提高质量指标预测与关键控制指标预测精度，实现先进控制和完整性管理的长期有效应用。河流穿越流固耦合作用分析过程和评估过程所完成的工作已经证明：

1）与普通CFD 相比，双向FSI 分析可能会产生不同的结果，因为在CFD 分析中无法捕捉到附加质量等影响。双向FSI 分析通过在时域内同时捕捉流体动力学和结构响应的相互作用，可以模拟接近实际情况。与普通CFD 相比，双向FSI 分析结果更具有参考价值。开发了一种双向FSI 分析的标准方法，该方法可以考虑现有经验模型所不包括的管道内部涡动的几何细节。

2）该技术可以识别管道中的不稳定振动和稳定振动。管道的完整性可以根据双向FSI 分析获得的振动类型进行评估——疲劳寿命（稳定振动）和最终失效（不稳定振动）。

3）该技术可以考虑一系列过河场景，例如：有或没有附加结构、管道尺寸、加压或非加压管道、流体流速、不同的河流宽度。

本文将FSI 分析技术与模型测试或全尺寸数据进行比较，将考虑不同流动角度的FSI 分析结果与传统半经验模型计算结果进行比较，揭示了FSI 结果对河流流动梯度和管道下方流道大小的敏感性。