空间变化海床阻力下的管道爬行：基于随机场的分析框架

《Marine Structures》：Pipeline walking under spatially variable seabed resistance: An analytical framework with random fields

　　

•提出了一种用于管线行走的分析-随机框架。
•轴向阻力以对数正态随机场表示。
•虚拟锚点由累积摩擦阈值求解。
•Abaqus/UEL模型验证了该公式，误差在0.036%以内。
•摩擦变异性增大分散度和上百分位响应。

1. 引言

海底管线对于油气输送至关重要。在深水环境中，管线通常直接铺设于海床上，因运营及维护检修所需的启停，承受高压高温（HPHT）变化的循环作用[1]。这些启停循环引起管线的轴向膨胀与收缩，导致管线行走的稳定性问题，即轴向位移逐渐累积的现象[2,3]。管线行走可能导致海底基础设施的损坏，例如汇管损坏或相连钢悬链式立管张力丧失，从而危及海底系统的完整性[4,5]。

管线行走可由不同机制引起，包括管线端部持续张力、海床坡度、热瞬变效应和多相流效应[6]。为量化这些机制，已开展了广泛的分析、数值和实验研究[4,[7], [8], [9], [10], [11], [12], [13], [14]]。在此框架下，Carr等人[4]提供了由上述四种激活机制驱动的管线行走解析解。Tian等人[12]在刚塑性土体框架中纳入了持续端部张力和海床坡度，Tian等人[14]进一步将该模型扩展至弹塑性土体阻力。

大多数现有研究采用确定性建模框架，假设管土相互作用参数沿管线长度方向空间均匀分布。管线沿线管土相互作用空间均匀性的假设与实际的南海床条件不一致。海底管线通常延伸数公里，且经常穿越非均质海床条件。沉积过程、应力历史、安装扰动、埋深和海床地形形态的变化导致土体强度和界面阻力存在显著的空间变异性[[15], [16], [17], [18]]。

现场观测和侧向屈曲的概率研究表明，空间变异性可根本性改变结构响应。White等人[19]表明，管线倾向于在局部薄弱区域屈曲，而非由整体上限阻力控制，凸显了局部变异性而非整体行为的重要性。Westgate等人[20]也指出，铺设时的埋深表现出多尺度空间变异性，Sicilia等人[21]强调了使用基于可靠性的评估的必要性，因为确定性评估可能导致错误预测或有偏设计。

管线沿线的土体空间变异性可建模为空间相关随机场，以同时捕捉局部变异性和空间连续性。随机场理论在基于可靠性的岩土工程中被广泛采用，为量化岩土不确定性提供了严格框架[17,[22], [23], [24]]。

本研究开发了一种用于管线行走的分析-随机框架，显式地纳入了轴向管土阻力的纵向变异性。等效轴向摩擦系数以对数正态随机场表示，给定均值和变异系数（COV），其空间依赖性通过底层高斯场的波动尺度δ来指定。Tian等人[12]的控制平衡关系被重构为累积摩擦阻力形式，使得加热和冷却阶段的虚拟锚点（VAP）成为依赖于特定实现、由唯一阻力阈值交点定义的随机变量。因此，主要贡献不仅在于向确定性解中添加随机场，而是从摩擦场统计到VAP重定位、VAP间距变异和行走速率概率分布的直接力学映射。该公式通过随机场和蒙特卡罗收敛性检验以及独立的Abaqus/UEL对比进行了评估。随后用于量化COV、波动尺度和自相关函数形式对全响应分布（包括其上尾）的影响，并阐明基于百分位的预测如何改变减灾决策。该框架旨在作为计算高效、面向机理的概率评估工具，与经现场标定的有限元分析和项目特定设计相辅相成，而非取而代之。

2. 管线行走解法回顾

Tian等人[12]提出了一种管线轴向行走的解析解，将海床阻力建模为刚塑性摩擦定律。管线理想化为位于倾斜海床上的轴向可变形杆件，如图1所示，并假设在加热-冷却循环过程中始终保持横向稳定。

下载：下载高分辨率图像（118KB） 下载：下载全尺寸图像
图1. 管线计算模型及符号约定。

在刚塑性假设下，单位长度轴向土体阻力表示为：

(1) $f(x) = \mu w \cos\theta \cdot \text{sign}(v(x))$

其中 $\mu$ 为摩擦系数，$w$ 为单位长度管线重量，$v(x)$ 为轴向速度，$\theta$ 为海床坡度。函数 $\text{sign}(\cdot)$ 表示滑动方向，取值为+1或?1，使摩擦力始终与运动方向相反。

考虑微分管线单元的轴向平衡，轴力 $N(x)$ 满足：

(2) $\frac{dN}{dx} = f(x) - w\sin\theta$

由此可得：

(3) $N(x) = \int_0^x (f(x) - w\sin\theta) \, dx + F_t = \int_x^L (f(x) - w\sin\theta) \, dx + F_b$

其中 $F_t$ 和 $F_b$ 分别表示管线顶部和底部端部的张力，$L$ 为管线总长度。

轴向应变 $\varepsilon(x)$ 由热膨胀分量和弹性分量组成：

(4) $\varepsilon(x) = \alpha \delta T + \frac{N(x)}{EA}$

其中 $\alpha$ 为线膨胀系数，$\delta T$ 为温度变化，$E$ 为杨氏模量，$A$ 为管线截面积。

轴向位移 $u(x)$ 由应变积分得到：

(5) $u(x) = \int_{x_S}^x \varepsilon(x) \, dx + u_S$

其中 $u_S$ 为任意选取的点S（坐标为 $x_S$）的位移。

对于整体平衡，总土体阻力平衡驱动力：

(6) $\int_0^L f(x) \, dx = wL\sin\theta + F_b - F_t = W\sin\theta + F_b - F_t = Q$

其中 $W = wL$ 为管线总重量。$Q = W\sin\theta + F_b - F_t$ 被称为"激活力"。

在循环热载荷作用下，稳态行走速率 $\Delta u$ 为：

(7) $\frac{\Delta u}{L} = \alpha\left(\delta T - \mu W\cos\theta - \frac{|Q|}{\alpha EA}\right) \cdot \frac{Q}{\mu W\cos\theta}$

其中 $\delta T$ 为温度变化量。

Tian等人[14]进一步将式(7)扩展为更一般的形式：

