海洋立管的作用是将油气资源从海底输送到海上平台，是海洋油气行业中的关键水下设备。立管内部通常为混合了多种气体与液体的多相流体，该种流动情况被称为段塞流。由于立管结构的几何非线性^[1-2]以及段塞流不均匀的质量分布，使得输送段塞流的立管动态响应预报十分复杂。

考虑到流固耦合方法计算复杂且耗时较长，一般在立管设计中需要建立合适的段塞流数学模型来开展海洋立管总体响应预报。段塞流具有速度和密度随时间变化的特征^[3]，当内部流动速度一致且每部分段塞长度相等时，即为稳态段塞流。早期研究往往将立管内部流动简化为均匀流考虑^[4]。Patel等^[5]提出了一种密度随位置和时间的变化而呈正弦变化的稳态段塞流模型。当前国内对于段塞流的研究主要集中在严重段塞流的成型机理以及影响^[6]，稳态段塞流模型相关研究较少。

Shabana^[7]提出了一种有限元建模方法，即绝对节点坐标法(absolute nodal coordinate formula, ANCF)。该方法中节点坐标在全局坐标系中定义，有效地避免了复杂的坐标变换。此外，ANCF采用斜率向量代替角坐标的方式来描述结构的局部旋转^[8]，可以准确地描述刚体位移与柔性变形之间的耦合关系。通过将ANCF与任意拉格朗日-欧拉(arbitrary Lagrangian-Eulerian, ALE)描述结合，可以精确高效地描述具有大变形特征的介质^[9]。

为预测内部段塞流引起的立管结构振动，本研究在ALE-ANCF框架下建立了海洋立管与内部段塞流的有限元分析模型。以钢悬链立管为例，分析了内部稳态段塞流对立管动态响应的影响。

1 ALE-ANCF有限元理论 1.1 控制体ALE-ANCF描述

ALE采用物质坐标来描述控制体的边界和沿介质的物质流动，介质上任一点的物质坐标为该点与起始点之间的介质长度。由于控制边界s₁和s₂是随时间变化的，控制区间的位置及长度随之变化。为便于后续控制区间的描述，定义局部坐标ξ为^[10]：

$ \xi=\frac{2 s-s_{1}(t)-s_{2}(t)}{s_{1}(t)-s_{2}(t)} $

(1)

在任意时刻t时控制体内部任意物质点的空间坐标为：

$ \boldsymbol{r}=\boldsymbol{r}(\xi, t) $

(2)

控制体的广义坐标在ALE-ANCF框架下可表示为：$ \boldsymbol{w}_{e}=\left[\boldsymbol{r}_{1}^{\mathrm{T}}~~\left(\boldsymbol{r}_{1}^{\prime}\right)^{\mathrm{T}} ~~\boldsymbol{r}_{2}^{\mathrm{T}}~~\left(\boldsymbol{r}_{2}^{\prime}\right)^{\mathrm{T}} ~~s_{1} ~~s_{2}\right]$。控制体内部任意点的空间位置向量和斜率向量表示为：

$ \mathit{\boldsymbol{r}} = \mathit{\boldsymbol{Nw}},{\mathit{\boldsymbol{r}}^\prime } = {\mathit{\boldsymbol{N}}^\prime }\mathit{\boldsymbol{w}} $

(3)

式中：N为形函数矩阵；$ \boldsymbol{w}=\left[\boldsymbol{r}_{1}^{\mathrm{T}}~~\left(\boldsymbol{r}_{1}^{\prime}\right)^{\mathrm{T}} ~~\boldsymbol{r}_{2}^{\mathrm{T}}~~\left(\boldsymbol{r}_{2}^{\prime}\right)^{\mathrm{T}}\right]$。

将式(4)中的空间位置向量直接对时间求导即可得到该点的速度$ {\mathit{\boldsymbol{\dot r}}}$与加速度$ {\mathit{\boldsymbol{\ddot r}}}$为：

