刚性立管涡激振动的数值模拟方法研究

Research on Numerical Simulation Methods for Vortex-Induced Vibration of Rigid Risers

摘要: 在复杂的海洋环境中，深海立管易发生涡激振动，进而导致立管发生疲劳损伤。针对这一问题，本文采用高精度的四阶龙格库塔法和有限差分法，对刚性圆柱体结构振子与尾流振子耦合的常微分方程组进行数值模拟。数值算例通过Matlab实现数值模拟过程，结果显示在较长时间区间内计算得出的无量纲位移的振幅与无量纲尾流振子的振幅均与现有文献中的结果保持了较高的一致性，证实了所采用的数值方法能够有效地模拟刚性圆柱体的涡激振动响应，为未来在涡激振动耦合的偏微分方程组数值模拟及相关参数研究方面提供了可靠的工具。

Abstract: In the complex marine environment, deep-sea risers are susceptible to vortex-induced vibrations (VIV), which can lead to fatigue damage. To address this issue, this paper employs the high-precision fourth-order Runge-Kutta method and the finite difference method to numerically simulate the ordinary differential equations that couple the rigid cylindrical structure oscillator with the wake oscillator. The numerical examples are executed through Matlab to simulate the numerical process. The results indicate that over an extended period, the computed amplitudes of the dimensionless displacement and the dimensionless wake oscillator maintain a high degree of consistency with the results found in existing literature. This confirms that the numerical methods used are effective in simulating the vortex-induced vibration response of rigid cylinders, providing a reliable tool for future numerical simulations of coupled partial differential equations related to vortex-induced vibration and the study of associated parameters.

1. 引言

深海立管是海洋石油开采中的重要设备，在海洋工程领域，立管系统经常受到复杂洋流环境的影响。在复杂海洋环境条件下，围绕立管的流场会不断地形成漩涡，并且这些漩涡会从立管表面脱离。这种连续的漩涡生成与脱落过程导致立管受到周期性的流体力作用。当流体力的频率与立管系统的固有震动频率接近时，立管将进入一种被称为“涡激共振”的动态状态，这种现象被称作涡激振动 [1] 。涡激振动会加剧立管的疲劳损伤，长期影响可能会导致海洋立管的结构破坏，进而导致严重的相关事故。因此，对涡激振动进行数值预报具有重要意义。

目前，涡激振动的研究主要有三种方法：实验方法 [2] [3] [4] [5] ，半经验模型方法 [6] - [11] 和CFD模拟方法 [12] [13] [14] 。实验方法是通过实际的实验来分析涡激振动特性，但隔水管本身的长径比很大，难以对隔水管进行原型实验。CFD模拟方法主要是通过计算机求解纳维尔–斯托克斯方程，但其精度极其依赖计算机的性能。半经验模型方法融合了实验数据与理论分析，通过在风洞实验结果上构建模型，并合理挑选模型参数，以实现模型特性与实际情况的高度一致性。目前使用最广泛的是Van der Pol尾流振子模型，最先是由Bishop和Hassan [15] 提出，Hartlen和Currie [16] 建立了最具代表性的模型。在应用半经验模型方法时，最常用的思路是先将尾部流场产生的流体力描述为一个非线性方程，随后将流体的振子方程与圆柱体的结构振动方程进行耦合，采用较高精度的数值方法对耦合方程进行求解，从而对立管的涡激振动响应进行预测。目前学者对耦合方程常用的数值求解方法包括：有限差分法 [17] [18] 、有限元法 [19] 、积分变换法 [20] 等。

为了提高计算的精度、稳定性以及效率，本文主要采用高精度四阶龙格–库塔法与有限差分法对刚性圆柱体结构振子与尾流振子的耦合常微分方程组进行数值模拟，并将两种方法作对比。通过数值实验结果与已发表论文中的无量纲位移和无量纲振子进行对比分析来表明该数值方法可以很好地模拟刚性圆柱体的涡激振动响应特性，在实际工程中具有重要的指导意义。