(8) $\frac{\Delta u}{L} = \alpha\left(\delta T - \mu W\cos\theta - \frac{|Q|}{\alpha EA}\right) \cdot \frac{(x_C - x_H)}{L}$

其中 $x_C$ 和 $x_H$ 分别表示冷却和加热阶段形成的锚点轴向位置。该表达式也与Carr等人[4]的结果一致。

3. 通过引入土体阻力空间变异性对解析解的扩展

在实际海洋条件下，海床起伏和土体固有非均质性导致轴向土体阻力沿管线呈非均匀分布。为表征这些条件的局部变异性和空间连续性，随机场蒙特卡罗方法已被用于岩土可靠性分析中[17,24,25]。在本公式中，轴向阻力的空间变异性通过空间变化的等效摩擦系数 $\mu(x)$ 引入，同时假设淹没状态下单位长度管线重量 $w$ 为常数。

3.1. 对数正态随机场的构建

轴向管土阻力受埋深、土体强度、海床地形、安装扰动和循环重塑效应等因素影响[15]。在工程实践中，由此产生的轴向阻力通常通过等效摩擦系数 $\mu$ 表示，定义为轴向阻力与有效淹没管线重量之比，如DNV-RP-F114[26]所采用。对于本研究中考虑的外露或浅埋管线，轴向阻力假设主要在管道埋深和接触区域附近的表层土中动员[15,26]。因此，等效摩擦场 $\mu(x)$ 代表该有效上部管土相互作用层沿纵向的变异性，而非完整的垂直地层。该简化假设认为相互作用层厚度足够，能够包含管土接触区和埋深区。如果管道埋深穿越层间界面，或沿管线路径遇到不同的近地表土层单元，则需要采用逐层特定的、分段平稳的或非平稳的随机场。

从概率角度而言，$\mu$ 代表多个乘性随机机制的综合结果。此类乘性过程通常用对数正态变量表示[27]，对数正态随机场也已被用于概率管线分析中描述空间变化的土体强度[28]。采用对数正态表示主要是为了确保 $\mu(x) > 0$；并非旨在暗示对数正态性普遍适用于轴向摩擦，也未通过特定场地的测量数据进行验证。

按照既有岩土随机场实践[22,23]，采用二阶平稳性作为海床条件名义均质路段的局部理想化。在此假设下，底层高斯场 $Y(x)$ 的均值和方差为常数，其自相关仅取决于空间间距。大尺度趋势和不同土层单元之间的过渡需在使用本公式之前予以去除或单独建模。

因此，将 $\mu(x)$ 建模为空间相关的对数正态随机场，提供了物理上合理且计算高效的管土阻力变异性表征，同时保留了管线设计实践中采用的等效摩擦表示方法[26]。

对数正态场 $\mu(x)$ 通过底层高斯随机场 $Y(x)$ 构建，使得：

(9) $\mu(x) = \exp(Y(x))$

等价地，$Y(x) = \ln \mu(x)$；此后 $Y(x)$ 表示底层高斯场。底层高斯场 $Y(x)$ 的参数与 $\bar{\mu}$ 和 $\sigma^2$ 的关系为：

(10) $\mu_Y = \ln\left(\frac{\bar{\mu}}{\sqrt{1+COV^2}}\right)$

(11) $\sigma_Y^2 = \ln(1+COV^2)$

其中 $\bar{\mu}$ 为摩擦系数均值，$\sigma$ 表示标准差，$COV = \sigma / \bar{\mu}$ 为摩擦系数的变异系数。这些变换确保所得对数正态场 $\mu(x)$ 具有预设的均值和方差。

底层高斯场 $Y(x)$ 的空间依赖关系用单指数自相关函数表示，这是岩土随机场分析中常用的简约模型[16,22]：

(12) $\rho_Y(\Delta x) = \exp\left(-\frac{2|\Delta x|}{\delta}\right)$

其中 $\Delta x = x_i - x_j$ 为两个位置之间的有符号间距，$\delta$ 为 $Y(x)$ 的波动尺度，此处定义为其自相关函数的双侧积分：$\delta = \int_{-\infty}^{\infty} \rho_Y(r) \, dr$。对于式(12)中的指数模型，当 $|r| = \delta/2$ 时，$\rho_Y(r) = 1/e$。第4.4节使用两种替代模型考察了预测的管线响应对所假设自相关结构的敏感性。由于 $\mu(x)$ 通过式(9)的非线性变换获得，其自相关函数与 $Y(x)$ 的自相关函数之间的关系为 $\rho_\mu(r) = \frac{\exp[\sigma_Y^2 \rho_Y(r)] - 1}{\exp(\sigma_Y^2) - 1}$。因此，变换后的场 $\mu(x)$ 通常不严格保持与 $Y(x)$ 相同的自相关函数或波动尺度。全文中 $\delta$ 特指为 $Y(x)$ 预设的波动尺度。

高斯场在 $x_i$ 和 $x_j$ 处的协方差因此为：

(13) $C_{ij} = \text{Cov}[Y(x_i), Y(x_j)] = \sigma_Y^2 \exp\left(-\frac{2|x_i - x_j|}{\delta}\right)$

在均匀空间网格上，指数协方差对应于一阶平稳高斯马尔可夫过程。因此，可以使用以下自回归递推精确生成该高斯场：

(14) $\begin{cases} Y_1 = \mu_Y + \sigma_Y \xi_1, \\ Y_i = \mu_Y + \varphi(Y_{i-1} - \mu_Y) + \sigma_Y \sqrt{1-\varphi^2} \, \xi_i, & i=2,\ldots,n_x \\ \varphi = \exp\left(-\frac{2h}{\delta}\right) \end{cases}$

其中 $h = L/(n_x - 1)$ 为均匀网格间距，$\xi_i \sim \mathcal{N}(0,1)$ 为独立标准正态随机变量。采用式(14)中的平稳初始化后，生成的场满足 $\text{Cov}[Y_i, Y_j] = \sigma_Y^2 \exp(-2|x_i - x_j|/\delta)$，从而精确再现预设的指数协方差。该递推表示既不需要全局协方差矩阵分解，也不需要模态截断，且每次实现的计算成本随空间点数线性增长。因此，特别适合本研究考虑的大空间网格和蒙特卡罗样本量。