$ \begin{array}{c} \dot{\boldsymbol{r}}=\boldsymbol{N} \dot{\boldsymbol{w}}+\left(\frac{\partial \boldsymbol{N}}{\partial s_{1}} \dot{s}_{1}+\frac{\partial \boldsymbol{N}}{\partial \dot{s}_{2}} \dot{s}_{2}\right) \boldsymbol{w}\\ \ddot{\boldsymbol{r}}=\boldsymbol{N} \ddot{\boldsymbol{w}}+\left(\frac{\partial^{2} \boldsymbol{N}}{\partial s_{1}{ }^{2}} \dot{s}_{1}^{2}+\frac{2 \partial^{2} \boldsymbol{N}}{\partial s_{1} \partial {s}_{2}} \dot{s}_{1} \dot{s}_{2}+\frac{\partial^{2} \boldsymbol{N}}{\partial {s}_{2}{ }^{2}} \dot{s}_{2}^{2} \right) \boldsymbol{w}+\\ \left(\frac{\partial \boldsymbol{N} }{\partial s_{1}}\ddot{s}_{1}+\frac{\partial \boldsymbol{N}}{\partial s_{2}} \ddot{s}_{2}\right) \boldsymbol{w}+2\left(\frac{\partial \boldsymbol{N}}{\partial s_{1}} \dot{s}_{1}+\frac{\partial \boldsymbol{N}}{\partial \dot{s}_{2}} \dot{s}_{2}\right) \dot{\boldsymbol{w}} \end{array} $

(4)

1.2 稳态段塞流单元

采用Patel等^[5]提出的稳态段塞流模型来模拟内部流体密度的时空变化，如图 1所示，其中t、n分别为管道切向与法向，D_i为立管内径，v_i为段塞速度，l_s为段塞长度。该数学模型中段塞流密度分为恒定分量与时变分量：

$ \rho_{f}=\rho_{0}\left(1+\eta \mathrm{e}^{\mathrm{i}(k s-\omega t+\psi)}\right) $

(5)

式中：ρ_f和ρ₀分别为为内部流体密度和恒定密度分量；k=2π/l_s为段塞流波数；η为波动幅值系数；ω为段塞圆频率；ψ为初始相位，用于指定初始时刻的段塞流质量分布。

段塞流单元可以在管道结构的切向上任意移动，形函数仅约束其在管道结构法向的位移。假定内部流动为稳态段塞流，段塞速度v_i = ω/k，此时对于任意段塞流单元，物质的流入速度与流出速度相等。段塞流单元边界的物质坐标为段塞速度的函数：

$ s_{1}=-v_{i} t, s_{2}=L_{e}-v_{i} t $

(6)

式中L_e为单元长度。将式(6)代入式(4)，即可得到段塞流单元内部任意点的速度与加速度。

段塞流单元的平衡方程为：

$ \boldsymbol{M}_{f} \ddot{\boldsymbol{w}}+\boldsymbol{Q}_{i}=0 $

(7)

式中：Q_i为广义惯性力；M_f为单元质量阵。根据达朗贝尔原理，Q_i与M_f可由虚功方程得到：

$ \begin{aligned} &\boldsymbol{Q}_{i}=\frac{L_{e}}{2} \int\limits_{-1}^{1} \rho_{f} A_{f} \boldsymbol{N}^{\mathrm{T}}\left[\frac{4 \boldsymbol{v}_{i}}{L_{e}} \cdot \frac{\partial \boldsymbol{N}}{\partial \xi} \dot{\boldsymbol{w}}+\right. \\ &\left.\left(\frac{4 \boldsymbol{v}_{i}^{2}}{L_{e}^{2}} \cdot \frac{\partial^{2} \boldsymbol{N}}{\partial \xi^{2}}+\frac{2 \dot{\boldsymbol{v}}_{i}}{L_{e}} \cdot \frac{\partial \boldsymbol{N}}{\partial \xi}\right) \boldsymbol{w}\right] \mathrm{d} \xi \end{aligned} $

(8)

$ \boldsymbol{M}_{f}=\frac{L_{e}}{2} \int\limits_{-1}^{1} \rho_{f} A_{f} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N} \mathrm{d} \xi $

(9)

1.3 立管单元

与段塞流单元不同，立管单元上的节点物质坐标保持恒定。将$ \dot{s}_{1}=\dot{s}_{2}=0$代入式(4)中可得到管道单元上任意节点的速度与加速度。采用欧拉-伯努利梁假设，管道单元的平衡方程表示为：

$ \boldsymbol{M}_{r} \ddot{\boldsymbol{w}}+\boldsymbol{Q}_{f}+\boldsymbol{Q}_{e}=0 $

(10)

式中：M_r、Q_f与Q_e分别为管道单元的质量矩阵、广义外力以及广义粘弹性力。根据达朗贝尔原理，M_r、Q_f与Q_e由虚功方程得：

$ \boldsymbol{M}_{r}=\frac{L_{e}}{2} \int\limits_{-1}^{1} {\rho}_{r} A_{r} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N} \mathrm{d} \xi $

(11)

$ \boldsymbol{Q}_{f}=-\frac{L_{e}}{2} \int\limits_{-1}^{1} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{f} \mathrm{d} \xi $

(12)