2. 数学模型

刚性圆柱体在深海中涡激振动模型如图1所示，考虑一个长度为 $L$ 直径为 $D$ 的刚性圆柱体，圆柱体在均匀来流 $U$ 作用下引起横流方向发生涡激振动响应。对此涡激振动模型，将与海洋流速方向相同的方向取作 $x$ 轴，称为顺流振动方向；与海洋流速方向垂直方向取作 $y$ 轴，称为横流振动方向；将立管重力方向取作 $z$ 轴，称为铅直方向，同时 $x、y、z$ 三个方向形成右手直角坐标系，具体方向如图1所示。

Figure 1. Model of vortex-induced vibration for a rigid cylinder

图1. 刚性圆柱体涡激振动模型

假设用 $Y$ 表示圆柱体在 $y$ 方向产生横向位移，于是刚性圆柱体的结构振动方程为：

$m\frac{d^{2}Y\left(T\right)}{d{T}^2}+r\frac{dY\left(T\right)}{dT}+hY\left(T\right)=\frac{1}{4}\rho U^{2}D{C}_{L0}q(\; T\; )$

上式中， $T$ 为时间， $h$ 为系统刚度， $r$ 为系统阻尼包括结构阻尼 $r_s$ 和流体阻尼 $r_f$ ，可表示为 $r=r_s+r_f$ ， $r_f=\gamma {\omega }_f\rho D^2$ ， $\frac{1}{4}\rho U^{2}D{C}_{L0}q\left(T\right)$ 为横流方向的升力， $\rho$ 为流体密度， $q\left(T\right)$ 为尾流振子的运动， $C_{L0}$ 为横向升力系数， $m$ 为整个系统质量包括单位长度结构质量 $m_s$ 及附加流体质量 $m_f$ ，可表示为 $m=m_s+m_f$ ， $m_f=C_M\rho D^2\pi /4$ ， $C_M$ 表示附加质量系数 [21] ，对于圆柱体取 $C_M=\mathrm{1.0.}$

采用改进的Van der pol方程 [22] 来表示尾流振子的非线性特性，可表示为：

$\frac{d^{2}q\left(T\right)}{d{T}^2}+\epsilon {\omega }_f\left(q^2\left(T\right)-1\right)\frac{dq\left(T\right)}{dT}+{\omega }_f^{2}q\left(T\right)=\frac{A}{D}\frac{d^{2}Y\left(T\right)}{d{T}^2}$

上式中， $\epsilon$ 为非线性项中的小参数， $A$ 为结构对流体的耦合动力参数。

于是得到刚性圆柱体涡激振动流固耦合方程如下：