管线用 $n_x = 5000$ 个均匀分布的点在 $L = 5000$ m上表示，$h = 1.0002$ m。对于每次实现，式(14)使用新采样的独立标准正态随机变量序列进行计算，所得高斯场通过式(9)逐点变换为预设的对数正态摩擦场。

累积阻力函数 $S(x)$ 通过对离散化摩擦场进行累积梯形积分获得，VAP位置通过线性插值在包含预设阻力阈值的相邻网格点之间求得。采用 $h = 5, 2, 1, 0.5$ m的网格间距进行了空间离散收敛性检验，各分辨率使用相同的随机数集以最小化蒙特卡罗抽样造成的差异。相对于 $h = 0.5$ m的参考解，$h = 1$ m离散化对 $\Delta u$ 的均值、标准差和95百分位数的影响分别仅为0.0005%、0.0006%和0.0227%。因此，后续所有计算均采用约1 m的网格间距。

3.2. 段平均摩擦系数

在完全动员的库仑摩擦下，单位长度轴向土体阻力可表示为 $f(x) = \mu(x) \, w\cos\theta \cdot \text{sign}(v(x))$。根据Carr等人[4]和Tian等人[12]，在加热或冷却阶段，管线响应由虚拟锚点（VAP）表征，该点将管线分为两个滑动方向相反的滑动区域。因此，在每个区域内 $v(x)$ 的符号恒定，该区域的合力摩擦仅取决于 $\mu(x)$ 在相应区间上的空间积分。

因此，与其沿管线积分摩擦系数 $\mu(x)$，不如采用段平均摩擦系数更为方便，其值等于各滑动区域的摩擦合力（见图2a）。

设管道长度为 $L$，坐标 $x \in [0, L]$。

**图2. 代表性实现的VAP判定：** (a) 摩擦场及分段平均值；(b) 归一化累积摩擦与阈值交点，用于定义VAP。

整条管道的平均摩擦系数定义为：

(15) $\overline{\mu}_{OL} = \frac{1}{L}\int_{0}^{L} \mu(x)\,dx$

在升温阶段，加热VAP $x_H$ 将管道分为 OH 和 HL 两段，对应的平均值定义为：

(16) $\overline{\mu}_{OH} = \frac{1}{x_H}\int_{0}^{x_H} \mu(x)\,dx$, $\overline{\mu}_{HL} = \frac{1}{L - x_H}\int_{x_H}^{L} \mu(x)\,dx$

类似地，在降温阶段，冷却VAP $x_C$ 将管道分为 OC 和 CL 两段，对应的平均值为：

(17) $\overline{\mu}_{OC} = \frac{1}{x_C}\int_{0}^{x_C} \mu(x)\,dx$, $\overline{\mu}_{CL} = \frac{1}{L - x_C}\int_{x_C}^{L} \mu(x)\,dx$

这些定义确保每段的摩擦合力被等价表示，例如：

$$\int_{0}^{x_H} \mu(x)\, w\cos\theta\, dx = \overline{\mu}_{OH}\, w\cos\theta\, x_H$$
$$\int_{x_H}^{L} \mu(x)\, w\cos\theta\, dx = \overline{\mu}_{HL}\, w\cos\theta\,(L - x_H)$$

冷却阶段有类似的关系式。由构造方式可知，全长平均值与加热VAP所诱导的分段一致：

(18) $\overline{\mu}_{OL}\, L = \overline{\mu}_{OH}\, x_H + \overline{\mu}_{HL}\,(L - x_H)$

式(16)-(18)表明，每一段的平均摩擦系数精确地保留了相应空间变化段的积分摩擦合力。因此，在每个实现中，该表示在轴向平衡和VAP判定层面具有严格的力等价性。这种等价性涉及段的积分阻力，并不意味着局部土壤阻力或位移场的逐点等价。

### 3.3 升温阶段与加热VAP的确定

在热膨胀过程中，如图3(a)所示，管道在加热VAP $x_H$ 两侧沿相反方向滑动，因此土壤阻力在该位置改变符号。

**图3. 由虚拟锚点(VAP)分隔的土壤阻力分布：** (a) 升温过程；(b) 降温过程。

假设完全动员的刚塑性摩擦，管道的轴向平衡为：

(19) $\overline{\mu}_{OH}\, w\cos\theta\, x_H + Q = \overline{\mu}_{HL}\, w\cos\theta\,(L - x_H)$

式(19)表达了由加热VAP分隔的两个滑动段之间的轴向力平衡。

联立式(18)和(19)得：

(20) $\frac{x_H}{L}\,\overline{\mu}_{OH} = \frac{1}{2}\left(\overline{\mu}_{OL} - \frac{Q}{W\cos\theta}\right)$

代入 $\overline{\mu}_{OH}$ 的定义式(16)，得到积分条件：

(21) $\int_{0}^{x_H} \mu(x)\,dx = \frac{L}{2}\left(\overline{\mu}_{OL} - \frac{Q}{W\cos\theta}\right)$

式(21)决定了摩擦场的任一实现中加热VAP $x_H$ 的位置。

### 3.4 降温阶段与冷却VAP的确定

在热收缩过程中，滑动方向与升温阶段相反，如图3(b)所示。此时出现冷却VAP $x_C$，将两个滑动区域分隔开。

按照相同的推导过程，轴向平衡为：

(22) $\frac{x_C}{L}\,\overline{\mu}_{OC} = \frac{1}{2}\left(\overline{\mu}_{OL} + \frac{Q}{W\cos\theta}\right)$

代入 $\overline{\mu}_{OC}$ 的定义得：

(23) $\int_{0}^{x_C} \mu(x)\,dx = \frac{L}{2}\left(\overline{\mu}_{OL} + \frac{Q}{W\cos\theta}\right)$

式(23)定义了冷却VAP $x_C$。

### 3.5 VAP的存在性与数值确定

3.3节"升温阶段与加热VAP的确定"和3.4节"降温阶段与冷却VAP的确定"中的推导假设在升温或降温过程中存在一个虚拟锚点(VAP)，将管道分隔为两个滑动区域。因此，有必要考察此类点存在的数学条件以及解是否唯一。

对于摩擦场 $\mu(x)$ 的给定实现，定义累积摩擦函数：

(24) $S(x) = \int_{0}^{x} \mu(x)\,dx$

由于摩擦系数沿管道满足 $\mu(x) > 0$，累积函数 $S(x)$ 在区间 $x \in [0, L]$ 上连续且严格递增。