$ \boldsymbol{Q}_{e}=\frac{L_{e}}{2} \int\limits_{-1}^{1}\left[\left(\frac{\partial \varepsilon}{\partial \boldsymbol{w}}\right)^{\mathrm{T}} E_{r} A_{r} \varepsilon+\left(\frac{\partial \kappa}{\partial \boldsymbol{w}}\right)^{\mathrm{T}} E_{r} I_{r} \kappa\right] \mathrm{d} \xi $

(13)

式中：ρ_r为立管材料密度；A_r为立管截面积；f为单元外力矢量；ε和κ分别为单元的轴向应变和曲率，定义为：

$ \varepsilon=\left|\boldsymbol{r}^{\prime}\right|-1, \kappa=\frac{\left|\boldsymbol{r}^{\prime} \times \boldsymbol{r}^{\prime \prime}\right|}{\left|\boldsymbol{r}^{\prime}\right|^{3}} $

(14)

单元的刚度矩阵可由Q_e对w直接求导得到^[11-12]:

$ \boldsymbol{K}=\frac{L_{e}}{2} \frac{\partial}{\partial \boldsymbol{w}} \int\limits_{-1}^{1}\left[\left(\frac{\partial \varepsilon}{\partial \boldsymbol{w}}\right)^{\mathrm{T}} E_{r} A_{r} \varepsilon+\left(\frac{\partial \kappa}{\partial \boldsymbol{w}}\right)^{\mathrm{T}} E_{r} I_{r} \kappa\right] \mathrm{d} \xi $

(15)

式(15)右端两部分分别表示刚度矩阵的拉伸分量和弯曲分量。

施加在立管单元上的外力包括重力、浮力、拖曳力和惯性力，后两项水动力通过莫里森公式计算：

