考虑跨肩处动态沉降的近海床悬跨管道耦合顺流向与横向涡激振动特性

《Marine Structures》：Coupled In-Line and Cross-Flow Vortex-Induced Vibration Characteristics of Near-Seabed Free-Spanning Pipelines Considering Dynamic Settlement at Span Shoulders

　　

摘要
海流可在跨流 (CF) 和顺流 (IL) 方向诱发自支撑海底管线的涡激振动 (VIV) 耦合现象，从而在悬跨肩部引发复杂的管线 - 土壤相互作用 (PSI)。管线振动会降低土壤强度并挖蚀悬跨肩部周围的土壤，从而加速管线沉降。管线沉降后，间隙比（定义为悬跨与海床之间的距离归一化于管径）减小，由于边界效应，改变了管线的 VIV 特性。本研究建立了一种综合方法，用于模拟自支撑管线的 CF 和 IL VIV，包含非线性 PSI 和边界效应。首先，使用耦合欧拉 - 拉格朗日 (CEL) 技术研究 PSI，并提出了一种适用于横向振荡管线的新型动态沉降模型。随后，基于矢量形式固有有限元 (VFIFE) 方法建立了自支撑管线模型，集成了非线性 PSI、动态沉降和尾流振荡器模型。参数研究评估了间隙比、轴向张力、动态沉降和悬跨长度的影响。结果表明，减小间隙比显著抑制了一阶和二阶模态下的 IL 振幅，同时显著降低了一阶和过渡模态下的 CF 振幅。悬跨肩部管线的横向振荡会导致额外的沉降，特别是在二阶模态下，振荡在更大范围内引起更深的沉降。

1. 引言
1.1 背景
复杂的海底地形可能导致海底管线出现悬跨，海流引起的冲刷可能在运行期间进一步加剧悬跨的发展。当海流流经自支撑管线时，会产生交替的涡旋脱落。当涡旋脱落频率接近管线的固有频率时，可能触发共振，导致严重的疲劳损伤。此外，间隙比、复杂的管线 - 土壤相互作用 (PSI) 和多个悬跨进一步复杂化了涡激振动 (VIV) 响应。在南中国海，地形和地貌特征非常复杂。遥控操作车 (ROV) 调查的检查数据显示，在 40 至 72 米水深处的软粘土海床上，悬跨段出现的频率日益增加 [1]。此外，悬跨长度和间隙比都在逐渐增加，对管线安全构成了严重威胁。因此，必须考虑非线性 PSI 来研究自支撑管线的 VIV 特性，以便为工程应用提供可靠数据。

研究人员通过数值模拟和实验研究调查了自支撑管线的 VIV 特性，重点关注 VIV 响应和抑制 [[2], [3], [4], [5]]、疲劳损伤 [[6], [7], [8]] 以及多悬跨 VIV 响应 [[9], [10], [11]]。研究表明，自支撑管线的 VIV 在顺流 (IL) 和跨流 (CF) 运动之间表现出复杂的耦合 [2]。在多悬跨管线中，不同悬跨之间存在显著相互作用，动态响应随悬跨长度变化 [9,10]。此外，悬跨比的变化显著影响相邻悬跨的疲劳损伤分布。随着相邻悬跨之间长度差异的增加，较长悬跨中的疲劳损伤显著加剧 [11]。

此外，悬跨肩部的非线性 PSI、间隙比和轴向张力共同导致了更复杂的振动特性。悬跨肩部的 CF 振动导致海床发生弹塑性变形，导致土壤强度退化和管线沉降，进而影响管线的固有频率和 VIV 特性。悬跨肩部不同的土壤条件对振动有显著影响 [12]。粗糙的海床会在自支撑管线中引发高水平的疲劳损伤，在 CF 和 IL 方向上均发生不可接受的损伤 [13]。非线性迟滞土壤模型也显著影响振动特性，特别是改变了悬跨肩部管线的弯曲应力 [14]。减小间隙比通常会导致圆柱体和管线的 CF 振幅降低 [[15], [16], [17], [18]]，并增加了发生一阶锁定的流速 [16]。对于柔性管线，减小间隙比导致模态转换的折算速度范围更宽，并增强了空间频率竞争 [18]。此外，轴向张力激励在较低折算速度下显著影响自支撑管线的疲劳损伤 [19]。考虑中面拉伸效应会增加系统的固有频率，当强迫振动达到临界阈值时，会发生倍周期分岔，影响管线疲劳损伤 [20]。

总之，自支撑管线的 VIV 特性受多因素相互作用支配，但当前研究主要侧重于孤立方面。例如，常见的数值模拟采用尾流振荡器模型，仅考虑 CF 振动而忽略 IL 分量。当通过有限差分法离散欧拉 - 伯努利梁时，这些模型通常仅计算恒定轴向张力下的 CF 振动，无法捕捉时变轴向张力效应。在悬跨肩部 PSI 建模中，简化方法仍然普遍，土壤通常使用线性弹性或理想弹塑性模型表示——甚至过度简化的简支/固定端条件。这种方法无法模拟沟槽形成、土壤软化或动态管线沉降。尽管一些实验研究使用粘土或沙子进行肩部支撑，但它们未能考虑 PSI 中的尺度效应。由于土壤强度和刚度受有效应力支配，单重力条件下的缩尺模型试验无法准确再现满尺寸行为 [21]。相反，离心实验虽然能够模拟真实的土壤应力，但在重现真实流场方面面临挑战。此外，自支撑管线 VIV 机制的可靠相似理论仍不完善，阻碍了实验模型与原型管线之间建立合理的对应关系 [10]。当前大多数实验进行无尺寸模型试验，在受控流条件下（特别是亚临界湍流）调查特定质量比下的振动特性，同时保持相同的斯特劳哈尔数并系统地改变折算速度。

1.2 研究目标
鉴于现有研究的局限性，本研究旨在开发一种计算高效的数值方法，用于自支撑管线的满尺寸 VIV 分析，能够准确捕捉 IL 和 CF 振动的耦合、PSI 以及边界效应。管线结构建模的有限元方法及相关离散求解算法应提供稳健的收敛性和可扩展性，以便能够轻松纳入其他影响因素。

已进行了广泛的圆柱体和管线 VIV 实验研究 [[22], [23], [24], [25], [26]]，并辅以采用计算流体动力学 (CFD) 和经验模型的数值研究 [[27], [28], [29], [30]]。实验方法提供了对 VIV 机制的直观见解，而 CFD 方法通过求解纳维 - 斯托克斯方程提供高保真度水动力学预测，尽管计算成本显著更高。经验模型——特别是广泛采用的尾流振荡器模型——由于计算效率高，在工程应用中备受青睐。Facchinetti [31] 基于实验数据对尾流振荡器模型中的经验参数进行了校准和优化，表明加速度耦合是范德波尔方程与结构振动方程之间的最佳策略。Xu [32] 通过圆柱体模型试验完善了这些参数。Qu [33] 开发了一种用于弹性支撑刚性圆柱体耦合 CF 和 IL VIV 的新型尾流振荡器模型，表明与实验数据吻合良好。Fehér [34] 通过三个新系数和耦合项进一步将模型的适用性扩展到边界效应，保持与实验基准的高度准确性。因此，本研究采用了 [34] 中增强型的尾流振荡器模型。

已进行了大量关于 PSI 的实验、理论和数值研究 [21,[35], [36], [37]]。满尺寸和缩尺实验主要关注钢悬链立管触地点 (TDZ) 的 PSI。还进行了理论和数值研究以调查 PSI 机制，包括基于弹塑性理论的方法和使用有限元方法 (FEM) 的土壤单元建模。然而，基于 FEM 的方法通常涉及显著更高的计算成本。离心测试对缩尺模型施加比重力加速度大几十倍的加速度，从而能够真实模拟满尺寸应力条件。该技术已被广泛用于调查圆柱体 - 土壤相互作用机制，并基于阻力 - 穿透曲线 (P-y 曲线) 开发经验 PSI 模型。最普遍的垂直 PSI 模型是 AB 模型 [38] 和 RQ 模型 [39]，两者均是基于大量实验数据经验得出的。尽管 AB 和 RQ 模型都包含了重新穿透阶段期间的土壤强度退化，但它们的能力不同：AB 模型模拟沟槽形成但未能捕捉由垂直循环加载引起的强度退化，而 RQ 模型保留原始海床高程作为每个重新穿透周期的起点，因此忽略了沟槽形成。You [40] 通过包含土壤强度退化扩展了 AB 模型，提出了重新穿透阶段的经验公式，称为 ABY 模型。Zargar [41] 开发了一种新型迟滞模型，捕捉土壤刚度退化和沟槽形成，这是通过离心试验完成的。Chen [42] 在 CABLE3D 程序的非线性 PSI 模块中实现了类似的机制。横向 PSI 的研究主要集中在大振幅横向运动期间的突破阻力、残余阻力和土壤棱体形成。在管线横向屈曲的研究中，通常采用三线性及类似的 P-y 曲线模型来模拟长管线中的 PSI [[43], [44], [45], [46]]。

此外，以前的研究表明，小振幅横向振荡可以显著增加管线沉降。Lund [47] 分析了 Zeepipe 2B 管线项目的现场数据，发现安装期间的横向运动导致沉降远远超过设计预测。实验证据进一步证实，振幅小至 0.05D 的横向振荡产生的沉降量可达静态加载引起沉降量的八倍。尽管有这些发现，关于横向振荡影响的研究仍然有限，主要侧重于大振幅运动。White [48] 提出了循环横向运动期间土壤棱体的经验模型，包含了前进管道挖走的土壤厚度。Cheuk [49] 基于离心试验和流动规则理论，开发了恒定垂直加载下管线的动态沉降模型。Dutta [50] 采用耦合欧拉 - 拉格朗日 (CEL) 技术模拟小振幅循环运动，通过离心试验比较验证了 CEL 在预测沉降方面的准确性。随后，Dutta [36] 确认带有用户子程序的 Abaqus CEL 能精确捕捉应变率和应变软化对土壤剪切强度的影响。小振幅横向振荡下管线沉降的演化是本研究的关键焦点。

采用矢量形式固有有限元 (VFIFE) 方法开发了自支撑管线的有限元模型。该方法最初由 Ding [51] 提出，采用矢量力学和点值描述将结构离散为质量粒子。动态响应使用平移和旋转微分方程求解，无需组装结构刚度矩阵——从而避免了矩阵奇异问题。由于其在处理大变形、碰撞、穿透和断裂等不连续行为方面的能力，VFIFE 在管线分析中显示出强大的适用性 [[52], [53], [54]]。此外，VFIFE 方法也已广泛应用于深海采矿垂直运输系统的动态分析 [[55], [56], [57]]。基于 VFIFE 的三维模型能够实现 IL-CF 耦合振动分析和时变轴向张力计算，同时允许实施柔性边界约束和计算管线端部的轴向位移。当与非线性 PSI 和涡激力集成时，该方法保持高度的数值稳定性和收敛性，并保持易于扩展以供未来研究。