因此，

(25) $S(0) = 0$, $S(L) = \int_{0}^{L} \mu(x)\,dx = \overline{\mu}_{OL}\, L$

加热和冷却VAP的控制方程，即式(21)和(23)，可紧凑地写为：

(26) $S(x_H) = \frac{L}{2}\left(\overline{\mu}_{OL} - \frac{Q}{W\cos\theta}\right) = R_H$

(27) $S(x_C) = \frac{L}{2}\left(\overline{\mu}_{OL} + \frac{Q}{W\cos\theta}\right) = R_C$

其中，$R_H$ 和 $R_C$ 分别表示升温和冷却情况下的累积轴向阻力阈值；因此，$x_H$ 和 $x_C$ 分别是满足 $S(x_H) = R_H$ 和 $S(x_C) = R_C$ 的唯一位置。

由于 $S(x)$ 严格递增，当且仅当对应方程的右端落在区间 $[0, \overline{\mu}_{OL}\, L]$ 内时，$x_H$ 和 $x_C$ 的解存在且唯一。该要求直接导出以下条件：

(28) $\frac{|Q|}{W\cos\theta} \leq \overline{\mu}_{OL}$

当满足此条件时，图2(b)中所示的归一化累积摩擦曲线 $S(x)/L$ 与归一化阈值 $R_H/L$ 和 $R_C/L$ 在唯一坐标 $x_H/L$ 和 $x_C/L$ 处相交。这些交点确定加热VAP $x_H$ 和冷却VAP $x_C$。

从物理角度看，条件(28)确保激活力的大小不超过管道的总可用摩擦承载力。当此条件成立时，管道响应由由内部VAP分隔的两个滑动区域所表征，如图3(a)和(b)所示。若 $\frac{|Q|}{W\cos\theta} > \overline{\mu}_{OL}$，则式(21)或式(23)在 $x \in [0, L]$ 内无法满足，不存在内部VAP。此时管道将沿由 $Q$ 的符号所规定方向发生整体滑动，而非形成两个相对的滑动区域。在蒙特卡罗实现中，该条件在评估VAP位置之前进行检查，任何违反式(28)的实现均被归类为整体滑动情况。

对数正态构造确保 $\mu(x)$ 始终严格为正；因此，随机场仅改变可用轴向阻力的大小和空间分布。摩擦力的符号由局部滑动方向通过 $\text{sign}[v(x)]$ 单独确定。海床坡度和端部张力通过 $Q = W\sin\theta + F_b - F_t$ 和法向力因子 $\cos\theta$ 进入公式。若不同的确定性荷载组合导致 $Q < 0$，则 $x_H$ 和 $x_C$ 的相对排序及行走方向一致地反转。这代表了一个物理上可接受的激活力方向反转，而非非物理的摩擦模式。

为清晰起见，逐实现的VAP判定按五个步骤进行。首先，生成分散值 $\mu(x)$ 并评估全长均值 $\overline{\mu}_{OL} = S(L)/L$。其次，通过累积梯形积分计算 $S(x)$。第三，评估升温和冷却阻力阈值：

$$R_H = \frac{L\left(\overline{\mu}_{OL} - \frac{Q}{W\cos\theta}\right)}{2}, \quad R_C = \frac{L\left(\overline{\mu}_{OL} + \frac{Q}{W\cos\theta}\right)}{2}$$

第四，检验式(28)以验证两个阈值均在 $[0, S(L)]$ 范围内。最后，$x_H$ 和 $x_C$ 分别作为满足 $S(x_H) = R_H$ 和 $S(x_C) = R_C$ 的唯一坐标获得。在离散实现中，每个VAP坐标通过在包围对应阈值的两个相邻采样点之间进行线性插值来确定。$S(x)$ 的严格单调性保证了每个允许的交点的唯一性。

确定性极限直接源于此构造。当摩擦场空间均匀时，即 $\mu(x) = \mu$，此时 $\overline{\mu}_{OL} = \mu$，$S(x) = \mu x$。式(26)和(27)则简化为：

$$\frac{x_H}{L} = \frac{1 - \frac{Q}{\mu W\cos\theta}}{2}, \quad \frac{x_C}{L} = \frac{1 + \frac{Q}{\mu W\cos\theta}}{2}$$

因此，$\frac{x_C - x_H}{L} = \frac{Q}{\mu W\cos\theta}$，这是Tian等人[12]获得的确定性VAP间距。

### 3.6 虚拟锚点的确定与后续响应评估

一旦确定了给定空间变化摩擦场 $\mu(x)$ 实现下的升温和冷却虚拟锚点(VAP)，即 $x_H$ 和 $x_C$，即可利用第2节开发的解析框架确定管道的完整力学响应。这些响应包括升温和冷却阶段中的轴向土壤阻力 $f(x)$、轴向力 $N(x)$、应变 $\varepsilon(x)$ 和位移 $u(x)$。

由于 $S(x)$ 在 $[0, L]$ 上连续且严格递增，它存在唯一的逆函数 $S^{-1}$。因此，虚拟锚点位置可显式写为：

(29) $x_H = S^{-1}(R_H)$, $x_C = S^{-1}(R_C)$

条件(28)保证两个阈值均在 $S^{-1}$ 的定义域内。对于离散化随机场，$S^{-1}(R_H)$ 和 $S^{-1}(R_C)$ 通过在包围对应阈值的相邻网格点之间进行线性插值来评估。

循环热载荷下的管道行走速率 $\Delta u$ 然后通过连续升温和冷却阶段之间的净位移累积来计算。结合式(8)与逐实现的平均摩擦系数得：

(30) $\frac{\Delta u}{L} = \alpha\left(\delta T - \frac{\overline{\mu}_{OL}\, W\cos\theta + |Q|}{\alpha E A}\right)\frac{x_C - x_H}{L}$

式(30)并非针对任意空间变化摩擦场的精确闭合解。段的积分摩擦合力和由此产生的VAP位置被精确保留，而后续的位移评估则通过逐实现的全长平均摩擦系数和VAP间距来表示空间场。为检验当摩擦变异性较小时此近似是否引入偏差，令：

$$\mu(x) = \mu_0 + \varepsilon \tilde{\mu}(x)$$

其中 $E[\tilde{\mu}(x)] = 0$，$\varepsilon$ 与COV成正比。对式(26)-(30)进行一阶展开得：