$ \left\{\begin{aligned} &\boldsymbol{f}_{D}= \frac{1}{2} \rho_{w} D\left[C_{D n} \boldsymbol{S}_{3 \times 3}(\boldsymbol{u}-\dot{\boldsymbol{r}})\left|\boldsymbol{S}_{3 \times 3}(\boldsymbol{u}-\dot{\boldsymbol{r}})\right|+\right.\\ &\ \ \ \ \ \ \ \ \ \ \left.C_{D t} \boldsymbol{T}_{3 \times 3}(\boldsymbol{u}-\dot{\boldsymbol{r}})\left|\boldsymbol{T}_{3 \times 3}(\boldsymbol{u}-\dot{\boldsymbol{r}})\right|\right] \\ &\boldsymbol{f}_{I}= \rho A_{o}\left[\boldsymbol{S}_{3 \times 3}\left(C_{I n} \dot{\boldsymbol{u}}+\left(C_{I n}-1\right) \ddot{\boldsymbol{r}}\right)+\right.\\ &\ \ \ \ \ \ \ \ \ \ \left.\boldsymbol{T}_{3 \times 3}\left(C_{It} \dot{\boldsymbol{u}}+\left(C_{t}-1\right) \ddot{\boldsymbol{r}}\right)\right] \end{aligned}\right. $

(16)

式中：D为立管外径；C_Dn、C_Dt分别为法向与切向的拖曳力系数；C_In、C_It分别为法向与切向的惯性力系数；u和${\mathit{\boldsymbol{\dot u}}} $分别为外部流体的速度与速度；S_3×3与T_3×3为与切矢量t₀相关的转换矩阵，表示为：

$ {\mathit{\boldsymbol{t}}_0} = \frac{{{\mathit{\boldsymbol{r}}^\prime }}}{{\left| {{\mathit{\boldsymbol{r}}^\prime }} \right|}},{\mathit{\boldsymbol{T}}_{3 \times 3}} = {\mathit{\boldsymbol{t}}_0}{\mathit{\boldsymbol{t}}_0}^{\rm{T}},{\mathit{\boldsymbol{S}}_{3 \times 3}} = {\mathit{\boldsymbol{I}}_{3 \times 3}} - {\mathit{\boldsymbol{T}}_{3 \times 3}} $

(17)

1.4 数值求解方法

静态平衡方程的求解采用Newton-Raphson方法，迭代形式为^[13]：

$ \begin{gathered} \boldsymbol{w}_{(n+1)}=\boldsymbol{w}_{(n)}+\Delta \boldsymbol{w}_{(n)} \\ \Delta \boldsymbol{w}_{(n)}=\left[\boldsymbol{K}\left(\boldsymbol{w}_{(n)}\right)\right]\left(\boldsymbol{Q}_{e}-\boldsymbol{Q}_{f}+\boldsymbol{Q}_{r}\right) \end{gathered} $

(18)

式中：n为迭代次数，当Δw满足：

$ \left\|\Delta \boldsymbol{w}_{(n)}\right\| \leqslant \boldsymbol{w}_{0} $

(19)

时迭代终止。

在t+Δt时刻，输送段塞流的立管的平衡方程为：

$ \boldsymbol{M}^{t+\Delta t} \ddot{\boldsymbol{w}}^{t+\Delta t}+\boldsymbol{K} \boldsymbol{w}^{t+\Delta t}=\boldsymbol{Q}_{f}^{t+\Delta}-\boldsymbol{Q}_{r}^{t+\Delta} $

(20)

式(20)可采用Newmark法进行求解，在该方法中，t+Δt时刻的位移、速度与加速度为：

$ \left\{\begin{array}{l} \boldsymbol{w}^{t+\Delta t}=\boldsymbol{w}^{t}+\Delta t \dot{\boldsymbol{w}}^{t}+\left[(0.5-\beta) \Delta t^{2}\right] \ddot{\boldsymbol{w}}^{t}+\beta \Delta t^{2} \ddot{\boldsymbol{w}}^{t+\Delta t} \\ \dot{\boldsymbol{w}}^{t+\Delta t}=\dot{\boldsymbol{w}}^{t}+[(1-\gamma) \Delta t] \ddot{\boldsymbol{w}}^{t}+\gamma \Delta t \ddot{\boldsymbol{w}}^{t+\Delta t} \\ \ddot{\boldsymbol{w}}^{t+\Delta t}=\frac{1}{\beta \Delta t^{2}}\left(\boldsymbol{w}^{t+\Delta t}-\boldsymbol{w}^{t}\right)-\frac{1}{\beta \Delta t} \dot{\boldsymbol{w}}^{t}-\left(\frac{1}{2 \beta}-1\right) \ddot{\boldsymbol{w}}^{t} \end{array}\right. $

(21)

将式(21)代入t+Δt时刻的平衡方程，可以得到递推公式：

$ \widetilde{\boldsymbol{K}}^{t+\Delta t} \boldsymbol{w}^{t+\Delta t}=\widetilde{\boldsymbol{Q}}^{t+\Delta t} $

(22)

其中：

$ \widetilde{\boldsymbol{K}}^{t+\Delta t} =\frac{1}{\beta \Delta t^{2}} \boldsymbol{M}^{t+\Delta t}+\boldsymbol{K} $

$ \begin{array}{c} \widetilde{\boldsymbol{Q}}^{t+\Delta t} =\boldsymbol{Q}_{f}^{t+\Delta}-\boldsymbol{Q}_{r}^{t+\Delta}+\\ \boldsymbol{M}^{t+\Delta t}\left[\frac{1}{\beta \Delta t^{2}} \boldsymbol{w}^{t}+\frac{1}{\beta \Delta t} \dot{\boldsymbol{w}}^{t}+\left(\frac{1}{2 \beta}-1\right) \ddot{\boldsymbol{w}}^{t}\right] \end{array} $

2 稳态段塞流对立管的影响分析 2.1 立管与段塞流参数

以钢悬链立管(SCR)为例，分析稳态段塞流对立管的影响。图 2给出了钢悬链立管模型的示意图，该立管总长1 260 m，工作水深为900 m，水平跨距700 m，顶部悬挂点位于水面下10 m。立管外径D=0.4 m与内径D_i=0.36 m，空气中单位长度的质量为189.1 kg/m，法向拖曳力系数与惯性力系数分别为1与2。

钢悬链立管内部输送流体为API 48原油和甲烷气体，内部流体平均密度选为632.135 kg/m³，液相密度与气相密度分别为790 kg/m³与0.675 kg/m³，波动幅值系数为0.25，段塞速度为2.5 m/s，段塞长度L_s=80 m。

为了探讨稳态段塞流对钢悬链立管的影响，本研究还对输送内部均匀流的立管进行了模拟，该内部均匀流模型具有与稳态段塞流相同的平均密度与流动速度，2种内部流动情况下的内部流体总质量相同。

2.2 静态分析

在分析稳态段塞流对立管的动态响应影响之前，先简要分析内部均匀流与段塞流情况下立管的静态特性。依据单个段塞的长度，将钢悬链立管划分为315个等长单元。采用土壤刚度来描述立管与海床之间的相互作用^[14]。图 3给出了2种内部流动情况下的立管张力与弯矩特性曲线。

通过观察图 3可以发现2种内部流动情况下的立管张立曲线几乎完全重合，这说明段塞流的不均匀密度分布对立管的张力影响很小。2种内部流动情况下的弯矩曲线具有相似的趋势，但是内部稳态段塞流情况下的弯矩曲线呈现出明显的局部波动。其原因是钢悬链立管的截面抗弯刚度极大，导致很小的曲率变形就会引起很大的弯矩变化。

2.3 动态分析

在动态仿真中，通过对立管的顶部节点施加X方向的强迫位移来模拟顶部平台的波频运动。激励位移r_x可以表示为：

$ r_{x}=\lambda \sin \left(\frac{2 {\rm{ \mathsf{ π} }}}{T_{m}} t\right) $

(23)

式中：顶部激励位移幅值λ和周期Τ_m分别取为0.5 m和9 s。动态仿真时长设置为2 000 s, 时间步长为0.05 s。

图 4给出了钢悬链立管顶部张力的时历曲线。可以看出在内部稳态段塞流的情况下，立管顶部张力的变化呈现出双频叠加的特性。与内部均匀流情况相比，稳态段塞流情况下顶部张力的波动幅值有所增大，2种情况下顶部张力幅值之间的差值为3.78 kN，仅为内部均匀流情况下顶部静态张力的0.27%，这说明稳态段塞流对钢悬链立管的动态张力影响很小。

图 5给出了2种内部流动情况下沿管长弯矩分布时历图。比较图 5(a)、(b)可以发现弯矩在稳态段塞流情况下的变化幅值明显大于内部均匀流情况。立管触地段的弯矩在2种流动情况下的变化趋势相近，而在立管悬挂段，弯矩的变化频率与段塞流频率相近。

为更明显地观察弯矩变化规律，分别选择触地点附近的节点35(距锚点140 m处)与悬挂段弯矩最大处附近的节点75(距锚点300 m处)作为观测点，图 6给出了上述2节点处弯矩变化(距锚点140 m处)的时历曲线。可以看出，稳态段塞流情况下上述两节点的弯矩变化均呈现出双频叠加的特性，节点35处的弯矩变化以顶部激励频率为主，而节点75处的弯矩变化以稳态段塞流的频率为主，这说明稳态段塞流对悬挂段弯矩的影响更为显著。

位移计算结果如图 7所示。结果显示，段塞流对节点35的X向位移几乎没有影响，而其Z向位移会随着段塞流的不均匀质量移动发生波动，2种内部流动情况下的最大Z向位移值相差0.01D。稳态段塞流情况下节点75的X与Z向位移均出现了双频叠加的特性，与均匀流动情况下的最大位移值相差分别为0.07D和0.08D。上述分析表明，稳态段塞流对立管触地段影响极小，但是会改变其Z向的运动特性。而对于立管悬挂段，稳态段塞流的影响显著，这也是该段弯矩变化特性发生改变的原因。

3 段塞长度参数分析

为分析段塞长度对立管动态响应的影响，对输送不同长度段塞流的钢悬链立管进行了时域仿真。段塞长度与管径之间的比值选取为L_s/D=100，160，220，280。内部流体平均密度与流动速度保持不变图 8和图 9分别给出了不同段塞长度下的节点35与节点75处的位移时历曲线。

由图 8(a)可知，段塞长度的改变对触地点X方向的位移没有明显影响。对于图 8(b)，可以发现随着段塞长度的逐渐增加，触地点的Z向位移幅值显著增大。此外，在液塞段经过时，触地点的Z向位移受到抑制，而气塞段经过时则与之相反，这一现象是由液塞与气塞交替流动导致的立管内部流体质量不断变化而引起的。由图 9(a)、(b)可以看出，段塞长度的改变会导致立管悬挂段位移幅值的变化，随着段塞长度的逐渐增加，悬挂段X与Z方向位移的双频叠加特性愈为明显。

不同段塞长度下的弯矩时历曲线如图 10所示。由图 10(a)可以看出，随着段塞长度的增加触地点弯矩明显增大，其变化的双频叠加特性更加明显，当L_s/D=100, 280时，弯矩平均值分别为101.06 kN与127.53 kN, 两者之间差值为触地点静态弯矩的36%。对于图 10(b)，随着段塞长度的增加，悬挂段弯矩的波动幅值逐渐增大，其变化的主导频率与段塞流频率相接近。

4 结论

1) 稳态段塞流的不均匀密度分布对立管的张力影响很小，但是对立管悬挂段的弯矩影响显著，弯矩变化主导频率与段塞流频率相近。

2) 稳态段塞流会导致钢悬链立管位移幅值的增大，液塞与气塞交替流动导致的内流质量变化会改变触地点垂向的运动特性。

3) 随着段塞长度的增加，立管悬挂段位移的双频叠加特性更为明显，触地点处的弯矩显著增大。