本研究始于 CEL 技术，以研究小振幅横向振荡下管线沉降的演化。推导了流动规则公式，提出了一种变垂直加载下的新型动态沉降模型，同时使用可用实验数据完善了垂直非线性 PSI 模型。在此基础上，基于 VFIFE 方法在 Fortran 中开发了管线模型，集成了非线性 PSI 模型和提出的动态沉降模型。使用尾流振荡器模型模拟 IL 和 CF 方向的 VIV。随后检查了自支撑管线的 VIV 特性，特别强调间隙比、轴向张力和悬跨长度的影响。提出的框架提供了一种综合方法，用于考虑非线性 PSI 和海床边界效应的 IL-CF 耦合 VIV 满尺寸模拟。论文内容概述如下：第 1 章（引言），第 2 章（理论基础），第 3 章（模型验证与参数拟合），第 4 章（结果与讨论），第 5 章（结论）。

2. 理论基础
2.1 用于管线建模的 VFIFE 梁单元理论
2.1.1 控制方程
VFIFE 梁单元的理论框架简要总结如下，详细公式见 [51]。

如图 1 所示，VFIFE 方法将结构离散化为通过规范单元相互连接的有限质量粒子。粒子运动和单元变形共同表征结构行为并计算内力。连续时间步之间粒子的运动由被称为路径单元的标准化方程组控制。采用中心差分格式，质量粒子的控制方程由牛顿第二定律推导得出：(1){xn+1=c1Δt2[M+Mf]?1(Psoil(t)+Phydro(t)+Pgravity(t)+Pbuoy(t)+f(t))+2c1xn?c2xn?1θn+1=c1Δt2I?1(Qsoil(t)+Qhydro(t)+Qgravity(t)+Qbuoy(t)+m(t))+2c1θn?c2θn?1c1=1/(1+μdΔt/2),c2=c1(1?μdΔt/2)其中，xn+1、xn 和 xn?1 表示粒子在三个连续时间点的坐标向量，θn+1、θn 和 θn?1 表示粒子在相应时刻的方向角向量，M 和 Mf 分别是结构质量矩阵和附加流体质量矩阵，f(t) 和 m(t) 分别是内力和力矩向量，μd 是阻尼比，Δt 是时间步长。Psoil(t)、Phydro(t)、Pgravity(t) 和 Pbuoy(t) 分别是由土壤阻力、水动力、重力和浮力引起的等效节点力向量。Qsoil(t)、Qhydro(t)、Qgravity(t) 和 Qbuoy(t) 分别是由土壤阻力、水动力、重力和浮力引起的等效节点弯矩向量。M 和 Mf 计算如下：(2)M+Mf=diag(m+mf,m+mf,m+mf)m=∑j=12(ρpipeAlj/2),mf=∑j=12(ρwaterCmAlj/2)其中，m 是管道单元质量，mf 是附加流体质量，A 是管道单元截面面积，lj 是单元长度，ρwater 和 ρpipe 分别是海水密度和管道密度，Cm 是附加质量系数。下载：下载高分辨率图像 (171KB)下载：下载完整尺寸图像图 1. 点值描述和路径单元。

惯性矩矩阵 I 计算如下：(3)I=∑j=12ΩTdiag(I1jI2jI3j)ΩI1j=(πDout4/32?πDin4/32)×ρpieplj/2,I2j=I3j=(πDout4/64?πDin4/64)×ρpieplj/2其中，Ω 是从全局坐标系到局部坐标系的转换矩阵。Dout 是外径，Din 是内径。

控制方程本身排除了刚度矩阵，从而避免了非线性分析中常见的收敛困难。为了实现迭代求解，必须显式计算载荷项。如公式 (1) 所示，平移位移和旋转位移均包含刚体运动分量和变形分量。该公式自然适应大位移、有限旋转和有限变形，无需专门的数值技术。内力 f(t) 和力矩 m(t) 通过反向运动过程确定，如 2.1.2 节所述，以提取纯变形。作用在管道单元上的垂直土壤阻力见 2.2 节，水动力见 2.4 节，横向土壤阻力见 2.5 节。获得作用在单元上的分布载荷后，使用 2.1.3 节描述的方法计算相应的等效节点力 Psoil(t)、Phydro(t)、Pgravity(t) 和 Pbuoy(t)，以及等效节点弯矩 Qsoil(t)、Qhydro(t)、Qgravity(t) 和 Qbuoy(t)。

2.1.2. 内力

图 2 说明了路径单元 ta?tb 内梁单元 1-2 的内力变化计算，其中 ta 和 tb 分别表示初始时间和最终时间。最终单元构型 1b-2b 受到反向运动以去除刚体运动，包括平移 ?u1 和旋转 ?γt。此变换得到纯变形，作为虚拟构型 1d-2d 与初始构型 1a-2a 之间的差值。平移向量定义为：u1=xb1?xa1，其中 xa1 和 xb1 分别表示时间 ta 和 tb 时节点 1 的坐标向量，对应于连接到单元 1–2 的粒子 1 的坐标向量。在全局坐标系中表示的旋转向量 γt 为：(4)γt=β^11ea1+θtretr=(β^1·ea1)ea1+sin?1(|ea1×eb1|)ea1×eb1|ea1×eb1|其中，β^11ea1 表示轴向旋转向量，θtretr 表示平移旋转向量。β^1=θb1?θa1={β^11β^21β^31}T 表示节点 1 的旋转向量。θa1 和 θb1 分别是时间 ta 和 tb 时节点 1 的方向角向量，对应于连接到单元 1–2 的粒子 1 的方向角向量。ea1 和 eb1 分别是时间 ta 和 tb 时轴 1-2 的单位向量，ea1=(xa2?xa1)/|xa2?xa1|，eb1=(xb2?xb1)/|xb2?xb1|。xa2、xb2 是时间 ta 和 tb 时节点 2 的坐标向量。下载：下载高分辨率图像 (192KB)下载：下载完整尺寸图像图 2. 反向运动。

在局部坐标系中表示为：(5)γ^t=θ^1ea1+θ^2ea2+θ^3ea3=β^11ea1+θtr(etr·ea2)·ea2+θtr(etr·ea3)·ea3其中，θ^1、θ^2 和 θ^3 是 γ^t 的分量；引入辅助参考点 3a 以定义平面 1a-2a-3a，且单位法向量 ea2 满足：ea2=(e13×ea1)/|e13×ea1|；e13=(xa3?xa1)/|xa3?xa1|，xa3 是节点 3a 的坐标向量；单位法向量 ea3 满足：ea3=ea1×ea2。正交归一基 ea1,ea2,ea3 构成路径单元 ta?tb 内单元 1-2 的局部坐标系。这些向量仅在初始路径单元中如此定义，对于后续路径单元，它们由旋转向量 γt 计算得出。

变形如下：(6)Δl=|xb2?xb1|?|xa2?xa1|(7)φ^j={φ^1jφ^2jφ^3j}T={β^1j?θ^1β^2j?θ^2β^3j?θ^3}T,j=1,2其中，Δl 是单元长度变化，φ^j 是节点 1 和节点 2 的角位移变化，β^2=θb2?θa2={β^12β^22β^32}T 是节点 2 的旋转向量，θa2 和 θb2 是节点 2 的方向角向量。

确定纯变形后，路径单元内的内力变化基于梁弯曲理论和静力平衡计算：(8)Δf^j={Δf^1jΔf^2jΔf^3j}T,Δm^j={Δm^1jΔm^2jΔm^3j}T(j=1,2)(9)Δf^12=EAljΔl^,Δm^12=GI^1ljφ^12,Δm^21=EI^2lj(4φ^21+2φ^22)Δm^22=EI^2lj(2φ^21+4φ^22),Δm^31=EI^3lj(4φ^31+2φ^32),Δm^32=EI^3lj(2φ^31+4φ^32)Δf^11=?Δf^12,Δm^11=?Δm^12,Δf^22=?1lj(Δm^31+Δm^32)Δf^32=1lj(Δm^21+Δm^22),Δf^21=?Δf^22,Δf^31=?Δf^32其中，Δf^j 和 Δm^j 是局部坐标系中表示的单元内力和力矩增量向量，E 是弹性模量，G 是剪切模量。

时间 tb 时的单元内力通过以时间 ta 时的参考内力 f^aj 和 m^aj 加上增量 Δf^j 和 Δm^j 进行增量累加，然后通过正向运动校正方向来确定。

2.1.3. 等效节点力和弯矩

在时间 tb，作用于梁单元 1b-2b 上的分布外力 p 必须在反向运动后施加于其虚拟构型 1d-2d，并将其分量分解到局部坐标系中作为 p^。(10)p^=Ω[Rt(γt)]Tp[Rt(γt)]T=Rt(?γt)=I3+sin(?γt)At+(1?cos(?γt))At2p^={p^1p^2p^3}T其中，Rt(γt) 是反向运动的转换矩阵，I3 是单位矩阵，At 是基于向量 γt 的矩阵，γt 是 γt 的模。

局部坐标系中的等效节点力和弯矩向量如下：(11)p^1=p^2=l2{p^1p^2p^3}T,?q^1=q^2=l212{0p^3p^2}T其中，p^1、p^2、q^1 和 q^2 是单元 1-2 施加给节点 1 和节点 2 的等效节点外力和力矩。

等效节点力和弯矩向量被转换到全局坐标系，然后通过正向运动返回到时间 tb 时的单元位置。

2.2. 垂直非线性管土相互作用建模

Li [58] 和 Chen [42] 比较了 AB、RQ 和 ABY 模型在循环载荷下的性能，均指出 ABY 模型中退化土壤刚度的极限值过低。因此，本研究采用 ABY 模型，并采用 Zargar [41] 基于实验数据提出的土壤强度退化曲线，该曲线更准确地表征了循环载荷下土壤强度的退化过程和极限状态。此外，文献 [21] 中的实验数据已用于修改最大土壤吸力公式。

当前模型将垂直管土相互作用 (PSI) 过程分为五个不同阶段（图 3）：O-A，初始穿透（骨架曲线）；A-B，弹性回弹；B-C，部分管土分离；C-O，完全分离；C-E，再穿透。z 是穿透深度，z1、z2、z3 和 zcn 分别是点 A、B、C、E 处的穿透深度。P 是单位长度土壤阻力，P1、P2、Pdn、Pcn 分别是点 A、B、D、E 处的土壤阻力。下载：下载高分辨率图像 (195KB)下载：下载完整尺寸图像图 3. 垂直非线性 PSI 模型。

骨架曲线基于 RQ 和 AB 模型使用的公式构建，并补充了浮力项：(12)Pb=NPsuDout+fbAsγsoilDoutDout2su=a(zDout)b(su0+sukz)Dout+fbAsγsoilDoutDout2(su0+sukz)其中，Pb 是骨架曲线上的单位长度土壤阻力，NP 是无量纲承载力因子，a 和 b 是拟合常数，su 是土壤剪切强度（在海床表面为 su0），suk 是强度梯度，γsoil 是有效土重，fb 是浮力因子，As 是位于泥面以下的管道名义面积。

第 n 次循环后的土壤阻力公式如下 [41]：(13)Pdn+1=,1?ρV1?e?εVn1+e?εVn,P1其中，n 是循环次数，参数 ρV 控制最大退化，εV 控制每次循环的退化梯度。当穿透深度超过 z1 时，循环次数 n 重置。

z1 和 zcn 之间的距离使用以下方程计算：(14)zcn?z1=αV(∑i=1n|Δzup|)βV其中，αV 和 βV 是控制 zcn?z1 的参数，Δzup 是增量挠度。

C-D 曲线：(15)P=Pd(2zr3.5zr6+1)zr=z?z3z1?z3D-E 曲线：(16)P=Pb((1?Pr)zr0.5+Pr)zr=z?z1zc?z1,Pr=Pd/Pb(z1)最大上拔吸力 [21]：(17)P2P1=φuplift(P2P1)deep=[1?(1?min(1,z^z^uplift,deep))Auplift]1/Buplift(P2P1)deep(18)z2=z1?(1+ωV)P1k01+φupliftωV?φuplift(z2?z3)=ψV(z1?z2)其中，φuplift 是控制上拔力的参数，(P2/P1)deep 是极限值，z^uplift,deep 是发生 (P2/P1)deep 时的归一化深度且等于 4，Auplift 和 Buplift 分别是拟合参数，等于 1.16 和 3.35，ωV 控制双曲线的渐近线，ψV 是控制路径 A-B 和 B-C 比例的参数，k0 是双曲线的初始斜率。

A-B 曲线：(19)P=P1+z?z11k0?z?z1(1+ωV)P1B-C 曲线：(20)P=P2(2zr2zr5+1)zr=z?z3z2?z3当从回滞环的边界或内部反转时，使用以下曲线：(21)Pre=Prev+z?zrev1k0+χVz?zrev(1+ω)P当穿透时 χV=1；当上拔时 χV=?1。

2.3. 动态管道沉降建模

基于离心机试验和流动法则理论，Cheuk [49] 提出了一个动态管道沉降模型。Cheuk [49] 将横向管道振荡犁出的土壤厚度称为“动态管道嵌入”。在本研究中，我们统一采用由横向管道振荡引起的动态管道沉降，Δz。基于线性化流动法则，Cheuk [49] 提出了一个包含穿透深度项的拟合公式（公式 (22)）。随着深度增加，在恒定 V/Pb 下，线性化流动法则的值逐渐减小，从而促进向最终沉降深度的收敛，如图 4 所示。(22)ΔzDout=pL(dzdu)f,L(ΔucycDout)qL下载：下载高分辨率图像 (197KB)下载：下载完整尺寸图像图 4. 随管道穿透变化的线性化流动法则 [49]。（在参考研究中，Vmax 对应本研究中的 Pb，w/D 对应 z/D）。

其中，pL 和 qL 是拟合参数，(dz/du)f,L 是线性化流动法则，Δucyc 是横向管道振荡幅度。(23)(dzdu)f,L=λL[VPb?(VPb)pp0?αL(zDout)]其中，(V/Pb)pp0 是代表海床表面值的常数，λL 和 αL 是拟合参数。V 是垂直载荷，Pb 是基于骨架曲线计算的土壤阻力，与 2.2 节中的定义相同。

然而，该模型存在两个主要问题。首先，初始循环中预测的沉降深度显著小于实验观察值。这种差异产生是因为，在 V/Pb 较高的初始循环期间，线性化流动法则产生的值远低于实际流动法则的值。其次，对于较大的初始穿透，预测结果与 CEL 模拟偏差很大。这是因为，在恒定 V/Pb 下，实际流动法则随深度增加而增加——这与包含穿透深度项的线性化流动法则相反（图 4）。因此，线性化流动法则在高和低 V/Pb 下均产生较低的预测值。因此，本研究从 [21] 的屈服面公式推导出实际流动法则，以取代线性化公式。

此外，该模型仅适用于恒定垂直载荷下的管道沉降。在管道涡激振动 (VIV) 和更一般的场景中，垂直载荷呈周期性变化，需要额外的参数来表征动态垂直载荷。在本研究中，引入了图 5 所示的参数来描述横向运动期间的动态载荷：Vmax 表示横向扫动内的最大垂直载荷，Vmean 表示横向扫动内的平均垂直载荷，uL 表示横向扫动上的横向位移距离，Δu 表示垂直载荷扫动面积质心到 uL/2 的距离，始终取正值。下载：下载高分辨率图像 (92KB)下载：下载完整尺寸图像图 5. 横向扫动期间垂直载荷的变化。

严格来说，在横向振荡期间的周期性垂直载荷下，沉降深度应在单次扫动内随横向位置变化，最大垂直载荷位置附近的沉降更大。然而，CEL 模拟表明，尽管存在周期性垂直载荷，穿透深度随横向位置的变化仅可忽略不计，甚至可能减小到零。因此，假设沉降深度在每次扫动内保持不变，并取最大值。

因此，单次横向扫动内的动态管道沉降增量 Δz 修改如下：(24)ΔzDout=pLexp(?|χLΔu|δL)(εLVmeanVmax)σL(dzdu)f(uLDout)qL其中，pL 和 qL 是与 Cheuk 模型中使用的类似拟合参数，uL 是单次扫动的横向位移，(dzdu)f 是流动法则。参数项 (?|χLΔu|δL)(εLVmeanVmax)σL 有效地量化了不同垂直载荷轮廓的影响，并在随后的参数拟合过程中显示出高精度，其中 χL、δL、εL 和 σL 是拟合参数。最大垂直载荷 Vmax 决定单次横向扫动期间的最大穿透深度。Vmean 和 Δu 用于考虑垂直载荷的平均大小和分布的影响。

因此，参数项以以上形式构建。由于流动法则与屈服面之间存在轻微的非关联性，需要对屈服面函数中的系数进行修改。屈服面可以写为[21]：
(25) f=H/h0 Pb?βfac(V/Pb?Vt/Pb)β1(1?V/Pb)β2=0
其中，H 为侧向土阻力，h0=Hmax/Pb 定义了峰值侧向阻力与垂直阻力的比值，Hmax 表示峰值侧向承载能力。Vt/Pb 是峰值垂直吸力与阻力的比值，随贯入深度变化。Vt/Pb 的计算公式与 (17) 相同，仅将 P2 替换为 Vt，P1 替换为 Pb。βfac、β1 和 β2 为参数。β1 和 β2 均等于 0.65。
本研究推导流动法则 (dz/du)f 如下：
(26) (dz/du)f=?f/?V / ?f/?H = ?h0(β1βfac(V/Pb?Vt/Pb)β1?1(1?V/Pb)β2?β2βfac(V/Pb?Vt/Pb)β1(1?V/Pb)β2?1)
h0 和 βfac 的计算公式如下：
(27) h0=h0,surface+?h(h0,deep?h0,surface)
βfac=(β1+β2)(β1+β2)β1β1β2β2(1?Vt/Pb)(β1+β2)
其中，h0,surface 是海床表面处的 h0 值，h0,deep 是 h0 的极限值，描述了在几个直径贯入深度下的水平承载能力（此时 h0 与贯入深度无关），?h 是一个转换因子，定义为：
(28) ?h=[1?(1?min(1,z^/z^h,deep))Ah]1/Bh
其中，z^h,deep 是发生 h0,deep 时的归一化深度，等于 3.5，Ah 和 Bh 是拟合参数。
在开发线性化流动法则时，Cheuk [49] 引入了一个与深度相关的偏移项。本研究随后的参数拟合表明，该偏移项显著提高了模型精度。因此，当前的流动法则公式结合了类似的基于贯入深度的偏移，将 V/Pb 替换为 (V/Pb)pp 如下：
(29) (V/Pb)pp=V/Pb?(αLz/Dout+(V/Pb)pp0)
其中，Pb 基于骨架曲线计算，但将剪切强度 su 替换为 su,op，后者包含了土体强度退化。
(30) su,op/su=δrem+(1?δrem)exp[?3ζrL((uL/Dout)0+∑((uL/Dout)+(zmax,L?zmin,L)/Dout))∑(uL/Dout)95]
ζ=min[1,(z/Dout)av/(z/Dout)]
其中，δrem=1/St，∑(uL/Dout)95=3N95，∑(uL/Dout) 表示之前的扰动，(uL/Dout)0 是表示初始扰动的常数。
方程 (30) 中引入了一个额外因子 ζrL，以考虑管道下沉时贯入“新鲜”土体的影响。因子 rL 是一个拟合参数。ζ 由之前扰动期间的平均贯入深度 (z/Dout)av 与当前贯入深度之比控制。虽然 Cheuk [49] 仅考虑了由侧向振荡引起的退化效应，但本研究通过引入一个额外项 (zmax,L?zmin,L)/Dout 进一步考虑了由垂直循环运动引起的退化，其中 zmax,L 是一次侧向扫描内的最大贯入深度，zmin,L 是最小深度。
动态管道沉降模型与垂直 PSI 模型之间的耦合方案遵循以下原则：当一次侧向扫描引起沉降增量 Δz 后，当维持或重新获得相同的土阻力时，从垂直 PSI 模型计算出的管道贯入深度 z 以及其他荷载必须增加 Δz，如图 6 所示。因此，实际贯入深度 z′ 满足 z′=z+Δz。当发生沉降增量 Δz 后，垂直 PSI 模型中使用的深度 z 更新为 z=z′?Δz。
下载：下载高分辨率图像 (699KB)
下载：下载全尺寸图像
图 6. 动态管道沉降模型与垂直非线性 PSI 模型之间的耦合方案及流程图。
当侧向扫描完成后，管道可能呈现两种运动状态之一：(1) 贯入/平衡状态，和 (2) 上拔状态。由于贯入状态和平衡状态表现出相似的行为，它们被归为一类。当管道处于贯入/平衡状态时，假定沉降增量 Δz 在侧向扫描后发生。因此，土反力点瞬间从 F 移至 F′。为了确保数值收敛和响应的连续性，沉降增量 Δz 在 n 个时间步长内施加，每步通过 z1=z′1?Δz/n 使贯入深度减少 Δz/n，实现转换 F-F′（图 6(a) 中的红色曲线）。
在上拔状态期间，沉降增量 Δz 影响随后的重新贯入阶段。当重新获得相同的土阻力时，土反力点 D 向右移动 Δz+Δz′（z′1=z1+Δz+Δz′），其中 Δz′ 代表由土体强度退化和垂直荷载引起的贯入增量，由垂直 PSI 模型计算。为了保持数值过程的一致性，采用相同的方法，在重新贯入阶段每步将贯入深度减少 Δz/n。步数 n 可根据振动频率、积分时间步长和收敛要求进行调整。本研究建议 1000–3000 个时间步长。
计算过程如图 6 (b) 所示。流程图的左侧显示了 VFIFE 计算过程。给定初始节点坐标和荷载，通过中心差分格式获得下一个时间步长的节点坐标。然后根据更新的节点坐标计算作用于节点上的内力和等效外荷载。更新节点力后，进行下一个积分步，重复此过程直到达到规定的总持续时间。流程图的右侧显示了动态沉降计算过程。由于动态沉降是在 n 个时间步长内逐渐施加的，因此首先在每个积分步之后确定由之前的侧向扫描引起的任何剩余动态沉降是否仍需要施加。如果已完全施加，垂直 PSI 计算中使用的贯入深度 z 直接从积分结果 zn+1 获得。否则，需要进行额外检查以确定当前的侧向扫描是否已完成。如果由之前的侧向扫描引起的动态沉降尚未完全施加，而当前的侧向扫描已完成，则减少时间步长数 n，然后应重新开始计算过程。如果当前的侧向扫描尚未完成，则垂直 PSI 计算中使用的贯入深度 z 减少 Δz/n。然后，计算出的垂直土阻力作为外荷载并入 VFIFE 计算中。
目前，对于侧向管道振荡引起的动态沉降后重新贯入曲线的后段（F’-G 或 D’-E’）仍缺乏公式。由于犁耕效应，垂直 PSI 模型中的方程 (14) 显然不适用于计算距离 z′c?z′1。侧向振荡后，最大贯入深度下方的土体强度同时受到垂直和水平扰动的影响而降低，导致 z′c 逐渐增大。此外，由于犁耕效应，z′1 和 z′c 同时增加。
如前所述，CEL 技术在 PSI 模拟中提供了更高的精度。因此，采用基于 ABAQUS 的 CEL 模型来研究重新贯入曲线的后段（F’-A 或 D’-E’）。分析结果表明，方程 D-E（式 (16)）能够令人满意地拟合这些重新贯入曲线；因此，保留其公式。距离 z′c?z′1 的方程给出在 (31) 中，详细的验证和参数拟合在第 3.3 节中呈现。
(31) z′C?z′1=ζVz′1(P′1/Pb)τV/Dout
其中，ζV 和 τV 是参数。
2.4. 水动力与尾流振荡器方程
Qu [33] 提出了一种尾流振荡器模型，用于预测两自由度 (2-DOF) 弹性安装刚性圆柱体的涡激振动 (VIV)。Fehér [34] 在 [33] 的尾流振荡器方程中引入了三个额外的经验系数，以考虑平面边界的影响。[33] 提出并在 [34] 修改后的尾流振荡器模型给出如下：
(32) q¨+εVIVωs(q2?1)q˙+ωs2q+ψVIV(e/Dout)ωnωsq?ηVIV(e/Dout)κVIVDouty¨q=βVIV(e/Dout)AVIVDoutz¨
其中，q 是尾流变量。ωs 是旋涡脱落频率，定义为 ωs=2πStVflow/Dout，St 为斯特劳哈尔数，Vflow 为均匀流速度。εVIV、κVIV 和 AVIV 是拟合实验数据的调谐参数。ωn 是圆柱体的固有振动频率。y¨ 和 z¨ 分别是圆柱体在 IL 和 CF 方向的加速度。ψVIV(e/Dout)、ηVIV(e/Dout) 和 βVIV(e/Dout) 是考虑平面边界影响的系数，e/Dout 是间隙比，e 是平面边界与圆柱体之间的距离。
尾流变量 q 通过中心差分格式求解，连续和离散的中心差分公式表达如下：
(33) q1=q0+Δtq˙0+q¨0Δt2/2
qn+1=(2?Δt2(ωs2+ωnωsψVIV?ηVIVκVIVy¨/Dout))qn+(0.5εVIVωsΔt(qn2?1)?1)qn?1+AVIVΔt2βVIVz¨/Dout / (1+0.5εVIVωsΔt(qn2?1))
文献 [34] 中用于确定经验系数 ψVIV、ηVIV 和 βVIV 的方程给出如下：
(34) {βVIV(e/Dout)=0.0083(e/Dout)3?0.1(e/Dout)2+0.3917(e/Dout)+0.5
ηVIV(e/Dout)=0.0083(e/Dout)3?0.1(e/Dout)2+0.3917(e/Dout)+0.5
ψVIV(e/Dout)=0, 2.0≤(e/Dout)≤5.0
{βVIV(e/Dout)=?0.2667(e/Dout)2+0.8667(e/Dout)+0.3
ηVIV(e/Dout)=?0.667(e/Dout)2+1.3667(e/Dout)+0.1
ψVIV(e/Dout)=0.5, 0.75≤(e/Dout)≤1.5
{βVIV(e/Dout)=?6.8(e/Dout)2+2.3(e/Dout)+1.45
ηVIV(e/Dout)=0.5
ψVIV(e/Dout)=?4.4(e/Dout)2+1.3(e/Dout)+0.45, 0.0≤(e/Dout)≤0.5
根据 [33] 中的方程，作用于单位长度圆柱体的瞬时阻力 FVD 和瞬时升力 FVL 可以写为：
(35) FVD=1/2 CDM ρwater Dout Un|Un|+1/2 αVIV CVL2 ρwater Dout Ux|Ux|
FVL=1/2 CVL ρwater Dout |Un|2 eaFLN
Ux=(Un·Vn|Vn|)·Vn|Vn|
CVL=qCL0/2
其中，CDM 是振荡平均阻力系数，CVL 是涡升力系数，CL0 是在静止圆柱体上测量的升力系数，ρwater 是流体密度，αVIV 是尾流振荡器系数，Vflow 是均匀流速度，U 是相对于圆柱体运动的流速度，U=Vflow?y˙，Vn 是 Vflow 垂直于圆柱体的分量，y˙ 是圆柱体在 IL 方向的速度，Un 是 U 垂直于圆柱体的分量，Ux 是 Un 在 Vn 方向的分量。
在 Fortran 程序中，Un 使用 ea1 和 U 的叉积计算（图 7）：
(36) |Un|=|ea1×U|
Un=ea1×U×ea1
下载：下载高分辨率图像 (69KB)
下载：下载全尺寸图像
图 7. Un 示意图。
eaFLN 计算如下（图 8）：
(37) eaFLN=ea1×Un|ea1×Un|
下载：下载高分辨率图像 (75KB)
下载：下载全尺寸图像
图 8. eaFLN 示意图。
2.5. 侧向非线性管土相互作用建模
本研究采用 [43] 中的侧向非线性 PSI 模型，该模型采用三线性公式，如图 9 所示。
下载：下载高分辨率图像 (67KB)
下载：下载全尺寸图像
图 9. 侧向 PSI 三线性模型。
不排水侧向突破阻力模型如下：
(38) PL,bre,u/suDout=1.7(z/Dout)0.6+0.23(V/suDout)0.8+0.6(γsoilDout/su)(z/Dout)2
其中，PL,bre,u 是不排水土条件下列向突破阻力。
不排水侧向残余阻力 PL,res,u 可计算为：
(39) PL,res,u/V=0.32+0.8(z/Dout)0.8
突破阻力的侧向动员距离 ybrk 估计为：
(40) ybrk/Dout=0.02+0.25z/Dout
残余阻力的侧向动员位移 yres 估计为：
(41) yres/Dout=1.5
3. 模型验证与参数拟合
3.1. CEL PSI 模型
在 ABAQUS 软件中使用 CEL 技术建立了一个 PSI 有限元模型。实现了用户子程序以考虑应变率和应变软化对土体剪切强度的影响。该模型与 [36,50] 中报告的模拟数据以及 [49] 中报告的实验数据进行了验证。首先，使用 [36] 中的案例 L00、LS0、L0R 和 LSR 验证土阻力，相应参数列于表 1。案例 L00、LS0、L0R 和 LSR 分别对应以下情况：忽略应变率和应变软化效应，仅考虑应变软化效应，仅考虑应变率效应，以及同时考虑两者效应。如图 10 所示，数值结果与参考数据吻合良好。因此，开发的 CEL PSI 模型准确捕捉了考虑应变率和应变软化效应的土阻力。
表 1. 有限元分析中使用的几何和土体参数。
空单元格 | 参数 | 案例 A | Dutta, 2015 [36] | 案例 B | 本研究
管道 | 外径，Dout | 0.8 m | 0.4572m
土体 | 泥线处土体不排水剪切强度，su0 | 2.3kPa | 1.2kPa
不排水剪切强度梯度，suk | 3.6kPa/m | 3.6kPa/m
不排水弹性模量，Eu | 500su | 500su
泊松比，vu | 0.495 | 0.495
有效单位重量，γsoil | 6.5kN/m3 | 6.5kN/m3
剪切强度增加率，μ | 0.1 | 0.1
参考剪切应变率，γ˙res | 3×10?6 | 3×10?6
重塑土敏感性，St | 3.2 | 3.2
土体强度退化 95% 时的累积绝对塑性剪切应变，ξ95 | 10 | 10
下载：下载高分辨率图像 (275KB)
下载：下载全尺寸图像
图 10. 垂直土阻力结果与文献数据的对比。
随后，验证了 CEL PSI 模型预测由侧向管道振荡引起的动态管道沉降的精度。验证使用了 [49] 中报告的案例 KC05 和 HP07 的第 1 阶段数据，参数列于表 2。如图 11 所示，40 个循环内计算的动态管道沉降与参考数据吻合良好，证实了该模型准确捕捉了侧向振荡期间的沉降。
表 2. [49] 中使用的几何和土体参数。
空单元格 | 参数 | KC05 | HP07
管道 | 外径，Dout | 0.8 m | 0.8 m
侧向振荡幅度，Δucyc | 0.05Dout | 0.05Dout
初始管道沉降，zinitial | 0.12Dout | 0.22Dout
垂直荷载，V | 2.17kN/m | 2.52kN/m
土体 | 泥线处土体不排水剪切强度，su0 | 0.75kPa | 0.40kPa
不排水剪切强度梯度，suk | 1.6kPa/m | 2.5kPa/m
不排水弹性模量，Eu | 500su | 500su
泊松比，vu | 0.495 | 0.495
有效单位重量，γsoil | 6.0kN/m3 | 3.0kN/m3
剪切强度增加率，μ | 0.1 | 0.1
参考剪切应变率，γ˙res | 3×10?6 | 3×10?6
重塑土敏感性，St | 4.0 | 1.7
土体强度退化 95% 时的累积绝对塑性剪切应变，ξ95 | 10 | 10
下载：下载高分辨率图像 (227KB)
下载：下载全尺寸图像
图 11. 动态沉降结果与文献数据的对比。
3.2. 基于 VFIFE 方法和尾流振荡器模型的管道模型
基于 VFIFE 方法开发了一个管道动力学模型，并集成了尾流振荡器模型。验证使用了 [34] 中报告的数据，参数列于表 3。对比结果展示在图 12 中。当折算速度低于 8 或高于 8.5 时，CF 归一化幅值和 IL 平均阻力系数均观察到了极好的一致性。仅在 8—8.5 折算速度范围内观察到轻微偏差。这些结果证实了该模型在预测 IL 和 CF 振动方面的准确性。
表 3.

[34] 中圆柱体和尾流振荡器模型的参数。

| 参数 | 值 |
| :--- | :--- |
| **圆柱体** | |
| 外径，Dout | 0.04 m |
| 阻尼比，ζdamp | 0.02 |
| 圆柱体固有频率，ωn,cylinder | 8.17 rad/s |
| 质量比，mcylinder* | 5.5 |
| **尾流振荡器模型** | |
| 振荡平均阻力系数，CDM | 1.1 |
| 静止圆柱体上测得的升力系数，CL0 | 0.3 |
| 生成 1.1 振荡平均阻力系数的系数，αVIV | 2.2 |
| 附加质量系数，Ca | 1 |
| 水的密度，ρwater | 1000 |
| 斯特劳哈尔数，St | 0.2 |
| AVIV | 16 |
| εVIV | 0.1 |
| κVIV | 3 |

下载：下载高分辨率图片 (333KB)
下载：下载全尺寸图片
图 12. VIV 模拟结果与 [34] 中报告的参考数据的比较。

为了进一步验证跨空管道 VIV 模拟的模型准确性，使用 [16] 中报告的实验数据进行了额外的验证，参数列于表 4。如图 13 所示，CF 归一化振幅比较表明该模型能准确模拟跨空管道的 VIV 响应。振幅峰值处观察到的轻微偏差可能归因于实验中使用的张重块的惯性效应，这影响了轴向张力，进而影响了管道的固有频率。

表 4. [16] 中管道和尾流振荡器模型的参数。

| 参数 | 值 |
| :--- | :--- |
| **管道** | |
| 外径，Dout | 0.016 m |
| 壁厚，tpipe | 0.003 m |
| 长度，L | 2.6 m |
| 弹性模量，E | 210 GPa |
| 阻尼比，ζ | 0.035 |
| 圆柱体固有频率，ωn,cylinder | 16.78 rad/s |
| 质量比，mcylinder* | 2.62 |
| **尾流振荡器模型** | |
| 振荡平均阻力系数，CDM | 1.1 |
| 静止圆柱体上测得的升力系数，CL0 | 0.3 |
| 生成 1.1 振荡平均阻力系数的系数，αVIV | 2.2 |
| 附加质量系数，Ca | 1 |
| 水的密度，ρwater | 1000 |
| 斯特劳哈尔数，St | 0.2 |
| AVIV | 16 |
| εVIV | 0.1 |
| κVIV | 3 |

下载：下载高分辨率图片 (154KB)
下载：下载全尺寸图片
图 13. VIV 模拟结果与 [16] 中报告的实验数据的比较。

### 3.3. 动态管道沉降模型的参数拟合

垂直 PSI、横向 PSI 和动态沉降模型被纳入基于 VFIFE 的管道动力学模型中，随后进行了参数校准。虽然文献 [21] 中报告了 (27)、(28) 中的参数 h0,surface、Ah 和 Bh 的参考值，但文献中使用的骨架曲线忽略了土壤浮力效应。由于本研究在骨架曲线中引入了浮力项，这些参数需要重新校准。尽管 Hodder 的骨架曲线没有明确考虑浮力——且由于管道直径较小，其影响有限——但在测试过程中管道不可避免地受到了浮力作用。因此，实验数据和屈服面方程仍具有重要的参考价值。Merifield [59] 将这些解进行了扩展以考虑土壤浮力效应，提供了 0.5Dout 范围内 H 和 V 的预测曲线。

图 14 展示了 Hodder 公式、Merifield 解和 ABAQUS 模拟的比较结果。模型参数与文献 [21] 保持一致，Dout=0.5m，su0=3.5kPa，suk=0.7kPa/m，而 CEL PSI 模型的其他参数与表 1 (Case B) 一致。结果显示，Hodder 方程在 0.3-0.4Dout 范围内略微低估了数值，而在 0-0.2Dout 范围内略微高估了数值。重新校准的参数 (h0,surface=0.024, Ah=1.049 和 Bh=2.505) 与参考数据具有更好的吻合度。

下载：下载高分辨率图片 (365KB)
下载：下载全尺寸图片
图 14. h0 与参考数据的比较。

ζV 和 τV 的参数拟合程序如下进行。建立了初始贯入深度分别为 0.12、0.24、0.50 和 0.75Dout 的基准场景。在振幅为 0.05Dout 和 0.15Dout 的条件下进行 40、50 和 60 次横向循环后，实施了位移控制的贯入。在横向振荡期间施加恒定的垂直荷载，四种初始贯入深度条件对应的垂直荷载分别为 2.05、2.80、4.25 和 5.75kPa。模型参数列于表 1 (Case B)。图 15 展示了骨架曲线、再贯入土阻力以及使用公式 D-E (式 (16)) 拟合的曲线（针对部分情况）。公式 (16) 和 (31) 的结合应用在 ζV=0.009 和 τV=1.209 的参数配置下，产生了精确的土阻力曲线。

下载：下载高分辨率图片 (464KB)
下载：下载全尺寸图片
图 15. 再贯入土阻力与 CEL 结果的比较。

图 16 展示了参数 ζV=0.009 和 τV=1.209 的式 (31) 响应面，以及从所有情况 CEL 模拟中获得的 z′c?z′1 值。

下载：下载高分辨率图片 (203KB)
下载：下载全尺寸图片
图 16. 式 (31) 的响应面。

同样，使用 CEL 技术在 ABAQUS 中建立了一个 CEL PSI 模型，采用显式动态分析。进行了模拟以研究不同垂直荷载施加方法下的动态管道沉降。模型配置和网格尺寸如图 17 所示。文献 [36] 中的最小单元尺寸为 0.05Dout (0.04 m)，而本研究基于文献数据和网格敏感性分析将其设置为 0.044Dout (0.02 m)。模型参数列于表 1 (Case B)。

下载：下载高分辨率图片 (394KB)
下载：下载全尺寸图片
图 17. CEL PSI 模型的模型配置和网格。

图 18 显示了五种不同的垂直荷载施加曲线，分别指定为 N、V、L、R 和 U。

下载：下载高分辨率图片 (199KB)
下载：下载全尺寸图片
图 18. 垂直荷载施加曲线。

总共使用了 60 个模拟案例进行参数拟合，相应的参数条件见表 5。

表 5. CEL 模拟的参数条件。

| 初始管道沉降，zinitial | 最大垂直荷载，Vmax | 横向振荡振幅，Δucyc | 垂直荷载施加曲线 |
| :--- | :--- | :--- | :--- |
| 0.12Dout | 2.05 kPa | 0.05Dout, 0.10Dout, 0.15Dout | N, V, L, R, U |
| 0.24Dout | 2.80 kPa | | |
| 0.50Dout | 4.25 kPa | | |
| 0.75Dout | 5.75 kPa | | |

鉴于参数范围，采用粒子群优化 (PSO) 来确定最佳参数组合。使用本研究开发的模型进行了相同的模拟。目标函数定义为两个模型在所有情况下的每次横向振荡循环后管道沉降的均方误差 (MSE)，如以下方程所示。

(42) MSE = 1 / (ncase * ncycle) * ∑(i=1 to ncase) ∑(j=1 to ncycle) (zCEL - z?)^2

其中，ncase 是案例数量，ncycle 是每个案例的横向振荡循环次数，zCEL 表示 CEL 模拟数据，z? 表示所提模型的预测值。

优化结果和参数组合如图 19 和表 6 所示。此外，图 19 还显示了使用 Cheuk 模型对 12 个恒定垂直荷载案例获得的模拟结果。对于初始贯入深度为 0.12D 和 0.24D 的案例，Cheuk 模型在前 20 个周期内的拟合性能劣于所提模型。对于初始贯入深度为 0.5D 和 0.75D 的案例，Cheuk 模型在后 20 个周期内的拟合性能也较差。对于这 12 个案例，两个模型的 MSE 值分别为 0.0070 和 0.0012。

下载：下载高分辨率图片 (1MB)
下载：下载全尺寸图片
下载：下载高分辨率图片 (933KB)
下载：下载全尺寸图片
图 19. CEL 模拟与当前模型获得的动态沉降结果比较 (红色曲线代表当前模型，黑色虚线曲线代表 Cheuk 模型，蓝色曲线代表 CEL)。

表 6. 动态管道沉降模型的参数组合。

| 参数 | 值 | 注释 |
| :--- | :--- | :--- |
| St | 2.5 | 目标土壤参数 |
| N95 | 4rL | |
| 2.62 | | 校准参数 |
| pL | 0.83 | |
| qL | 1.01 | |
| (uL/Dout)0 | 1.90 | |
| αL | 0.10 | |
| (V/Pb)pp0 | -0.03 | |
| χL | 5.26 | |
| δL | 3.77 | |
| σL | 1.54 | |
| εL | 0.65 | |

需要注意的是，此参数组合专用于本研究；对于其他场景，此处描述的计算程序可用于根据模拟或实验数据确定特定的最佳参数。

本研究还调查了横向振荡速度对动态管道沉降的影响。模拟矩阵包含 16 个参数条件，结合了两种初始贯入深度 (0.12Dout, 0.75Dout)、两种横向振荡振幅 (0.05Dout, 0.15Dout) 和四种频率 (0.5Hz, 1Hz, 2Hz 和 3Hz)。50 个周期后，管道沉降深度的比较如图 20 所示，表明横向振荡引起的动态沉降深度对振荡速度不敏感。因此，所提动态沉降模型排除了振荡速度效应，并在所有 CEL 模拟中采用统一的 1Hz 频率。

下载：下载高分辨率图片 (174KB)
下载：下载全尺寸图片
图 20. 不同横向振荡频率下的管道沉降结果比较。

### 3.4. 动态沉降模型的讨论

为了进一步验证动态沉降模型的准确性，模拟了文献 [49] 中报告的四个案例，即 KC04、KC05、HP06 和 HP07。所提模型中使用的参数见表 7 中的组合 1。如图 21 所示，所提模型能准确模拟管道的动态沉降，且其预测精度显著高于 Cheuk 模型。初始阶段的预测和最终沉降值均与实验数据非常吻合。此外，两组相同类型的土壤模拟使用了相同的参数集，这进一步证明了所提模型的通用性和适用性。

表 7. 所提模型中使用的参数。

| 参数 | 组合 1 | 组合 2 | 组合 3 |
| :--- | :--- | :--- | :--- |
| rL | 5.0 | 2.86 | 1.0 |
| pL | 1.0 | 1.0 | 1.0 |
| qL | 1.26 | 1.20 | 1.25 |
| (uL/Dout)0 | 5.0 | 2.29 | 5.0 |
| αL | 0.85 | 1.14 | 0.37 |
| (V/Pb)pp0 | -0.18 | -0.26 | 0.50 |
| χL | 5.26 | 5.26 | 5.26 |
| δL | 3.77 | 3.77 | 3.77 |
| σL | 1.54 | 1.54 | 1.54 |
| εL | 0.65 | 0.65 | 0.65 |

下载：下载高分辨率图片 (466KB)
下载：下载全尺寸图片
图 21. 针对 Cheuk 研究中报告的案例的两个模型比较 (红色曲线代表当前模型，蓝色虚线曲线代表 Cheuk 模型，黑色虚线曲线代表实验结果)。

如果针对每种类型的土壤单独进行参数优化，可以实现对实验数据的更好拟合 (图 22)。表 7 中的参数组合 2 用于案例 KC04 和 KC05，而参数组合 3 用于案例 HP06 和 HP07。

下载：下载高分辨率图片 (442KB)
下载：下载全尺寸图片
图 22. 使用参数组合 2 和 3 的比较 (红色曲线代表当前模型，蓝色虚线曲线代表 Cheuk 模型，黑色虚线曲线代表实验结果)。

模型参数的适用性讨论如下。St 和 N95 是固有的土壤性质。经过广泛的计算比较后，我们发现直接使用其实际值也能实现高预测精度。因此，在本节中，使用了 St 和 N95 的实际值。

需要拟合的参数包括与沉降深度相关的 rL、pL 和 qL；与流动法则相关的 (uL/Dout)0、αL 和 (V/Pb)pp0；以及与垂直荷载相关的 χL、δL、σL 和 εL。使用的管道和土壤参数组合列于表 8，其余参数与表 1 中的 Case B 相同。对于每种管道 - 土壤组合，计算了 60 个具有不同条件的动态贯入案例。案例设置与表 5 相同，除了 Vmax 的值。Vmax 设置为管道达到初始贯入深度时作用在管道上的土阻力。

表 8. 土壤和管道参数。

| 组合 | 参数 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| **管道** | 外径，Dout | 0.4572 m | 0.4572 m | 0.4572 m | 0.4572 m | 0.4572 m | 0.60 m | 0.80 m |
| **土壤** | 泥面处土壤的不排水抗剪强度，su0 | 0.6 kPa | 1.2 kPa | 1.2 kPa | 1.8 kPa | 2.4 kPa | 1.2 kPa | 1.2 kPa |
| | 不排水抗剪强度梯度，suk | 2.4 kPa/m | 2.4 kPa/m | 3.6 kPa/m | 3.6 kPa/m | 4.8 kPa/m | 3.6 kPa/m | 3.6 kPa/m |

在之前的分析中，与垂直荷载相关的参数 χL、δL、σL 和 εL 保持不变，同时获得了合理的准确预测。在本节中，我们发现这四个参数的变化范围非常小。因此，我们认为与垂直荷载相关的参数具有相对较广的适用性，并且这四个参数保留与之前使用的相同值。

表 9 列出了用于参数拟合的案例组合以及这些案例获得的相应 MSE。当使用所有七个管道 - 土壤组合进行参数优化时，最佳参数集产生的 MSE 为 0.0032。当使用案例 1-5（具有相同管道直径但不同土壤参数）进行参数优化时，最佳参数集产生的 MSE 为 0.0028。当使用案例 3、6 和 7（具有相同土壤参数但不同管道直径）进行参数优化时，最佳参数集产生的 MSE 为 0.0019。因此，与管道直径相比，土壤参数对拟合参数的适用性有更为显著的影响。

表 9. 参数集和 MSE。

| 案例 | 参数 | 案例 1～7 | 案例 1～5 | 案例 3,6,7 | 案例 1 | 案例 2 | 案例 3 | 案例 4 | 案例 5 | 案例 6 | 案例 7 |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| | rL | 5.45 | 5.86 | 4.34 | 1.00 | 1.54 | 3.03 | 4.43 | 6.68 | 1.74 | 1.00 |
| | pL | 1.00 | 1.14 | 0.92 | 0.49 | 0.60 | 0.64 | 0.73 | 0.84 | 0.64 | 0.533 |
| | qL | 0.94 | 0.93 | 0.85 | 0.76 | 0.86 | 0.92 | 0.99 | 1.04 | 0.91 | 0.83 |
| | (uL/Dout)0 | 2.32 | 0.10 | 1.32 | 5.10 | 2.83 | 2.65 | 1.95 | 2.19 | 3.50 | 4.54 |
| | αL | 0.025 | 0.030 | 0.021 | 0.010 | 0.024 | 0.059 | 0.062 | 0.091 | 0.037 | 0.01 |
| | (V/Pb)pp0 | 0.053 | 0.042 | 0.090 | 0.275 | 0.191 | 0.049 | -0.006 | -0.065 | 0.175 | 0.266 |
| | χL | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 | 5.26 |
| | δL | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 | 3.77 |
| | σL | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 | 1.54 |
| | εL | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 |
| | MSE | 0.0032 | 0.0028 | 0.0019 | 0.0008 | 0.00066 | 0.00041 | 0.00026 | 0.00018 | 0.00087 | 0.0017 |

还对案例 1-7 中的每个管道 - 土壤参数组合分别进行了参数拟合。当针对每个管道 - 土壤参数组合单独进行参数优化时，所提模型的精度可以显著提高。因此，除与垂直荷载相关的参数外，其他参数具有一定的适用性。然而，针对特定案例的参数优化当然能产生更准确的结果。

### 3.5. 垂直非线性管土相互作用模型

使用文献 [21] 中的实验数据验证了垂直 PSI 模型。关键参数列于表 10。图 23 中展示的比较结果表明，计算预测与实验数据在三个关键阶段（初始贯入、上抬和再贯入）均具有极好的吻合度。采用循环贯入轨迹来专门检查循环退化行为。如图 23(b) 所示，该模型准确捕捉了循环荷载下的土强度退化。

表 10.

竖向非线性 PSI 模型的模型参数

| 类别 | 参数 | 数值 |
| :--- | :--- | :--- |
| 管道 | 外径，Dout | 0.5 m (原型尺度) |
| 土壤 | 泥面处土壤不排水抗剪强度，su0 | 3.5 kPa |
| | 不排水抗剪强度梯度，suk | 0.7 kPa/m |
| | 水下重度，γsoil | 6kN/m3 |
| 土壤模型 | a | 6.15 |
| | b | 0.15 |
| | ω | 0.9 |
| | ψ | 3 |
| | k0 | 1 |
| | | 200kPa |
| | αV | 1.0 |
| | βV | 0.25 |
| | ρV | 0.4 |
| | εV | 0.8 |

图 23. 竖向非线性 PSI 模型的验证结果。

4. 结果与讨论

4.1. 模拟工况

图 24 描绘了悬跨管道的配置，突出了其关键参数。L3 表示单跨长度。L1 和 L2 表示两侧肩部段的长度，这些段承受横向和竖向 PSI。两个端点均承受轴向拉力 T，且所有三个旋转自由度和纵向 (IL) 位移均完全约束。e 表示跨段与海床之间的间隙。

图 24. 悬跨管道示意图。

表 11 中列出的管道参数指的是南海的一个海上石油项目。竖向 PSI 模型的参数与表 10 一致，而尾流振荡器模型的参数与表 3 一致。

表 11. 管道和土壤参数

| 类别 | 参数 | 数值 |
| :--- | :--- | :--- |
| 管道 | 海水密度，ρwater | 1025kg·m?3 |
| | 外径，Dout | 0.4572m |
| | 内径，Din | 0.381m |
| | 弹性模量，E | 210GPa |
| | 剪切模量，G | 80.7GPa |
| | 管道密度，ρpipe | 7850kg·m?3 |
| | 原油密度，ρoil | 850kg·m?3 |
| | 泊松比，μ | 0.3 |
| | 左肩长度，L1 | 170m |
| | 右肩长度，L2 | 170m |
| | 跨段长度，L3 | 60m |
| | 轴向拉力，T | 1300kN |
| | 悬跨管道在 CF 方向的固有频率，ωn,z(fn,z) | 3.33rad·s?1(0.53Hz) |
| | 悬跨管道在 IL 方向的固有频率，ωn,y(fn,y) | 2.59rad·s?1(0.41Hz) |
| | 阻尼比，μd | 0.05 |
| 土壤 | 泥面处土壤不排水抗剪强度，su0 | 1.2 kPa |
| | 不排水抗剪强度梯度，suk | 3.6 kPa/m |
| | 水下重度，γsoil | 6.5kN/m3 |

4.2. 悬跨管道的涡激振动响应

图 25(a) 展示了沿整个管道在 CF 方向上的最大均方根 (RMS) 折减振动振幅以及最大折减振动振幅。所有结果均为折减速度 Vr=πV/(ωn,z+ωn,y)Dout 的函数。由于两个方向的固有频率不同，因此采用两个固有频率的平均值来计算折减速度 Vr。当 Vr=4 时，涡脱落频率接近管道的固有频率，悬跨段开始表现出一阶振动模态。CF 振动频率稳定在 0.81fn,z-1.32fn,z 范围内。当 Vr=7 时，折减振幅达到最大值 0.56。随后，随着涡脱落频率与固有频率之间的偏差继续增加，振动进入过渡区，振幅逐渐减小。当 Vr>10 时，涡脱落频率接近固有频率的两倍，导致振动频率发生跳跃。CF 振动随后进入第二个锁定区，频率范围为 2.09fn,z 到 2.32fn,z，振幅逐渐增加。在研究范围内，最大振幅达到 0.45。RMS 振幅跟随最大振幅的趋势相似，但数值较低。

图 25. 最大均方根 (RMS) 折减振动振幅、最大折减振动振幅和主频率。

与 CF 振动相比，IL 振动振幅要小得多 (图 25(b))。当一阶或二阶模态被完全激发时，最大 IL 振幅仅达到相应 CF 振幅的 10–25%。在低折减速度 (Vr=2) 下，涡脱落频率接近固有频率的两倍。因此，波动阻力频率变得接近固有频率，从而激发小振幅的 IL 共振。IL 振动也表现出两个锁定区，频率分别稳定在 2.07fn,y-3.02fn,y 和 5.37fn,y-5.93fn,y 范围内。值得注意的是，IL 振动频率始终近似为 CF 振动频率的两倍。

图 26 展示了时空变化的折减位移云图，以及 RMS 和最大振幅的沿跨演化。纵坐标值表示折减位移，已去除管道初始位置的平均偏移量。仅显示中心 200m 段，因为剩余段的振动振幅可忽略不计。后续图表中的显示范围可根据横坐标值确定，横坐标值表示距离左端点的距离。

图 26. 时空变化折减位移云图、RMS 振动振幅沿跨演化及最大振动振幅。

在 CF 方向，当振动处于一阶模态 (Vr=7) 时，显著振幅仅发生在跨中。在肩部，峰值振幅仅达到 0.08，相当于跨中值的 10.71%。在 IL 方向，尽管振幅较小——峰值为 0.065——但肩部段表现出 0.030 的振幅，相当于跨中值的 46.15%。此外，跨中两侧出现了驻波和行波的混合特征，表明 IL 振动包含多个频率成分。

当 Vr=10 时，CF 振动仍以一阶模态为主，而在 IL 方向出现了模态竞争。IL 振幅的空间分布在“双峰”和“单峰”模式之间交替，逐渐以行波特征为主。肩部段出现了另外两个振幅峰值，肩跨与跨中振幅比增加到 58.14%。

当 Vr=11 时，CF 振动转变为二阶模态，表现出典型的驻波特征。肩部振幅显著增加到 0.11，相当于在 Vr=10 条件下观察到的值的 3.2 倍，并接近最大悬跨振幅的 44%。在 IL 方向，肩部段表现出比悬跨段更大的振幅 (高达 0.029)，并具有显著的行波特征。相比之下，跨中区域显示出相对稳定的位移波峰和波谷分布，表明具有局部驻波特征。IL 振动同时显示多个波峰和波谷，波峰的数量和空间分布随时间动态演变。

这表明 IL 方向存在多模态叠加和频率竞争，导致更宽的频带带宽和更复杂的频谱内容。频谱展宽可归因于多个因素。首先，在过渡区，一阶和二阶模态共存并竞争。其次，当 CF 振动振幅较大时，管道相对于来流的瞬时速度和攻角会周期性地变化。这可能会在 IL 响应中引入一个 CF 振动频率分量以及谐波分量。此外，随着 Vr 的增加，尾流剪切层的不稳定性增强，这也可能在 IL 方向诱发低频分量。因此，与主要表现出驻波特征的 CF 振动相比，IL 振动更为复杂。

当 Vr=16 时，肩部 CF 振动振幅增加到 0.22，相当于悬跨段内最大振幅的 48.54%。在 IL 方向，肩部峰值振幅上升到 0.0746，远超一阶模态下观察到的值。两个肩部的峰值振幅同时出现并同步过渡到波谷。类似于 Vr=11 的情况，这些波峰和波谷的空间位置随时间不断移动，表明具有明显的行波特征。此外，双峰表现出数值差异，主峰在左右肩部之间周期性交替，而波谷表现出类似的空间动态。

总之，在所有振动模态下，包括过渡模态，肩部均存在可测量的 IL 振动。在二阶模态下，IL 振动的行波特征增大了跨段肩部的振荡振幅。这些区域的 IL 振幅显著超过悬跨段，表明管道沉降可能大幅增加。相比之下，肩部 CF 振动仅在二阶模态下变得显著。

图 27、图 28、图 29、图 30 进一步证实了上述讨论中 IL 振动表现出更宽的频带带宽。这些图还展示了在 Vr=7、10、11 和 16 时，悬跨段沿多个节点的 CF 和 IL 频率谱、位移时间历程和轨道轨迹。在不同的振动模态下，所有节点的 CF 振动主要表现出正弦特征。此外，CF 振动的频率内容相对简单，并保持在第一或第二固有频率锁定状态。这表明 CF 振动主要由涡脱落引起的交变升力主导。

图 27. 折减速度为 7 时各节点的 CF 和 IL 频率谱、位移时间历程和轨道轨迹。

图 28. 折减速度为 10 时各节点的 CF 和 IL 频率谱、位移时间历程和轨道轨迹。

图 29. 折减速度为 11 时各节点的 CF 和 IL 频率谱、位移时间历程和轨道轨迹。

图 30. 折减速度为 16 时各节点的 CF 和 IL 频率谱、位移时间历程和轨道轨迹。

当 Vr=7 时，IL 振动也主要表现出正弦行为，具有两个较小的振幅分量，导致所有节点都有明确的“8”字形轨道轨迹。低频分量为 0.53 Hz，与 CF 振动频率相同。因此，在第一锁定区，IL 振动不仅受波动阻力的倍频激励控制，还受到 CF 振动的调制，从而产生一个 CF 振动频率分量。这可归因于两个因素。一方面，大幅度的 CF 振动导致相对来流速度和攻角周期性变化。另一方面，由重力和沉降引起的管道曲率导致 CF 方向的结构 - 几何不对称，这也导致了 CF 频率分量。

当 Vr=10 时，IL 振动中出现显著的低频分量 (0.227Hz)，其振幅向跨中逐渐增加。该分量不是 CF 频率耦合分量。相反，它代表振动进入过渡区后多模态竞争引起的低频漂移。一阶和二阶 IL 模态交替出现，如图 26 (g) 所示，并表现出显著的行波特征。这也导致跨中轨道轨迹出现显著的周期性拉伸和失真。

当 Vr=11 时，IL 振动已完成模态转换，再次受倍频波动阻力和 CF 振动频率分量控制。远低于 CF 振动频率的低频分量几乎消失。在跨中，仅观察到一个 CF 振动频率两倍的分量，而在其他节点出现明显的 CF 振动频率分量。远离跨中的节点表现出新月形轨迹，这在之前的折减速度下未观察到。

随着 Vr 进一步增加，尾流不稳定性持续加剧，导致 IL 振动中出现更多低频分量。在 Vr=16 时，除跨中外，所有其他节点的 IL 振动均表现出两个低频分量，并具有“8”字形轨道轨迹。值得注意的是，靠近跨中的轨道轨迹演变成类似水滴的形状。

由于管道横截面上弯矩方向的不断变化，正应力分析不能局限于四个基点 (顶/底/左/右)。本研究采用两个关键点来表征管道应力行为：最大 RMS 正应力点和最大正应力波动幅度点。

图 31 显示了所有横截面上最大 RMS 应力和最大应力波动幅度的空间分布。在一阶模态下，最大 RMS 应力和最大应力幅度的峰值均发生在跨中和肩部段。具体而言，跨中在 Vr=8 时达到峰值 RMS 应力 139.10MPa，在 Vr=7 时达到峰值应力幅度 68.24MPa。在肩部，峰值 RMS 应力 99.30MPa 发生在距离 TDP 30m 处且 Vr=7，并在同一折减速度下距离 TDP 3m 处观察到峰值应力幅度 74.85MPa。分析证实，在一阶模态下，跨中承受最高的应力水平 (139MPa) 且具有相当大的波动幅度，而肩部段尽管应力水平适中 (30-70MPa)，却表现出最大的应力波动幅度 (74.85MPa)。

在空间分布上，TDP 附近 50m (0.83L3) 范围内的肩部段保持高应力水平，而 TDP 附近 10m (0.17L3) 范围内的区段表现出显著的波动幅值。下载：下载高分辨率图片 (938KB) 下载：下载全尺寸图片 图 31. 管道沿线的最大均方根应力和最大应力波动幅值分布。在二阶模态中，均方根应力和应力幅值均显著超过一阶模态中的观测值。在跨中附近，峰值均方根应力达到 217.45MPa，峰值应力幅值达到 169.09MPa，两者均发生在 Vr=16 时。在肩部，峰值均方根应力达到 180.09MPa，伴随峰值应力幅值 109.93MPa。肩部区段中均方根应力峰值的位置比一阶模态中更靠近 TDP。此外，在较高的折减速度下，即使是在距离 TDP 70m 处（对应 1.17 倍跨长 L3），仍可观察到量级为 10 MPa 的应力波动幅值。图 32 通过展示管道垂直位移的包络线、水平和垂直平面内的曲率包络线以及三维曲率包络线，解释了应力特征的根本原因。与 IL 方向相比，较高的 CF 振动幅值会引起垂直平面内更大的曲率变化。因此，应力波动的幅值主要由垂直平面内的曲率控制，例如悬空跨和过渡区的前段。下载：下载高分辨率图片 (909KB) 下载：下载全尺寸图片 图 32. 管道的位移和曲率包络线。在一阶模态中，垂直平面内的曲率显著超过水平平面内的曲率，表明管道应力主要由 CF 振动和垂直 PSI 控制。在肩部，均方根应力峰值发生在中间段。这是因为由管道横向振荡引起的显著动态沉降仅限于前段，导致中间段的垂直曲率更大。相比之下，在二阶模态中，较高的 CF 幅值和 IL 频率在两个区段中都产生了深度的管道沉降，从而降低了垂直曲率。二阶模态中较高的流速导致整个悬空跨和过渡区产生更大的水平曲率，从而产生较高的均方根应力。因此，在过渡区中，一阶模态下的均方根应力主要由垂直曲率控制，而二阶模态下主要由水平曲率控制，导致峰值位置不同。如果仅考虑 CF 方向振动，二阶模态中由增强的水平曲率引起的应力峰值位置偏移无法被准确捕捉。在以往采用简化土壤模型并仅考虑 CF 振动的研究 [14,60] 中，跨中肩部段的均方根应力通常被认为远高于跨中。然而，在充分考虑到 PSI 的非线性后，本研究得出了相反的结论。这是因为频繁的管道 - 土壤接触导致肩部段产生显著沉降，从而降低了曲率和应力。然而，这也导致跨中肩部段的应力幅值与跨中相当，并且在管道的更宽区域延伸出可观的应力幅值。因此，跨中肩部段的 PSI 对管道的应力分布有显著影响。在工程分析中，应仔细考虑 PSI 以获得准确的应力分布和疲劳损伤预测。本研究考察了一阶模态 (Vr=7) 和二阶模态 (Vr=16) 中表现出最大截面应力波动幅值的悬空跨截面，分析其均方根应力和波动幅值分布，如图 33 所示。关键观察表明，最大均方根应力位置和最大应力幅值位置并不重合。前者发生在朝向入流偏移的管顶处（一阶模态为 135.4MPa，二阶模态为 205.5MPa），对应的应力幅值分别为 63.68MPa 和 98.22MPa。相比之下，最大应力幅值观察于管底处（一阶模态为 68.2MPa；二阶模态为 169.1MPa），伴随较低的均方根应力值 80.04MPa 和 87.18MPa。因此，在评估由 VIV 引起的管道疲劳损伤时，评估应不仅考虑高应力水平区域，还应考虑应力波动幅值大的区域，即使平均应力相对较低。应特别关注朝向入流偏移的管顶和管底的疲劳损伤。下载：下载高分辨率图片 (764KB) 下载：下载全尺寸图片 图 33. 均方根应力和应力波动幅值分布。TDP 处的垂直 P-y 曲线如图 34 所示。在一阶模态中，归一化管道贯入深度相对较小，最大为 4.79。在二阶模态期间，贯入深度增加至 7.27，与一阶模态相比上升了 52%。下载：下载高分辨率图片 (353KB) 下载：下载全尺寸图片 图 34. 肩部处的垂直 P-y 曲线。距离 TDP 10m 范围内的肩部区段的垂直 P-y 曲线如图 35 所示。在一阶模态中，最大土壤抗力达到 10.19kN/m。在二阶模态中，峰值土壤抗力增加至 16.92kN/m，伴随高达 12.51kN/m 的吸力。在管道上抬阶段，10m 范围内的所有管道区段均受到显著的土壤吸力。下载：下载高分辨率图片 (443KB) 下载：下载全尺寸图片 图 35. Vr=7 和 Vr=16 时的垂直 P-y 曲线。4.3. 轴向张力的影响不同水平的轴向张力会改变管道的固有频率，从而改变其 VIV 特性。分析了四种轴向张力情况——1300、2300、3300 和 4300kN——以量化其对振动幅值、频率响应和应力分布的影响。最大幅值和主导频率如图 36 所示。随着轴向张力增加，一阶和二阶模态中的 CF 振动幅值均增加，而 IL 振动表现出最小的变化。在 Vr=7 时，最大 CF 幅值上升高达 0.10，增幅为 22.08%。此外，承受较高轴向张力的管道表现出更高的频率和更早地过渡到更高阶模态。下载：下载高分辨率图片 (730KB) 下载：下载全尺寸图片 图 36. 不同轴向张力下的最大幅值和主导频率。图 37 展示了在完全发展的一阶和二阶模态响应下，管道沿线最大均方根应力和最大应力波动幅值的分布。结果表明，较低的轴向张力由于 VIV 幅值较高而导致更大的应力波动。具体而言，降低轴向张力会使一阶模态中的两个峰值分别增加高达 38.06% 和 50.88%，二阶模态中分别增加高达 38.06% 和 30.01%。尽管均方根应力通常随较高的轴向张力而增加，但 1300kN 轴向张力情况由于较大的 VIV 幅值而表现出最大的均方根应力。下载：下载高分辨率图片 (858KB) 下载：下载全尺寸图片 图 37. 不同轴向张力下的最大均方根应力和最大应力波动幅值分布。4.4. 间隙比的影响如前分析所示，VIV 触发肩部深处的管道沉降。一旦 PSI 和 VIV 响应稳定下来，垂直贯入深度和侧向位移均停止发展，管道围绕其平均偏移位置振动。然而，以往的研究忽略了海床边界效应。本节研究在管道沉降稳定后，间隙比如何调节 VIV 行为。如图 38(a) 所示，最大 CF 振动幅值表现出三个不同的状态：(1) 当 Vr<6 时，幅值随间隙比的变化可忽略不计。(2) 在 7-10 的关键折减速度范围内，对应于一阶和过渡模态，最大幅值随较小的间隙比而减小。在 Vr=7、8、9 和 10 时，最大减小量分别为 0.0763、0.1285、0.1381 和 0.0745，对应的减小幅度分别为 13.53%、24.37%、34.29% 和 29.69%。(3) 当 Vr>10 时，间隙比的影响再次变得可忽略不计。在 IL 方向，当 Vr>5 时，间隙比对幅值有显著影响（图 38(b)）。一阶和二阶模态的最大幅值减小量分别高达 0.0179 (减小 27.29%) 和 0.0186 (减小 24.72%)。下载：下载高分辨率图片 (409KB) 下载：下载全尺寸图片 图 38. 不同间隙比情况下的最大振动幅值。主导频率的比较如图 39 所示，揭示出仅当一阶和二阶模态完全激发时，频率差异才可辨别。总体而言，主导频率随间隙比的减小而略微增加。下载：下载高分辨率图片 (469KB) 下载：下载全尺寸图片 图 39. 不同间隙比情况下的主导频率。在完全发展的一阶和二阶模态响应下，管道沿线最大均方根应力和最大应力波动幅值的比较如图 40 所示。在一阶模态期间，间隙比对均方根应力影响最小，但使肩部处的应力幅值峰值减小 12.66%，跨中处减小 13.03%。在二阶模态中，将间隙比从 5.0 减小到 0.5 会使肩部处的均方根应力峰值减小 9.22%，而肩部处的幅值峰值减小 7.89%，悬空跨处减小 7.82%。下载：下载高分辨率图片 (764KB) 下载：下载全尺寸图片 图 40. 不同间隙比情况下的最大均方根应力和最大应力波动幅值分布。总之，减小间隙比显著抑制了一阶和二阶模态中的 IL 振动幅值，减小量高达 27.29%。一阶和过渡模态中的 CF 幅值也有所减小，减小量高达 34.29%。此外，减小间隙比会略微增加主导振动频率并减小应力波动幅值。在一阶模态中，肩部段和跨中的应力幅值可减小 10% 以上，从而减轻这些区段的疲劳损伤。然而，在二阶模态中，应力幅值的减小不太明显。4.5. 动态管道沉降的影响本节展示了使用相同参数的 VIV 计算，排除了由管道横向振荡引起的动态管道沉降，并假设足够大的间隙比以消除边界效应。最大振动幅值如图 41 所示。动态沉降对 CF 振动的影响可忽略不计，仅在二阶模态中引起 IL 幅值的微小减小。两种情景下的主导振动频率几乎相同（图 42）。下载：下载高分辨率图片 (448KB) 下载：下载全尺寸图片 图 41. 排除沉降 (SE) 条件下的最大振动幅值。下载：下载高分辨率图片 (210KB) 下载：下载全尺寸图片 图 42. 排除沉降 (SE) 条件下的主导频率。贯入包络线如图 43 所示。当 Vr>4 时，动态管道沉降显著加深了肩部处的贯入深度，而对距离 TDP 30m 以外的管道区段影响很小。在 Vr=7 时，TDP 处的最大贯入深度比排除沉降 (SE) 情况高出 45.17%。下载：下载高分辨率图片 (752KB) 下载：下载全尺寸图片 图 43. 过渡区 (TDZ) 内的贯入包络线。在过渡模态中，当 CF 振动保持在一阶模态 (Vr<11) 时，考虑动态沉降主要加深了 TDP 附近的贯入。然而，一旦 CF 振动过渡到二阶模态 (Vr=11)，动态沉降就会显著加深整个过渡区的贯入。在距离 TDP 70m 的位置，最大贯入深度从 0.41Dout 增加至 1.81Dout。当 Vr>11 时，与排除沉降的情景相比，TDP 处的贯入深度增加了 66.26% 至 68.12%。因此，忽略动态沉降会严重低估二阶 CF 振动下的贯入深度。因此，当两端的轴向张力保持不变时，即使肩部段的沉降深度增加了 45.71%–68.12%，自由悬空管道仍可能表现出幅值和频率非常相似的 VIV。换句话说，完全考虑土壤强度退化和塑性变形的土壤边界条件对自由悬空管道的固有频率影响有限。连续振动最终会将肩部段的土壤重塑为不再改变的稳定状态。这种土壤状态仅修改肩部段的平衡配置和悬空跨的平衡位置。悬空跨的 VIV 特性仍主要由轴向张力、跨长、质量比、阻尼和弯曲刚度等因素主导。基于之前的分析，当 Vr 从 10 增加到 11 时，肩部段的 IL 振动幅值几乎不变，而 CF 振动经历了完全的模态转换。

因此，尽管动态沉降模型受横向振荡影响，但可以得出结论，最终贯入深度由垂直管线 - 土壤相互作用（PSI）控制——包括横向（CF）振动幅度、垂直载荷和土壤强度退化。图 44 显示了最大均方根应力和最大应力波动幅度。仅在一阶模态下的均方根应力和二阶模态下的幅度中观察到明显差异。在一阶模态下，肩部处的峰值均方根应力向触地点（TDP）移动并减少了 7.16%，同时悬跨部分的应力也略有下降。这种现象产生的原因是动态沉降仅增加了 TDP 附近的贯入深度，从而减小了局部曲率并增加了远端管段段的曲率。在二阶模态下，距离 TDP 15m 以外的区段应力波动幅度被大幅高估，增加了 5.20—17.40 MPa (44.2%—81.0%)。下载：下载高分辨率图像 (463KB) 下载：下载全尺寸图像 图 44. 排除沉降（SE）条件下的最大均方根应力和最大应力波动幅度分布。尽管忽略动态沉降对较大间隙比下的预测 VIV 幅度影响有限，但它会大幅降低对悬跨肩部沉降深度的预测并影响应力评估。具体而言，它会低估一阶模态下悬跨肩部和跨中的均方根应力，同时高估二阶模态下过渡区（TDZ）的应力波动幅度。更重要的是，当间隙比降至 2 以下时，VIV 特性发生显著变化，表明在小间隙比情况下必须考虑动态沉降。

4.6. 悬跨长度的影响

针对间隙比为 5 和 0.5 的情况分析了三种悬跨长度（40, 50, 60 m），而其他所有参数保持与前述分析一致。

图 45 显示了最大振动幅度。随着悬跨长度减小，CF 和 IL 幅度均呈现下降趋势。在一阶模态下，60m 跨度的峰值 CF 和 IL 幅度分别为 0.56 和 0.065，而 40m 跨度对应的数值分别降至 0.41 和 0.050。在二阶模态下，60m 跨度的峰值幅度分别为 0.45 和 0.075，而 40m 跨度则为 0.25 和 0.032。这些结果表明，幅度减小在二阶模态下更为显著。下载：下载高分辨率图像 (449KB) 下载：下载全尺寸图像 图 45. 不同悬跨长度情况下的最大振动幅度。

随着间隙比的减小，最大幅度的减小情况如图 46 所示。对于一阶模态下的 CF 振动：60m 跨度在 Vr=9 时实现了最大的抑制量 0.13 (32.61%)，而 40m 跨度在 Vr=11 时的最大减小量为 0.06 (24.37%)。在二阶模态下，所有三种情况的幅度减小量均保持在 3%–9% 的稳定范围内。因此，在一阶模态下，对间隙比的敏感性取决于跨度长度，较长的悬跨表现出更显著的 CF 幅度抑制。这种敏感性在二阶模态下可以忽略不计。对于 IL 振动，减小百分比总体上保持一致，无论 Vr、悬跨长度或振动模态如何，范围均在 15% 到 30% 之间。下载：下载高分辨率图像 (387KB) 下载：下载全尺寸图像 图 46. 不同悬跨长度情况下的最大幅度减小量。

随着悬跨长度的减小，主频显著增加，而一阶模态的锁定区明显变宽（图 47）。对于 40m 的情况，频率跳跃延迟到 Vr=12 时才发生。下载：下载高分辨率图像 (365KB) 下载：下载全尺寸图像 图 47. 不同悬跨长度情况下的主频。

图 48 比较了完全激发振动模态时沿管线方向的最大均方根应力和最大应力波动幅度。在一阶模态下，缩短跨度使跨中处的峰值均方根应力最多增加 25.89%，而肩部处最多减少 20.75%。相比之下，在二阶模态下，较短的跨度使跨中和肩部处的峰值均方根应力分别最多增加 88.94% 和 33.54%。均方根应力的这种显著放大是因为较短的跨度具有更高的固有频率，导致在相同的折减速度下具有更高的速度和更大的水平曲率。然而，如此极端的海底流速在实际海洋条件下很少达到。下载：下载高分辨率图像 (855KB) 下载：下载全尺寸图像 图 48. 不同悬跨长度情况下的最大均方根应力和最大应力波动幅度分布。

将悬跨长度从 60m 减小到 40m 会使一阶模态下跨中处的峰值应力波动幅度增加 35.45%，而肩部处仅产生轻微增加。然而，在二阶模态下，趋势相反，肩部处的峰值幅度减少了 38.01%。因此，改变悬跨长度会在两种模态下产生相反趋势的应力幅度变化。

对于所有三种情况，减小间隙比会持续降低一阶模态下的应力波动幅度。尽管如此，间隙比对二阶模态下的应力波动幅度以及两种模态下的均方根应力影响可忽略不计。

5. 结论

本研究推导了流动规则公式，并提出了一种针对横向振荡管线的新型动态沉降模型。此外，开发了一种综合方法，用于考虑非线性 PSI 和海底边界效应的大跨度悬跨管线 IL 和 CF 方向的 VIV 模拟。基于 VFIFE 方法和非线性 PSI 模型建立了悬跨管线模型。随后研究了悬跨管线的 VIV 特性，以及流速、轴向张力、间隙比和动态沉降的影响。主要结论如下：

1. 与 CF 振动相比，IL 振动幅度明显较低，在锁定区内，最大幅度仅达到相应 CF 幅度的 10%–25%。在所有振动模态下，IL 振动频率始终约为 CF 振动频率的两倍。此外，IL 振动表现出宽的频率带宽和复杂的光谱成分。在过渡模态期间，IL 振动中发生空间频率竞争，表现出明显的行波特性和沿管线方向轨道轨迹的显著变化。
2. 二阶模态表现出比一阶模态显著更高的截面应力。评估 VIV 引起的管线疲劳时，不仅应考虑高应力水平区域，还应考虑应力波动幅度较大的区域，即使平均应力相对较低。非线性垂直 PSI 模型减小了肩部区段的均方根应力，但增加了应力波动幅度。
3. 减小间隙比显著抑制了一阶和二阶模态下的 IL 振动幅度，减小量高达 27.29%。一阶和过渡模态下的 CF 幅度也有所减小，多达 34.29%。此外，较小的间隙比会略微增加主振动频率，并将峰值应力波动幅度减少 7.82% –13.03%。
4. 当间隙比超过边界效应阈值时，忽略横向管线振荡引起的动态沉降对 IL 和 CF 方向的最大幅度影响较小。然而，动态沉降显著加深了悬跨肩部的管线沉降，增加量高达 68.12%，且在二阶模态下其影响范围延伸至管线的更大部分。因此，在小间隙比条件下，边界效应至关重要。因此，在小间隙比条件下，必须考虑动态沉降及其导致的边界效应增强。
5. 本研究提出的动态沉降模型可以准确预测由横向振荡引起的管线沉降。然而，当土壤或管线参数改变时，某些模型参数需要重新拟合。与垂直载荷相关的参数表现出相对较广的适用性。

未来的工作将集中在更小的间隙比和多跨管线。随着间隙比减小到本研究考虑的范围（0.5–5.0）之外，可能会发生管线 - 海床接触。计划进行数值和实验研究以检查较小间隙比下的 VIV 响应，所得数据将用于改进尾流振子模型。由于 VFIFE 方法能够有效处理碰撞，当前模型结合改进的尾流振子模型，可以为研究悬跨管线中的碰撞问题提供一个计算高效的框架。此外，将开发适用于多跨管线 TDZ 的 PSI 模型，以支持未来对多跨管线 VIV 的研究。

CRediT 作者贡献声明

余建兴：写作 - 审阅与编辑，项目管理，资金获取，调查，方法论。刘晓伟：写作 - 原始草稿，写作 - 审阅与编辑，概念化，数据管理，形式分析，调查，方法论，验证，软件。于洋：项目管理，资源，软件，监督。金自航：调查，验证，可视化。李真面：调查，验证，形式分析。杨正龙：调查，形式分析，资金获取，项目管理，监督。

余建兴 | 刘晓伟 | 于洋 | 金自航 | 李真面 | 杨正龙

天津大学水利工程智能建造与运维全国重点实验室，天津 300072，中国