$$E[\Delta u] = \Delta u^{(0)} + O(\varepsilon^2)$$

因为 $\overline{\mu}_{OL}$、$x_H$ 和 $x_C$ 的一阶摄动均为零均值场 $\tilde{\mu}(x)$ 的线性泛函，因此其集合均值为零。因此，当变异性趋于零时，力等价程序对集合平均行走速率不引入一阶偏差；任何均值偏移从 $O(\text{COV}^2)$ 开始。对于较大的变异性，高阶效应变得越来越重要，当行走响应主要由累积阻力和VAP间距而非单个滑动段内的高度局部化波动所支配时，该近似预计最为准确。

在概率分析中，上述过程对大量 $\mu(x)$ 实现进行重复，以获得行走速率和VAP位置的统计分布。在零变异性极限下，$\overline{\mu}_{OL} = \mu$ 且 $(x_C - x_H)/L = Q/(\mu W\cos\theta)$。将这些关系代入式(30)可恢复式(7)，证实当空间变异性消失时，本公式精确地退化为确定性解析解。

### 3.7 评估管道行走的计算框架

本研究采用的整体计算框架如图4所示。对于空间变化摩擦场 $\mu(x)$ 的每个蒙特卡罗实现，该程序将确定性管道参数、热载荷和海床坡度与以对数正态随机场建模的随机空间变异性（以 $\mu(x)$ 的均值和COV以及 $Y(x)$ 的波动尺度 $\delta$ 为特征）进行整合。首先评估该实现的有效平均摩擦系数 $\overline{\mu}_{OL}$，然后从累积阻力条件中数值确定虚拟锚点 $x_H$ 和 $x_C$。在确定VAP位置后，利用第2节推导的闭合形式解析关系计算轴向土壤阻力 $f(x)$、轴向力 $N(x)$、应变 $\varepsilon(x)$ 和位移 $u(x)$。随后，从升温和冷却阶段之间的位移差获得净行走速率 $\Delta u$。在大量实现的集合上重复此过程，得到VAP位置和行走速率 $\Delta u$ 的概率分布，从而在统一的蒙特卡罗框架内将随机空间变异性与确定性解析力学联系起来。

**图4. 基于蒙特卡罗的随机管道行走分析框架流程图。**

## 4. 概率分析与讨论

本节对空间变化轴向摩擦引起的管道行走进行随机评估的案例研究。采用对数正态随机场来建模管道长度方向上轴向摩擦系数 $\mu(x)$ 的空间变异性。通过蒙特卡罗模拟评估行走速率 $\Delta u$ 的概率分布及相应的虚拟锚点位置。

随机分析中使用的基准计算参数如表1所示。在每个蒙特卡罗集合中，所有管道、材料、热载荷和荷载参数，包括海床坡度和端部张力，均保持固定，仅采样空间摩擦场 $\mu(x)$。$\mu(x)$ 的COV和 $Y(x)$ 的波动尺度 $\delta$ 在参数组之间变化，而非在单个集合内部变化。因此，每个案例报告的响应变异性完全归因于轴向摩擦的规定空间变异性。对于所采用的值 $\theta = 2°$、$F_t = 100\,\text{kN}$ 和 $F_b = 300\,\text{kN}$，激活力为 $Q = 409.4\,\text{kN} > 0$，在研究的参数范围内未发生基准行走方向的反转。

**表1. 随机管道行走分析的计算参数。**

| 参数 | 数值 |
|------|------|
| 管道截面积，$A$ | 0.02 m2 |
| 管道长度，$L$ | 5000 m |
| 管道浮重，$w$ | 1200 N/m |
| 温度变化，$\delta T$ | 100°C |
| 杨氏模量，$E$ | 2.07 × 1011 Pa |
| 泊松比，$\nu$ | 0.3 |
| 热膨胀系数，$\alpha$ | 1.165 × 10?? °C?1 |
| 海床坡度，$\theta$ | 2° |
| 上端张力，$F_t$ | 100 kN |
| 下端张力，$F_b$ | 300 kN |
| 平均摩擦系数，$\bar{\mu}$ | 0.5 |
| $\mu$ 的变异系数，COV | 0.20；0.05–0.50（参数范围） |
| $Y$ 的波动尺度，$\delta$ | 150 m；30–500 m（参数范围） |
| 每个COV-$\delta$组合的蒙特卡罗实现数，$N_{MC}$ | 10,000 |

确定性管道和荷载参数对应于典型高温海底管道，而所研究的COV和 $\delta$ 范围参考了既有的岩土和管道随机场研究[20,23,28]，用于参数敏感性分析，而非针对特定 offshore 场地的校准统计量。对于本文考虑的基准情况和参数范围，没有任何实现违反条件(28)；因此，所有报告的蒙特卡罗统计量均对应VAP控制区域，排除了整体滑动情况。

### 4.1 随机场生成与响应公式的验证

首先检验使用式(14)中平稳自回归递推生成的轴向摩擦随机场的统计一致性。如图5(a)所示，代表性位置处的经验概率密度函数与规定的对数正态分布高度吻合，表明边缘统计量被准确再现。

图5(b)中的多个代表性实现示例说明了μ(x)沿管线长度的空间波动情况，而图5(c)中的累积积分S(x)则反映了空间变化的阻力逐步累积的过程，该过程主导着虚拟锚点的形成。

下载：下载高分辨率图像（577KB）
下载：下载完整尺寸图像
图5. 生成的对数正态轴向摩擦随机场的统计验证：(a)边缘概率密度函数对比；(b)μ(x)的样本实现；(c)累积积分S(x)；(d)lnμ(x)的自相关：经验值与指数模型。

空间相关结构进一步通过Y(x)的经验自相关函数进行评估，该函数由200个独立随机场实现估计得到，如图5(d)所示。经验平均自相关函数紧密跟随Y(x)的预设指数相关模型，该模型以δ为特征，而理论相关曲线在主相关范围内保持在经验估计值±1个标准差的包络线内。这些结果证实，生成的场能够准确重现μ(x)的预设边缘统计特性和Y(x)的目标自相关结构，从而为后续的随机行走分析提供了统计可靠的基础。