$\{\begin{array}{c}m\frac{d^{2}Y\left(T\right)}{d{T}^2}+r\frac{dY\left(T\right)}{dT}+hY\left(T\right)=\frac{1}{4}\rho U^{2}D{C}_{L0}q\left(T\right)\\ \frac{d^{2}q\left(T\right)}{d{T}^2}+\epsilon {\omega }_f\left(q^2\left(T\right)-1\right)\frac{dq\left(T\right)}{dT}+{\omega }_f^{2}q\left(T\right)=\frac{A}{D}\frac{d^{2}Y\left(T\right)}{d{T}^2}\end{array}$ (1)

立管的初始条件为

$\{\begin{array}{c}{Y|}_{T=T_0}=0,{\frac{dY}{dT}|}_{T=T_0}=0\\ {q|}_{T=T_0}=O\left({10}^{-3}\right),{\frac{dq}{dT}|}_{T=T_0}=0\end{array}$ (2)

3. 求解方法

首先对原方程组进行无量纲化转化

$y=Y/D$ ， $t=T\cdot {\omega }_f$

于是耦合方程(1)转化为

$\{\begin{array}{c}\frac{d^{2}y(t)}{d{t}^2}+\left(2\delta \xi +\frac{\gamma }{\mu }\right)\frac{dy(t)}{dt}+{\delta }^{2}y(t)=Mq(t)\\ \frac{d^{2}q(t)}{d{t}^2}+\epsilon (q^2(t)-1)\frac{dq(t)}{dt}+q(t)=A\frac{d^{2}y(t)}{d{t}^2}\end{array}$ (3)

上式中， $\delta =\frac{1}{St{U}_r},\mu =\frac{m}{\rho D^2},M=C_{L0}/\left(16{\pi }^{2}S{t}^2\mu \right).$

初始条件(2)式转化为

$\{\begin{array}{c}{y|}_{t=t_0}=0,{\frac{dy}{dt}|}_{t=t_0}=0\\ {q|}_{t=t_0}=O({10}^{-3}),{\frac{dq}{dt}|}_{t=t_0}=0\end{array}$ (4)

假设 $y\left(t\right)、q\left(t\right)$ 是耦合方程(3)的解，计算时间区间为 $\left[0,T\right]$ ，将时间区间进行 $N$ 等分，于是时间步长为 $\Delta t=T/N$ ，时间节点可表示为 $t_n=n\Delta t\left(n=0,1,\cdots ,N\right).$

3.1. 四阶龙格–库塔法

四阶龙格库塔法是一种高精度的单步算法，基本思想是利用当前点的信息和斜率的加权平均来预测下一个点的值。采用四阶龙格库塔方法对耦合方程(3)式进行求解，令 $\tilde{y}\left(t\right)=\frac{dy\left(t\right)}{dt}$ ， $\tilde{q}\left(t\right)=\frac{dq\left(t\right)}{dt}$ ，于是(3)可化为如下一阶微分方程组

$\{\begin{array}{c}\frac{dy\left(t\right)}{dt}=\tilde{y}\left(t\right)\\ \frac{d\tilde{y}\left(t\right)}{dt}=-a\tilde{y}\left(t\right)-by\left(t\right)+cq\left(t\right)\\ \frac{dq\left(t\right)}{dt}=\tilde{q}\left(t\right)\\ \frac{d\tilde{q}\left(t\right)}{dt}=-\epsilon \left(q^2\left(t\right)-1\right)\tilde{q}\left(t\right)-q\left(t\right)+A\left(-a\tilde{y}\left(t\right)-by\left(t\right)+cq\left(t\right)\right)\end{array}$ (5)

其中 $a=2\delta \xi +\frac{\gamma }{\mu }$ ， $b={\delta }^2$ ， $c=M.$

初始条件(4)式化为

$\{\begin{array}{c}y\left(t_0\right)=0\\ \tilde{y}\left(t_0\right)=0\\ q\left(t_0\right)=O\left({10}^{-3}\right)\\ \tilde{q}\left(t_0\right)=0\end{array}$ (6)

令 $P\left(t\right)={\left[y\left(t\right),\tilde{y}\left(t\right),q\left(t\right),\tilde{q}\left(t\right)\right]}^T$ ， $P^{\prime }\left(t\right)={\left[y^{\prime }\left(t\right),{\tilde{y}}^{\prime }\left(t\right),q^{\prime }\left(t\right),{\tilde{q}}^{\prime }\left(t\right)\right]}^T$ ， $P_0={\left[0,0,O\left({10}^{-3}\right),0\right]}^T$ ，则方程(5)与(6)可用以下向量形式表示：

$\{\begin{array}{c}{P}^{\prime }\left(t\right)=F\left(t,P\left(t\right)\right)\\ P\left(t_0\right)=P_0\end{array}$ (7)

上述(7)式对应的四阶龙格库塔向量形式为

$P_{n+1}=P_n+\frac{\Delta t}{6}\left(K_1+2K_2+2K_3+K_4\right)$ ， $n=0,1,2,\cdots$

$\{\begin{array}{c}{K}_1=F\left(t_n,P_n\right)\\ K_2=F\left(t_n+\frac{\Delta t}{2},P_n+\frac{\Delta t}{2}{K}_1\right)\\ K_3=F\left(t_n+\frac{\Delta t}{2},P_n+\frac{\Delta t}{2}{K}_2\right)\\ K_4=F\left(t_n+\Delta t,P_n+\Delta t{K}_3\right)\end{array}$

以上 $K_i$ $\left(i=1,2,3,4\right)$ 为向量，时间步长 $\Delta t$ 和 $t_n\left(n=0,1,\cdots ,N\right)$ 均为实数，截断误差为 $O\left(\Delta t^4\right).$

3.2. 有限差分法

有限差分方法作为一种经典的数值计算方法，广泛应用于微分方程组的求解。对(3)式中的 $y\left(t\right)$ 与 $q\left(t\right)$ 微分采用标准中心差商来代替，于是得到

$\{\begin{array}{c}\begin{array}{l}\frac{y\left(t_n+\Delta t\right)-2y\left(t_n\right)+y\left(t_n-\Delta t\right)}{\Delta t^2}+\left(2\delta \xi +\frac{\gamma }{\mu }\right)\left(\frac{y\left(t_n+\Delta t\right)-y\left(t_n-\Delta t\right)}{2\Delta t}\right)+{\delta }^{2}y\left(t_n\right)-Mq\left(t_n\right)\\ =\left[\frac{d^{2}y\left(t_n\right)}{d{t}^2}+\left(2\delta \xi +\frac{\gamma }{\mu }\right)\frac{dy\left(t_n\right)}{dt}+{\delta }^{2}y\left(t_n\right)-Mq\left(t_n\right)\right]+\left(2\delta \xi +\frac{\gamma }{\mu }\right)\frac{\Delta t^2}{6}{y}^{\left(3\right)}\left({\xi }_1\right)+\frac{\Delta t^2}{12}{y}^{\left(4\right)}\left({\xi }_3\right)\end{array}\\ \begin{array}{c}\frac{q\left(t_n+\Delta t\right)-2q\left(t_n\right)+q\left(t_n-\Delta t\right)}{\Delta t^2}+\epsilon (q^2\left(t_n\right)-1)\left(\frac{q\left(t_n+\Delta t\right)-q\left(t_n-\Delta t\right)}{2\Delta t}\right)+q\left(t_n\right)-A\frac{y\left(t_n+\Delta t\right)-2y\left(t_n\right)+y\left(t_n-\Delta t\right)}{\Delta t^2}\\ =\left[\frac{d^{2}q\left(t_n\right)}{d{t}^2}+\epsilon \left(q^2\left(t_n\right)-1\right)\frac{dq\left(t_n\right)}{dt}+q\left(t_n\right)-A\frac{d^{2}y\left(t_n\right)}{d{t}^2}\right]+\epsilon \left(q^2\left(t_n\right)-1\right)\frac{\Delta t^2}{6}{q}^{\left(3\right)}\left({\xi }_2\right)-A\frac{\Delta t^2}{12}{y}^{\left(4\right)}\left({\xi }_3\right)+\frac{\Delta t^2}{12}{q}^{\left(4\right)}(\; \xi \; 4\; )\end{array}\end{array}$

进而得到耦合方程(3)的差分方程为

$\{\begin{array}{c}\left(\frac{1}{\Delta t^2}+\frac{2\delta \xi +\gamma /\mu }{2\Delta t}\right)y_{n+1}\\ =M{q}_n+\left(\frac{2}{\Delta t^2}-{\delta }^2\right)y_n+\left(\frac{2\delta \xi +\gamma /\mu }{2\Delta t}-\frac{1}{\Delta t^2}\right)y_{n-1}\\ \left[\frac{1}{\Delta t^2}+\frac{\epsilon (q_n^2-1)}{2\Delta t}\right]q_{n+1}\\ =A\frac{y_{n+1}+y_{n-1}-2y_n}{\Delta t^2}+\left(\frac{2}{\Delta t^2}-1\right)q_n-\frac{1}{\Delta t^2}{q}_{n-1}+\epsilon \left(q_n^2-1\right)\frac{q_{n-1}}{2\Delta t}\end{array}$ (8)

其中， $y_n$ 、 $q_n$ 为函数 $y\left(t\right)、q\left(t\right)$ 在节点 $t_n$ 处的近似值，此格式的截断误差为 $O\left(\Delta t^2\right).$

同理，对初始条件(4)式使用标准中心差商可得

$\{\begin{array}{c}{y}_0=0,y_{-1}=y_1\\ q_0=O\left({10}^{-3}\right),q_{-1}=q_1\end{array}$

4. 数值算例

验证龙格–库塔法与有限差分法对刚性圆柱体涡激振动模型数值求解的准确性，采用文献 [23] 中的深海立管参数，即 $\xi =3.1\times {10}^{-3}$ ， $\mu$ = 250，A = 12， $\epsilon =0.3$ ， $C_D=2.0$ ，St = 0.2， $C_{L0}=1$ ， $Ur=5$ ，将龙格–库塔法与有限差分法的数值结果与其他文献中的数据进行对比。龙格–库塔法与有限差分法的计算时间步长均取 $\Delta t=0.1$ ，计算的无量纲时间区间为 $\left[0,2000\right]$ ，通过Matlab实现上述数值算法。

4.1. 无量纲位移

图2是刚性涡激振动耦合方程利用龙格库塔法与有限差分法求解，在运行 $t=2000$ 个无量纲时间时求得的立管无量纲位移时程图。从图像上看，当 $t$ 较小时， $y$ 先呈现缓慢上升的趋势，随后趋于稳定振动的状态。龙格库塔法与有限差分法计算得到的稳定振动状态临界点分别为 $t=1350$ 与 $t=1250$ ，无量纲位移最大值分别为0.054007、0.054131，与文献 [23] 数值结果非常吻合。

Figure 2. Variation of the dimensionless displacement y with dimensionless time t

图2. 无量纲位移y随无量纲时间t的变化

4.2. 无量纲振子

图3是刚性涡激振动耦合方程利用龙格库塔法与有限差分法求解，在运行 $t=2000$ 个无量纲时间时求得的立管无量纲振子时程图。从图像上看， $q$ 的振幅迅速增大，然后很快保持稳定振动。龙格库塔法与有限差分法计算得到的无量纲位移最大值分别为2.674921、2.676359，与文献 [23] 中数值结果非常吻合。

Figure 3. Variation of the dimensionless oscillator q with dimensionless time t

图3. 无量纲振子q随无量纲时间t的变化

5. 结论

在本文中，我们采用了高精度的四阶龙格库塔法和有限差分法两种数值方法对刚性圆柱体结构振子与尾流振子耦合的常微分方程组进行了数值模拟，探索了这些方法在模拟刚性圆柱体涡激振动响应特性时的有效性与准确性。通过数值算例结果表明，无论是采用四阶龙格库塔法还是有限差分法，在较长时间区间内计算出的无量纲位移与无量纲振子振幅均与现有文献中的结果保持了较高的一致性，证实了所采用的数值方法能够有效地模拟刚性圆柱体的涡激振动响应，而且为未来在涡激振动耦合的偏微分方程组数值模拟及相关参数研究方面提供了可靠的数值工具。在比较这两种方法时，四阶龙格–库塔法虽然在计算过程中需要更大的计算量，但它提供了更高的精度。这对于需要非常精确结果的应用场景是一个重要的优势。然而，有限差分法在处理具有复杂几何或边界条件的问题时更为方便，并且对于初步分析和快速求解具有一定的吸引力。最终，在考虑实际问题时，研究人员可以根据具体的应用需求和可用资源，选择适当的数值方法以实现最佳的模拟效果。

基金项目

中国石油大学(北京)油气资源与工程全国重点实验室课题资助(编号PRE/DX-2404)。