针对COV=0.2且δ=150 m的基线指数相关案例，考察了蒙特卡洛估计值的收敛性。如图6所示，Δu的运行均值、标准差、P5和P95在大约4000至5000次实现后趋于稳定。以全部10000次实现的估计值作为参考，在4000次实现时，均值、标准差、P5和P95的相对偏差分别为0.052%、0.306%、0.149%和0.129%。随着进一步抽样，仅出现微小的残余波动。因此，为每个研究的COV-δ组合独立采用10000次实现，以提供稳定的中心响应和分布尾部估计值。

下载：下载高分辨率图像（522KB）
下载：下载完整尺寸图像
图6. 蒙特卡洛响应统计量的收敛性：(a)Δu的运行均值；(b)运行标准差；(c)运行P5；(d)运行P95。

所提出公式在确定性极限下的一致性已在第3.6节中从理论上确立。其在有限摩擦变异性条件下的精度进一步通过图7(a、b)与使用Abaqus/Standard获得的独立时域有限元解的对比加以检验。有限元模型遵循Tian等[12]中采用的刚塑性轴向土阻力理想化假设。使用了与解析计算相同的5000 m管线和确定性加载参数，管线沿轴向以1 m的间距离散化。为每个区段分配了一个轴向阻力UEL，其极限阻力指定为Tmax,i = μi w cosθ ΔxFE，其中μi为对应随机场区段的摩擦系数，ΔxFE = 1 m。针对基线随机场案例解析行走速率分布的P10、P50和P90附近选取了三个代表性实现。

下载：下载高分辨率图像（196KB）
下载：下载完整尺寸图像
图7. 解析-随机公式与Abaqus/UEL的验证对比：(a)归一化的加热和冷却VAP位置；(b)所选随机场实现的行走速率。

分析了三个完整的加热-冷却循环，采用初始动员循环之后的响应进行对比。FE行走速率直接从全周期位移增量获得，而xH和xC则分别从加热和冷却位移增量场的零点独立识别。对于P10、P50和P90实现，解析行走速率分别为0.3155、0.3682和0.4239 m/cycle，与UEL结果0.3156、0.3682和0.4239 m/cycle相比，相应的相对偏差分别为0.036%、0.003%和0.009%，归一化VAP位置的最大偏差为L的0.0093%。高度一致性证实，所提出的随机场VAP公式能够重现独立增量有限元解所得到的实现依赖型VAP位置以及时域行走累积。

4.2. 管线行走的随机行为

图8、图9展示了所提出的基于随机场的框架在加热-冷却循环分析中的应用。图8展示了空间摩擦场的一个代表性实现的全周期响应。与确定性解析基准[12]相比，随机响应表现出由摩擦系数局部波动引起的明显空间振荡。这些振荡在土阻力分布中最为明显，随后通过空间积分传递至轴力、应变和位移场。对比验证了所提出的公式在传播空间摩擦不确定性贯穿整个热循环的同时，保留了平衡一致的场响应。

下载：下载高分辨率图像（481KB）
下载：下载完整尺寸图像
图8. 空间摩擦场代表性实现的加热-冷却循环响应：(a)土阻力；(b)轴力；(c)轴向应变；(d)轴向位移。

下载：下载高分辨率图像（458KB）
下载：下载完整尺寸图像
图9. 加热-冷却响应的随机包络线与确定性预测对比：(a)土阻力；(b)轴力；(c)轴向应变；(d)轴向位移。

图9进一步展示了从多个随机场实现中获得的随机包络线。灰色带量化了加热-冷却响应的空间不确定性。轴力和轴向应变的变异性在虚拟锚点xH和xC之间的区间内最为显著，而该区域之外的响应离散相对较小。与确定性解相比，随机包络线表明空间摩擦变异性主要放大加热和冷却VAP之间有效动员区内的响应离散。这些结果表明，随机框架不仅捕捉了局部波动，还捕捉了其对整体力学响应的累积影响。

基线随机响应中行走速率Δu的结果在图10中呈现，该结果针对COV=0.2且δ=150 m的指数相关案例，基于10000次蒙特卡洛实现。图10(a)中的经验概率密度函数以接近确定性预测值为中心：平均和中间行走速率分别为0.369和0.367 m/cycle，与确定性值0.368 m/cycle相比。然而，中心90%区间从P5=0.302 m/cycle延伸至P95=0.440 m/cycle。因此，P95值比确定性预测值大约高19.3%。图10(b、c)中的经验累积分布和箱线图进一步展示了不可忽视的响应分散。

下载：下载高分辨率图像（151KB）
下载：下载完整尺寸图像
图10. 基线行走速率的概率分布：(a)经验PDF；(b)经验CDF；(c)箱线图。

这些结果表明，尽管确定性解提供了合理的中心趋势估计，但单一确定性值无法量化由空间变异性摩擦引起的较大行走响应的可能性。因此，尾部百分位数以及（在项目特定许可行走限值可用的情况下）相应的超越概率应考虑在风险导向的管线评估中。

4.3. 摩擦变异性与波动尺度的影响

图11量化了摩擦系数的COV对行走响应的影响。如图11(a)所示，平均行走速率Δu随COV的增加而稳步上升，但该增幅的幅度相对温和。在研究的COV=0.05至0.50范围内，平均行走速率Δu从0.3684增加到0.3706 m/cycle，增幅仅为0.60%。

下载：下载高分辨率图像（191KB）
下载：下载完整尺寸图像
图11. 摩擦系数COV对行走响应的影响：(a)Δu的均值；(b)Δu的标准差；(c)Δu的第95百分位数。

在分散度度量方面观察到明显更强的敏感性。图11(b)显示，在同一COV范围内，Δu的标准差从0.0106增加到0.1020 m/cycle，增幅为866.32%。同时，图11(c)中的上尾统计量P95(Δu)从0.3856增加到0.5483 m/cycle，增幅为42.19%。均值弱增与分散度和上尾的强烈放大之间的这种差异表明存在非线性不确定性放大机制：低摩擦区段促进大的VAP偏移和大位移累积的效果远优于高摩擦区段抑制它们的效果。

从物理角度看，增加COV会强化管线沿线低摩擦区段与高摩擦区段之间的对比，从而提高出现促进热循环中更大累积位移的延伸弱区的概率。这些结果表明，摩擦变异性对响应中心趋势的影响有限，但对响应分散和上尾行为具有主导性影响，突出了COV作为随机管线行走评估中的关键不确定性参数。

波动尺度δ被指定给底层高斯场Y(x)，并通过对数正态变换控制μ(x)的沿线持续性，而非代表单个海床层的垂直厚度。对于等式(12)中采用的Y(x)的指数自相关函数，当|r|=δ/2时，ρY(r)=1/e；因此，所研究的δ=30–500 m范围代表Y(x)中15–250 m的相关距离。下限描述了浅层土强度和管线嵌入深度的相对局部化变化，而上限则允许与沉积单元转换、海床形态和安装相关嵌入趋势相关的空间一致变化[15,20]。这种多尺度沿线变异性已在HPHT管线的铺设后嵌入数据中被观察到[20]。因此，所选择范围用于本5 km管线的宽参数敏感性包络，不应解释为固有土强度的通用或场地校准范围。在项目特定应用中，δ应根据去除趋势的空间数据或反分析进行校准，以与底层高斯场Y(x)一致[22,23]；当存在不同土单元或突变的砂/泥转换时，应为管线路线的相应区段分别分配独立的平稳随机场，或采用非平稳场。

图12展示了δ对行走速率Δu统计特征的影响。如图12(a)所示，平均行走速率Δu随δ的增加单调递增，尽管增幅较小。在所研究的δ=30 m至500 m范围内，平均行走速率Δu从0.3688增加到0.3710 m/cycle，变化仅为0.60%。

下载：下载高分辨率图像（277KB）
下载：下载完整尺寸图像
图12. 波动尺度δ对行走响应的影响：(a)Δu的均值；(b)Δu的标准差；(c)Δu的第95百分位数。

相比之下，图12(b)展示了响应分散度更强的大幅增长。Δu的标准差从0.0194增加到0.0702 m/cycle，增幅为262.18%，而图12(c)显示P95(Δu)从0.4004增加到0.4906 m/cycle，增幅为22.55%。这些结果表明，波动尺度主要影响行走响应的分散和上尾，而非其均值。

较大的δ值增加了Y(x)的空间持续性，从而促进了延伸低摩擦区段的形成，减少了空间平均效应，并允许局部弱点更有效地影响全局响应。总体而言，δ对平均行走速率Δu的影响相对于其对分散度和上尾响应的影响仍然温和。这确定了波动尺度作为风险导向评估中的关键参数，特别是在极端行走而非平均行走主导工程决策时。表2总结了两个参数分析的典型数值统计量。对比进一步证明，COV和δ的变化对平均Δu的影响很小，而显著放大了其标准差和上尾统计量P95。

表2. 代表性随机场参数集的Δu统计摘要。

COV | δ (m) | Δu的均值(m/cycle) | Δu的标准差(m/cycle) | Δu的P95(m/cycle)
0.05 | 150 | 0.3684 | 0.0106 | 0.3856
0.20 | 150 | 0.3686 | 0.0418 | 0.4397
0.50 | 150 | 0.3706 | 0.1020 | 0.5483
0.20 | 30 | 0.3688 | 0.0194 | 0.4004
0.20 | 100 | 0.3692 | 0.0347 | 0.4264
0.20 | 500 | 0.3710 | 0.0702 | 0.4906

注：每个参数组使用10000次蒙特卡洛实现进行评估。

4.4. 对自相关结构的敏感性

尽管单指数自相关函数在岩土工程随机场建模中被广泛使用[16,22]，但相关结构通常无法从有限的场地数据中唯一确定。不同的自相关函数还意味着不同的样本路径平滑度，高斯和Whittle–Matérn模型已被用于表征土的空间变异性[29]。因此，通过对比指数模型、ν=1的Matérn模型和平方指数（高斯）模型，考察了预测的管线行走响应对此建模选择的敏感性。

为隔离自相关函数形状的影响，三个模型均被赋予相同的μ(x)对数正态边缘统计特性，而其底层高斯场均被赋予相同的波动尺度δ=150 m，遵循基于COV和波动尺度的随机场表征[22,23]。底层高斯场的自相关函数定义为：

?{ ρE(Δx) = exp(?2|Δx|/δ)
?{ ρM(Δx) = z·K?(z)，z = π|Δx|/δ
?{ ρG(Δx) = exp[?π(|Δx|/δ)2]

其中K?(·)为第二类修正贝塞尔函数。

这些定义满足2^∞∫_0 ρ(Δx)dΔx = δ，确保三个底层高斯场尽管具有不同的局部光滑度和衰减特性，但具有相同的波动尺度。对于本次敏感性对比，指数场是使用第3.1节所述的精确平稳AR(1)递推关系生成的，而Matérn和平方指数高斯场则是使用其各自协方差序列的快速傅里叶变换循环嵌入方法生成的。

如图13(a)所示，三个模型产生了明显不同的自相关结构。尽管如此，它们预测的行走速率累积分布却高度一致（图13(b)）。与指数模型相比，Matérn模型的均值、标准差和95百分位的Δu差异分别为0.045%、1.694%和-0.105%；高斯模型的对应差异分别为-0.025%、1.278%和-0.138%。加热和冷却虚拟锚点的均值差异小于0.04%；最大差异为2.09%，出现在平方指数模型的虚拟锚点间距标准差中。双样本柯尔莫哥洛夫-斯米尔诺夫检验发现，指数模型的累积分布与Matérn模型（KS统计量=0.0132，p=0.346）或平方指数模型（KS统计量=0.0120，p=0.465）之间不存在统计学上的显著差异。因此，在本文所考察的匹配的对数正态边际统计量和高斯场波动尺度条件下，自相关函数形状对累积阻力、虚拟锚点重定位及由此产生的管道行走响应的影响有限。

下载：下载高分辨率图像（219KB） 下载：下载全尺寸图像

图13. 管道行走响应对自相关模型的敏感性：(a) δ=150 m时，指数、Matérn (ν=1)和平方指数（高斯）自相关函数；(b) 10,000次实现得到的Δu累积分布。

4.5 控制虚拟锚点变异性机理

图14阐明了虚拟锚点（VAP）的随机变异性及其对行走响应的直接影响。归一化加热和冷却虚拟锚点位置x_H/L和x_C/L表现出明确的概率分布（图14(a,b)），表明空间摩擦随机性会引起管道沿线有效约束位置的系统性偏移。对于基准工况，x_H/L的均值为0.4317，标准差为0.0170，P5-P95范围为0.4038-0.4598。类似地，x_C/L的均值为0.5681，标准差为0.0170，P5-P95范围为0.5396-0.5956。这些数值表明虚拟锚点位置并非固定不变，而是随空间变化摩擦场的累积积分在管道长度的可观范围内发生变化。

下载：下载高分辨率图像（347KB） 下载：下载全尺寸图像

图14. 虚拟锚点位置变异性及其对行走响应的影响：(a) x_H/L的分布；(b) x_C/L的分布；(c) 随机Δu与确定性预测；(d) 累积摩擦力S(x)及对应的虚拟锚点位置。

这种位置变异性的后果反映在图14(c)所示的行走响应中。行走速率Δu的均值为0.3686 m/周期，标准差为0.0418 m/周期，P5-P95范围为0.3020-0.4397 m/周期，而确定性预测值为0.3685 m/周期。同时，归一化虚拟锚点间距(x_C - x_H)/L的均值为0.1365，标准差为0.0120，P5-P95范围为0.1169-0.1563，确定性值为0.1365。因此，确定性解接近随机响应的均值，但无法捕捉由虚拟锚点重定位引起的完整响应分散范围。

其内在机理在图14(d)中加以说明，图中累积摩擦力积分S(x)与定义x_H和x_C的激活水平相交。μ(x)的空间波动改变了S(x)的斜率，从而移动了决定虚拟锚点的交点。由于行走速率Δu由这些锚点位置之间的间距控制，摩擦的随机性通过虚拟锚点重定位传播，产生了观测到的随机响应。

综上所述，这些结果表明，管道行走的不确定性并非直接来源于局部摩擦变异性本身，而是来源于其对控制全局轴向变形的虚拟锚点位置和间距的集成效应。

4.6 缓解设计的示意性启示

为说明如何将预测的概率分布转化为缓解决策，将行走速率按管道长度归一化为Δu* = Δu / L。考虑示意性允许值Δu*_allow = 8.0×10??/周期，对应于本工况下0.4 m/周期的允许行走速率。

对于选定的设计百分位p，引入一个额外的有效轴向阻力F_A,p，其方向与激活力相反，使得Q_eff = Q - F_A,p。所需的最小阻力由P_p[Δu(Q_eff)/L] ≤ Δu*_allow确定。对于每个试验值F_A,p，使用10,000次蒙特卡罗实现重复基准随机场分析。所得缓解需求汇总于表3。确定性和P70响应在无需额外阻力的情况下满足示意性准则。相比之下，P80和P95响应超出允许值，分别需要额外的有效轴向阻力2.5和33.3 kN。因此，尽管确定性评估不会触发缓解措施，但采用更高的概率设计水平可以改变决策并量化相应的阻力需求增加。

表3. 基于允许归一化行走速率的示意性缓解需求。

| 设计依据 | Δu (m/周期) | Δu/L (×10??/周期) | F_A,p (kN) | Q_eff (kN) |
|---|---|---|---|---|
| 确定性 | 0.3685 | 7.37 | 0.0 | 409.4 |
| P70 | 0.3896 | 7.79 | 0.0 | 409.4 |
| P80 | 0.4029 | 8.06 | 2.5 | 406.9 |
| P95 | 0.4397 | 8.79 | 33.3 | 376.1 |

如果每个相同的锚固单元提供折减后的有效阻力R_A，则所需数量可估计为n_p = F_A,p / R_A。如果需要完全抑制行走，锚固系统必须提供F_A = Q = 409.4 kN，此时Q_eff = 0。

4.7 适用范围与局限性

本研究旨在作为基于虚拟锚点的解析解的初次随机场扩展，其范围由保持解析可行性的假设所定义。等效轴向摩擦力被表示为不随周期变化的一维平稳对数正态场，同时管道-土响应被理想化为刚塑性且纯轴向的，在预设的纵向海底坡面上。因此，该公式无法分辨非平稳地层过渡或三维海底地形、循环重塑导致的嵌入和阻力演化、弹塑性、软化或速率相关土体行为，以及轴向行走与侧向屈曲之间的相互作用。先前的研究[13,14]表明，其中一些效应会实质性地影响预测响应。此外，尽管UEL对比提供了对解析实现的独立数值验证，但假设的随机场统计量尚未针对现场特定数据进行标定，且目前尚无物理模型或现场观测数据可供参考。因此，结果应被解读为机理层面的概率性洞察，而非项目特定的设计预测。未来工作将把现场标定的非平稳或多变量随机场与能够在连续周期内更新管道-土体状态并解析真实三维海底上轴向-侧向耦合行为的增量有限元/UEL公式相结合，随后与管道-土体试验和现场监测进行验证。

5 结论

本研究建立了一个解析-随机框架，将空间变化的轴向管道-土阻力纳入基于虚拟锚点的管道行走解。累积阻力公式保留了控制加热和冷却虚拟锚点的分段积分摩擦力合力，同时利用实现相关的平均摩擦力和虚拟锚点间距来评估行走速率。当变异性趋于零时，该公式可恢复确定性解。对于有限变异性，与独立的Abaqus/UEL模型对比，行走速率最大差异为0.036%，虚拟锚点位置最大偏差为管道长度的0.0093%。

结果表明，空间变异性主要通过实现相关的虚拟锚点重定位影响管道行走。对于基准工况，平均行走速率接近确定性预测，但P95高出了19.3%。增加变异系数显著放大了响应分散性和上百分位响应，而波动尺度主要控制薄弱区域的空间相关性和由此产生的分散范围。在保持μ(x)相同边际统计量和底层高斯场相同波动尺度的条件下，Matérn和平方指数自相关模型使报告的响应统计量变化不超过2.09%，表明在所考察的平稳情况下自相关形状具有次要影响。在示意性允许行走准则下，P80和P95设计水平分别需要额外的有效轴向阻力2.5和33.3 kN，而确定性预测满足该准则。因此，该框架应被视为一种面向机理的概率工具，用于筛选和可靠性解读；项目特定的应用需要现场标定的随机场统计量、必要时采用更真实的本构模型和轴向-侧向建模，并与物理模型或现场观测进行验证。

贡献声明

吴万超：原稿撰写、可视化、软件、方法论、调查、形式分析、概念化。沈智超：审阅与编辑、监督、项目管理、方法论、资金获取、概念化。田英辉：审阅与编辑、验证、监督、方法论、概念化。

吴万超 | 沈智超 | 田英辉