刚性立管顺流向与横向涡激振动的数值模拟分析

Numerical Simulation and Analysis of Vortex-Induced Vibration of Rigid Risers in the In-Line and Cross-Flow Directions

摘要: 涡激振动是深海立管在复杂海洋环境中常见的动力响应之一，常导致立管结构疲劳损坏。本研究针对刚性立管的涡激振动问题展开分析，重点讨论顺流向与横向振动特性。在数值模拟过程中，构建耦合常微分方程组来描述结构与尾流振子之间的相互作用，并采用四阶龙格–库塔法和有限差分法对其进行求解。模拟结果显示，两种数值方法所得结果高度一致，验证了所用方法的准确性与有效性。进一步地，基于该方法研究了一系列质量比与结构阻尼比条件下，位移振幅随折合速度的变化规律。

Abstract: Vortex-induced vibration (VIV) is one of the common dynamic responses of deep-sea risers in complex ocean environments, often leading to fatigue damage of the riser structure. This study analyzes the VIV of rigid risers, with a focus on the in-line and cross-flow vibration characteristics. In the numerical simulation process, a coupled system of ordinary differential equations is constructed to describe the interaction between the structure and the wake oscillator. The fourth-order Runge-Kutta method and the finite difference method are employed to solve the system. The simulation results show a high degree of consistency between the two numerical methods, verifying the accuracy and effectiveness of the proposed approach. Furthermore, based on this method, the variation of displacement amplitude with reduced velocity is investigated under a series of mass ratio and structural damping ratio conditions.

1. 引言

海洋工程结构在设计和操作中面临众多挑战，尤其是位于深水环境中的立管系统[1]。当流体以一定速度流过立管时，会在其两侧形成规律的涡旋脱落，这种脱落过程会产生周期性变化的力，进而导致立管产生振动，这一现象称作涡激振动[2]-[4]。这种振动会导致立管结构出现疲劳损伤，缩短使用寿命，甚至造成断裂事故。例如，在深海油气开发中，立管的疲劳失效可能导致油气泄漏，造成巨大的经济损失和环境污染。因此，准确预测刚性立管在涡激振动下的动力响应特性，对于优化立管设计、提高其抗疲劳性能具有重要意义。

涡激振动响应的研究方法可分为实验方法[5]-[8]和数值模拟方法[9]-[13]，数值模拟方法又可分为计算流体动力学方法以及半经验模型方法。深海立管涡激振动的流固耦合问题本质上是高阶微分耦合方程数值模拟问题，学者们采用的数值解法主要为有限元方法[14]-[16]与有限差分法[17]-[19]。

目前，刚性圆柱体涡激振动研究主要集中于单自由度横流(Cross-flow, CF)方向的振动特性，但顺流(In-line, IL)方向的振动特性同样不可忽视。本文基于尾流振子模型，研究双自由度涡激振动，探究CF与IL两个方向上的涡激振动响应特性。构建耦合常微分方程组，并通过四阶龙格-库塔法和有限差分法这两种方法进行数值求解，对结果进行对比分析。进一步地，研究了一系列质量比与结构阻尼比条件下，位移振幅随折合速度的变化规律，为深海立管设计提供更全面的参考。

2. 数学模型

本文以一根单位长度、直径为$D$ 的刚性圆柱体为研究对象，其在均匀流速$U$ 作用下同时发生横流向与顺流向的涡激振动。为建立数学模型，定义如下坐标系：以海流方向为$x$ 轴，表示顺流振动方向；与海流垂直的方向为$y$ 轴，表示横流振动方向；沿立管重力方向为$z$ 轴，表示铅直方向，如图1所示。

设$X$ 和Y分别表示圆柱体在顺流$x$ 方向和横流$y$ 方向的位移，则其结构振动方程为[20]：

$\{\begin{array}{c}m\frac{d^{2}X(T)}{d{T}^2}+r\frac{dX\left(T\right)}{dT}+h[X\left(T\right)+{\alpha }_x^{\ast }{X}^3\left(T\right)+{\beta }_x^{\ast }X\left(T\right)\cdot Y^2\left(T\right)]=F_x\left(T\right),\\ m\frac{d^{2}Y\left(T\right)}{d{T}^2}+r\frac{dY\left(T\right)}{dT}+h[Y\left(T\right)+{\alpha }_x^{\ast }{Y}^3\left(T\right)+{\beta }_x^{\ast }Y\left(T\right)\cdot X^2\left(T\right)]=F_y\left(T\right),\end{array}$

上式中，$T$ 为时间；$m$ 为系统总质量；$r$ 为阻尼系数，包括结构阻尼系数$r_s$ 与流体阻尼系数$r_f$ 两项，其中$r_f={\omega }_f\gamma \rho D^2$ ，${\omega }_f$ 为漩涡脱落频率，$\gamma$ 为粘滞力系数，取0.8，$\rho$ 为流体密度。$h$ 为系统刚度，${\alpha }_x^{\ast }、{\alpha }_y^{\ast }、{\beta }_x^{\ast }、{\beta }_y^{\ast }$ 为结构的几何非线性系数。

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

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

$F_x、F_y$ 分别表示作用在结构上的$x$ 方向以及$y$ 方向水动力载荷，表达式如下：

$F_x\left(T\right)=\overline{F_D}\left(T\right)\text{+}{F}_D\left(T\right)-\frac{F_L\left(T\right)}{U}\cdot \frac{dY\left(T\right)}{dT},$

$F_y\left(T\right)\text{=}{F}_L\left(T\right)+\frac{F_D\left(T\right)}{U}\cdot \frac{dY\left(T\right)}{dT},$

$\{\begin{array}{l}\overline{F_D}\left(T\right)=\frac{1}{2}\rho D{U}^2\overline{C_D}\left(T\right),\\ F_D\left(T\right)=\frac{1}{2}\rho D{U}^2{C}_D\left(T\right),\\ F_L\left(T\right)=\frac{1}{2}\rho D{U}^2{C}_L\left(T\right),\end{array}$

上式中，$\overline{C_D}、C_D、C_L$ 分别表示平均拖曳力系数、振荡拖曳力系数以及升力系数。$\overline{C_D}\text{=1}\text{.2}$，$C_D\left(T\right)\text{=}{C}_{D\text{0}}\cdot p\left(T\right)/\text{2}$ ，$C_L\left(T\right)\text{=}{C}_{\text{L0}}\cdot q\left(T\right)/\text{2}$ ，取$C_{D\text{0}}=0.2$ ，$C_{\text{L0}}=0.3$ 。

为了描述尾流振子的非线性特性，本研究采用改进的Van der Pol方程[21]，其数学表达式如下：

$\{\begin{array}{c}\frac{d^{2}p\left(T\right)}{d{T}^2}+2{\epsilon }_x{\omega }_f\left[p^2\left(T\right)-1\right]\frac{dp\left(T\right)}{dT}+4{\omega }_f^{2}p\left(T\right)=\frac{A_x}{D}\frac{d^{2}X\left(T\right)}{d{T}^2},\\ \frac{d^{2}q\left(T\right)}{d{T}^2}+{\epsilon }_y{\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_y}{D}\frac{d^{2}Y\left(T\right)}{d{T}^2},\end{array}$

上式中，${\epsilon }_x、{\epsilon }_y、A_x、A_y$ 为经验参数。

因此可得涡激振动的流固耦合方程为：

$\{\begin{array}{c}m\frac{d^{2}X\left(T\right)}{d{T}^2}+r\frac{dX\left(T\right)}{dT}+h\left[X\left(T\right)+{\alpha }_x^{\ast }{X}^3\left(T\right)+{\beta }_x^{\ast }X\left(T\right)\cdot Y^2\left(T\right)\right]\\ =\frac{\rho U^{2}D\overline{C_D}}{2}+\frac{\rho U^{2}D{C}_{D0}p\left(T\right)}{4}-\frac{\rho UD{C}_{L0}q\left(T\right)}{4}\cdot \frac{dY\left(T\right)}{dT},\\ m\frac{d^{2}Y\left(T\right)}{d{T}^2}+r\frac{dY\left(T\right)}{dT}+h[Y\left(T\right)+{\alpha }_y^{\ast }{Y}^3\left(T\right)+{\beta }_y^{\ast }Y\left(T\right)\cdot X^2\left(T\right)]\\ =\frac{\rho U^{2}D{C}_{L0}q\left(T\right)}{4}+\frac{\rho UD{C}_{D0}p\left(T\right)}{4}\cdot \frac{dY\left(T\right)}{dT},\\ \frac{d^{2}p\left(T\right)}{d{T}^2}+2{\epsilon }_x{\omega }_f\left[p^2\left(T\right)-1\right]\frac{dp\left(T\right)}{dT}+4{\omega }_f^{2}p\left(T\right)=\frac{A_x}{D}\frac{d^{2}X\left(T\right)}{d{T}^2},\\ \frac{d^{2}q\left(T\right)}{d{T}^2}+{\epsilon }_y{\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_y}{D}\frac{d^{2}Y\left(T\right)}{d{T}^2}.\end{array}$ (1)

立管的初始条件为：

$\{\begin{array}{l}{X|}_{T=T_0}=0,\begin{array}{cc}& \end{array}\begin{array}{c}\end{array}{\frac{dX}{dT}|}_{T=T_0}=0,\\ {Y|}_{T=T_0}=0,\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}{\frac{dY}{dT}|}_{T=T_0}=0,\\ {p|}_{T=T_0}=O\left({10}^{-3}\right),\begin{array}{c}\end{array}{\frac{dp}{dT}|}_{T=T_0}=0,\\ {q|}_{T=T_0}=O\left({10}^{-3}\right),\begin{array}{c}\end{array}{\frac{dq}{dT}|}_{T=T_0}=0.\end{array}$ (2)

3. 求解方法

首先对原方程组进行无量纲化处理，引入以下无量纲变量：

$x=X/D,\begin{array}{c}\end{array}y=Y/D,\begin{array}{c}\end{array}t=T\cdot {\omega }_n,$

其中，${\omega }_n$ 为系统固有频率，其数学表达式为${\omega }_n=\sqrt{h/m}$ 。

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

$\{\begin{array}{l}\frac{{\text{d}}^{2}x\left(t\right)}{{\text{dt}}^2}+\left(2\xi +\frac{\gamma \Omega }{\mu }\right)\frac{\text{d}x\left(t\right)}{\text{dt}}+\left[x\left(t\right)+{\alpha }_x^{}{x}^3\left(t\right)+{\beta }_x^{}x\left(t\right)\cdot y^2\left(t\right)\right]\\ =\overline{M_D}{\Omega }^2+M_D{\Omega }^2\cdot p\left(t\right)-M_L{\Omega }^2\cdot q\left(t\right)\frac{2\pi }{Ur}\cdot \frac{dy}{dt},\\ \frac{{\text{d}}^{2}y\left(t\right)}{{\text{dt}}^2}+\left(2\xi +\frac{\gamma \Omega }{\mu }\right)\frac{\text{d}y\left(t\right)}{dt}+[y\left(t\right)+{\alpha }_y^{}{y}^3\left(t\right)+{\beta }_y^{}y\left(t\right)\cdot x^2\left(t\right)]\\ =M_L{\Omega }^2\cdot q\left(t\right)+M_D{\Omega }^2\cdot p\left(t\right)\frac{2\pi }{Ur}\cdot \frac{dy}{dt},\\ \frac{{\text{d}}^{2}p\left(t\right)}{{\text{dt}}^2}+2{\epsilon }_x\Omega \left[p^2\left(t\right)-1\right]\frac{\text{d}p(t)}{dt}+4{\Omega }^{2}p\left(t\right)=A_x\frac{{\text{d}}^{2}x\left(t\right)}{{\text{dt}}^2},\\ \frac{{\text{d}}^{2}q\left(t\right)}{{\text{dt}}^2}+{\epsilon }_y\Omega \left[q^2\left(t\right)-1\right]\frac{\text{d}q(t)}{dt}+{\Omega }^{2}q\left(t\right)=A_y\frac{{\text{d}}^{2}y\left(t\right)}{{\text{dt}}^2},\end{array}$ (3)

上式中，$\xi$ 表示结构阻尼比，表示为：$\xi =r_s/2m{\omega }_n$ ；${\alpha }_x^{}、{\alpha }_y^{}、{\beta }_x^{}、{\beta }_y^{}$ 表示无量纲几何非线性系数；$\Omega$ 表示频率比，表示为：$\Omega ={\omega }_f/{\omega }_n=St\cdot Ur$ ，$St$ 为斯特劳哈尔数，$St=\text{0}\text{.2}$ ，$Ur$ 表示折合速度，表示为：$Ur=2\pi U/{\omega }_{n}D$ ；$\mu$ 表示质量比，表示为：$\mu =m/\rho D^2$ 。${\overline{M}}_D、M_D、M_L$ 表示系统无量纲参数，表达式如下：

$\{\begin{array}{c}{\overline{M}}_D=\frac{{\overline{C}}_D}{8{\pi }^{2}S{t}^2\mu },\\ M_D=\frac{C_{D0}}{16{\pi }^{2}S{t}^2\mu },\\ M_L=\frac{C_{L0}}{16{\pi }^{2}S{t}^2\mu }.\end{array}$

初始条件(2)式转化为：

$\{\begin{array}{l}{x|}_{t=t_0}=0,\begin{array}{cc}& \end{array}\begin{array}{c}\end{array}{\frac{dx}{dt}|}_{t=t_0}=0,\\ {y|}_{t=t_0}=0,\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}{\frac{dy}{dt}|}_{t=t_0}=0,\\ {p|}_{t=t_0}=O\left({10}^{-3}\right),\begin{array}{c}\end{array}{\frac{dp}{dt}|}_{t=t_0}=0,\\ {q|}_{t=t_0}=O\left({10}^{-3}\right),\begin{array}{c}\end{array}{\frac{dq}{dt}|}_{t=t_0}=0.\end{array}$ (4)

假设$x\left(t\right)$ 、$y\left(t\right)$ 、$p\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{x}\left(t\right)=\frac{dx\left(t\right)}{dt}$ ，$\tilde{y}\left(t\right)=\frac{dy\left(t\right)}{dt}$ ，$\tilde{p}\left(t\right)=\frac{dp\left(t\right)}{dt}$ ，$\tilde{q}\left(t\right)=\frac{dq\left(t\right)}{dt}$ ，于是(3)可化为如下一阶微分方程组：

$\{\begin{array}{c}\frac{dx(t)}{dt}=\tilde{x}(t),\\ \frac{d\tilde{x}(t)}{dt}=-\left(2\xi +\frac{\gamma \Omega }{\mu }\right)\tilde{x}(t)-[x\left(t\right)+{\alpha }_x^{}{x}^3\left(t\right)+{\beta }_x^{}x\left(t\right)y^2\left(t\right)]+{\overline{M}}_D{\Omega }^2+M_D{\Omega }^{2}p(t)-M_L{\Omega }^{2}q(t)\frac{2\pi }{Ur}\tilde{y}(t),\\ \frac{dy\left(t\right)}{dt}=\tilde{y}\left(t\right),\\ \frac{d\tilde{y}\left(t\right)}{dt}=-\left(2\xi +\frac{\gamma \Omega }{\mu }\right)\tilde{y}\left(t\right)-[y\left(t\right)+{\alpha }_y^{}{y}^3\left(t\right)+{\beta }_y^{}y\left(t\right)x^2\left(t\right)]+M_L{\Omega }^{2}q\left(t\right)+M_D{\Omega }^{2}p\left(t\right)\frac{2\pi }{Ur}\tilde{y}\left(t\right),\\ \frac{dp\left(t\right)}{dt}=\tilde{p}\left(t\right),\\ \frac{d\tilde{p}\left(t\right)}{dt}=-2{\epsilon }_x\Omega \left[p^2\left(t\right)-1\right]\tilde{p}(t)+4{\Omega }_{}^{2}p\left(t\right)+A_x\frac{d\tilde{x}\left(t\right)}{dt},\\ \frac{dq\left(t\right)}{dt}=\tilde{q}\left(t\right),\\ \frac{d\tilde{q}\left(t\right)}{dt}=-{\epsilon }_y\Omega \left[q^2\left(t\right)-1\right]\tilde{q}\left(t\right)-{\Omega }^{2}q\left(t\right)+A_y\frac{d\tilde{y}\left(t\right)}{dt}.\end{array}$ (5)

初始条件(4)式化为：

$\{\begin{array}{c}x\left(t_0\right)=0,\\ \tilde{x}\left(t_0\right)=0,\\ y\left(t_0\right)=0,\\ \tilde{y}\left(t_0\right)=0,\\ p\left(t_0\right)=O\left({10}^{-3}\right),\\ \tilde{p}\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[x\left(t\right),\tilde{x}\left(t\right),y\left(t\right),\tilde{y}\left(t\right),p\left(t\right),\tilde{p}\left(t\right),q\left(t\right),\tilde{q}\left(t\right)\right]}^T,$

$P^{\prime }\left(t\right)={\left[x^{\prime }\left(t\right),{\tilde{x}}^{\prime }\left(t\right),y^{\prime }\left(t\right),{\tilde{y}}^{\prime }\left(t\right),p^{\prime }\left(t\right),{\tilde{p}}^{\prime }\left(t\right),q^{\prime }\left(t\right),{\tilde{q}}^{\prime }\left(t\right)\right]}^T,$

$P_0={\left[0,0,\text{0},0,O\left({10}^{-3}\right)，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)式中的$x\left(t\right)$ 、$y\left(t\right)$ 、$p\left(t\right)$ 与$q\left(t\right)$ 可以采用标准中心差商进行离散化近似，定义以下中心差商：

$\begin{array}{l}{\frac{dx}{dt}|}_{t=t_n}=\frac{x\left(t_n+\Delta t\right)-x\left(t_n-\Delta t\right)}{2\Delta t}+O\left(\Delta t^2\right),{\frac{d^{2}x}{d{t}^2}|}_{t=t_n}=\frac{x\left(t_n+\Delta t\right)-2x\left(t_n\right)+x\left(t_n-\Delta t\right)}{\Delta t^2}+O\left(\Delta t^2\right),\\ {\frac{dy}{dt}|}_{t=t_n}=\frac{y\left(t_n+\Delta t\right)-y\left(t_n-\Delta t\right)}{2\Delta t}+O\left(\Delta t^2\right),{\frac{d^{2}y}{d{t}^2}|}_{t=t_n}=\frac{y\left(t_n+\Delta t\right)-2y\left(t_n\right)+y\left(t_n-\Delta t\right)}{\Delta t^2}+O\left(\Delta t^2\right),\end{array}$

$\begin{array}{l}{\frac{dp}{dt}|}_{t=t_n}=\frac{p\left(t_n+\Delta t\right)-p\left(t_n-\Delta t\right)}{2\Delta t}+O\left(\Delta t^2\right),{\frac{d^{2}p}{d{t}^2}|}_{t=t_n}=\frac{p\left(t_n+\Delta t\right)-2p\left(t_n\right)+p\left(t_n-\Delta t\right)}{\Delta t^2}+O\left(\Delta t^2\right),\\ {\frac{dq}{dt}|}_{t=t_n}=\frac{q\left(t_n+\Delta t\right)-q\left(t_n-\Delta t\right)}{2\Delta t}+O\left(\Delta t^2\right),{\frac{d^{2}q}{d{t}^2}|}_{t=t_n}=\frac{q\left(t_n+\Delta t\right)-2q\left(t_n\right)+q\left(t_n-\Delta t\right)}{\Delta t^2}+O\left(\Delta t^2\right).\end{array}$

于是耦合方程(3)的差分方程为：

$\{\begin{array}{c}\left(\frac{1}{\Delta t^2}+\frac{2\xi +\gamma \Omega /\mu }{2\Delta t}\right)x_{n+1}\\ =\left(\frac{2}{\Delta t^2}-1\right)x_n+\left(\frac{2\xi +\gamma \Omega /\mu }{2\Delta t}-\frac{1}{\Delta t^2}\right)x_{n-1}+{\overline{M}}_D{\Omega }^2+M_D{\Omega }^2{p}_n-M_L{\Omega }^2{q}_n\frac{2\pi }{U_r}\frac{y_{n+1}-y_{n-1}}{2\Delta t}-{\alpha }_x{x}_n^3-{\beta }_x{x}_n{y}_n^2，\\ \left(\frac{1}{\Delta t^2}+\frac{2\xi +\gamma \Omega /\mu }{2\Delta t}-M_D{\Omega }^2{p}_n\frac{2\pi }{U_r}\frac{1}{2\Delta t}\right)y_{n+1}\\ =\left(\frac{2}{\Delta t^2}-1\right)y_n+\left(\frac{2\xi +\gamma \Omega /\mu }{2\Delta t}-M_D{\Omega }^2{p}_n\frac{2\pi }{U_r}\frac{1}{2\Delta t}-\frac{1}{\Delta t^2}\right)y_{n-1}+M_L{\Omega }^2{q}_n-{\alpha }_y{y}_n^3-{\beta }_y{y}_n{x}_n^2，\\ \left[\frac{1}{\Delta t^2}+\frac{2{\epsilon }_x\Omega \left(p_n^2-1\right)}{2\Delta t}\right]p_{n+1}=A_x\frac{x_{n+1}+x_{n-1}-2x_n}{\Delta t^2}+\left(\frac{2}{\Delta t^2}-4{\Omega }^2\right)p_n-\frac{1}{\Delta t^2}{p}_{n-1}+2{\epsilon }_x\Omega \left(p_n^2-1\right)\frac{p_{n-1}}{2\Delta t},\\ \left[\frac{1}{\Delta t^2}+\frac{{\epsilon }_y\Omega \left(q_n^2-1\right)}{2\Delta t}\right]q_{n+1}=A_y\frac{y_{n+1}+y_{n-1}-2y_n}{\Delta t^2}+\left(\frac{2}{\Delta t^2}-{\Omega }^2\right)q_n-\frac{1}{\Delta t^2}{q}_{n-1}+{\epsilon }_y\Omega \left(q_n^2-1\right)\frac{q_{n-1}}{2\Delta t}.\end{array}$ (8)

以上$x_n$ 、$y_n$ 、$p_n$ 与$q_n$ 分别表示$x\left(t\right)$ 、$y\left(t\right)$ 、$p\left(t\right)$ 与$q\left(t\right)$ 在t = t_n处的近似值，截断误差为$O\left(\Delta t^2\right)$ 。

对于初始条件(4)式，同样可以采用标准中心差商进行离散化处理：

$\{\begin{array}{c}\begin{array}{l}{x}_0=0,x_{-1}=x_1,\\ y_0=0,y_{-1}=y_1,\\ p_0=O\left({10}^{-3}\right),p_{-1}=p_1,\end{array}\\ q_0=O\left({10}^{-3}\right),q_{-1}=q_1.\end{array}$

3.3. 数值验证

为了验证数值模拟，本文采用四阶龙格–库塔法和有限差分法分别对耦合方程进行求解，并将模拟结果进行对比分析。立管参数来自文献[20]，即$m^{\ast }=2.36$ ，$U_r=6$ ，$\xi =6\times {10}^{-3}$ ，$\gamma =0.8$ ， ${\alpha }_x^{}={\alpha }_y^{}={\beta }_x^{}={\beta }_y^{}=0.7$ ，$A_x=A_y=12$ ，$St=0.2$ ，${\epsilon }_x=0.3$ 。根据参考文献[5]，得${\epsilon }_y=0.00234e^{0.2283m*}$ 。数值模拟采用固定时间步长$\Delta t=0.1$ ，无量纲时间取$\left[0,1000\right]$ 。

图2和图3分别展示了刚性涡激振动耦合方程采用龙格–库塔法和有限差分法求解后，在运行$t=1000$ 个无量纲时间时立管无量纲位移的时程曲线。从图中可以看出，当无量纲时间$t$ 较小时，$x$ 先呈现明显的振荡特性，随后振幅逐渐衰减并趋于稳定状态。$y$ 表现为初始阶段的上升趋势，之后同样趋于稳定振动。通过对比两种数值方法的结果，龙格–库塔法计算得到的$x$ 和$y$ 最大值分别为0.3085和1.3075，而有限差分法得到的相应值分别为0.3085和1.3083。两种方法的结果高度吻合，进一步验证了数值模拟的准确性和可靠性。

图4和图5分别呈现了刚性涡激振动耦合方程采用龙格–库塔法和有限差分法求解后，在运行$t=1000$ 个无量纲时间时立管无量纲振子的时程曲线。从图中可以看出，振子$p$ 的振幅表现出先增长后趋于平缓，随后再次增长并最终达到稳定状态；而振子$q$ 的振幅则呈现快速增大并迅速趋于稳定的特征。通过数值计算，龙格–库塔法得到的$p$ 和$q$ 最大值分别为3.1367和24.7851，而有限差分法得到的相应值分别为3.1547和24.7966。两种数值方法的结果高度一致，进一步验证了数值模拟的准确性和可靠性，表明所采用的数值方法能够有效捕捉涡激振动的非线性动力学特性。

(a) 四阶龙格–库塔法 (b) 有限差分法

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

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

(a) 四阶龙格–库塔法 (b) 有限差分法

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

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

(a) 四阶龙格–库塔法 (b) 有限差分法

Figure 4. Variation of the dimensionless oscillator p with dimensionless time t

图4. 无量纲振子p随无量纲时间t的变化

(a) 四阶龙格–库塔法 (b) 有限差分法

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

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

4. 数值算例

4.1. 主要参数介绍

4.1.1. 质量比与结构阻尼比

在涡激振动的精确预测和分析中，质量比和结构阻尼比是与结构物相关的两个主要参数。在分析中大多时使用的质量比为$m^{\ast }$ ，其与前文所述$\mu$ 之间的关系为$m^{\ast }=4\mu /\pi -C_M$ ，4，其中$C_M$ 为附加质量系数，对于圆柱体通常取为1。其含义为结构质量与被结构排开的液体质量的比。本文质量比取1到10之间。

结构阻尼比$\xi$ 表达式为：$\xi =r_s/2m{\omega }_n$ ，$r_s$ 为结构阻尼系数，可以通过振幅衰减实验来测量，${\omega }_n$ 为振动系统的固有频率。结构阻尼比用于衡量结构系统在振动过程中能量耗散的能力，通常用于描述结构如何衰减振动能量。本文结构阻尼比取0.01、0.05、0.1和0.2来研究。

4.1.2. 响应振幅与折合速度

振幅大小是评估涡激振动对结构疲劳损伤程度的重要参数。本文研究了刚性立管在顺流向与横向上的运动，并计算了相应的响应振幅。记无量纲顺流向和横向振幅分别为：$\max \left(x\right)$ 、$\max \left(y\right)$ 。

在涡激振动研究中，折合速度$Ur$ 是一个重要的无量纲参数，其表达式为：$Ur=2\pi U/{\omega }_{n}D$ ，它用于表征流体流速与结构固有振动频率之间的关系，通常用于分析不同流速条件下的VIV响应特性。本文折合速度取值范围取1~14。

4.2. 不同质量比分析

本节采用龙格–库塔方法对刚性立管的涡激振动响应进行详细分析。通过数值模拟，探讨了质量比以及结构阻尼比对立管振动行为的影响。在模拟过程中，立管所用参数与文献[20]中所述一致，时间步长$\Delta t$ 保持为0.01。

在研究质量比时，保持其他参数不变，对五种不同的质量比$m^{\ast }$ 进行数值模拟，分别为2、4、6、8、10，对于折合速度$Ur$ ，选择1.0至14.0的范围，并以0.1的间隔进行变化。

图6清晰展示了在不同质量比(分别为2、4、6、8、10)条件下，圆柱体在各种折合速度$Ur$ 下，其无量纲位移振幅$\max \left(x\right)$ 与无量纲位移振幅$\max \left(y\right)$ 的变化规律。可以明显看出，当质量比逐渐变大时，最大无量纲位移振幅$\max \left(x\right)$ 与$\max \left(y\right)$ 均呈现出递减的趋势，并且达到最大位移振幅所对应的折合速度$Ur$ 也随之减小。这表明在较高的质量比下，圆柱体在较低的流速下即可达到最大的振动响应。

Figure 6. Variation of dimensionless displacement amplitudes max(x) and max(y) with reduced velocity Ur under different mass ratios

图6. 不同质量比下，无量纲位移振幅max(x)、max(y)随折合速度Ur的变化情况

4.3. 不同结构阻尼比分析

进一步探讨阻尼比的作用，在固定质量比$m^{\ast }=2.36$ 的条件下，深入分析结构阻尼比$\xi$ 对圆柱体涡激振动响应的影响。选取了四组不同的结构阻尼比$\xi$ ，具体数值为0.01、0.05、0.1和0.2。

图7清晰展示了结构阻尼比分别取值为0.01、0.05、0.1和0.2时，刚性圆柱体的无量纲位移振幅$\max \left(x\right)$ 、$\max \left(y\right)$ 随折合速度$Ur$ 的变化情况。可以明显看出，随着阻尼比的逐步增加，刚性圆柱体的涡激振动所表现出的最大无量纲位移幅值$\max \left(x\right)$ 、$\max \left(y\right)$ 均呈现出递减的趋势。

Figure 7. Variation of dimensionless displacement amplitudes max(x) and max(y) with reduced velocity Ur under different structural damping ratios

图7. 不同结构阻尼比下，无量纲位移振幅max(x)、max(y)随折合速度Ur的变化情况

结果表明，增加结构阻尼比$\xi$ 能够有效降低圆柱体在流体动力作用下的振动响应。其物理机制主要包括：能量耗散与抑制共振。增加结构阻尼比会提高系统的能量耗散能力，使振动能量更快地转化为热能或其他形式的能量；同时在涡激振动中，增加结构阻尼比可以显著降低共振幅值，避免结构因共振而导致的疲劳损伤。从数学角度看，增加结构阻尼比$\xi$ 会增大运动方程中的$\left(2\xi +\frac{\gamma \Omega }{\mu }\right)\frac{\text{d}x\left(t\right)}{\text{dt}}$ 一项，从而改变系统的动力学特性。

5. 结论

本研究针对刚性立管的涡激振动问题展开了系统性研究，重点分析了顺流向和横向的振动响应特性。为了准确描述刚性圆柱体结构振子与尾流振子的耦合动力学行为，本文建立了相应的常微分方程组，并采用两种高精度数值方法——四阶龙格–库塔法和有限差分法进行求解。数值模拟结果表明，两种方法所得结果具有高度一致性，验证了数值方法的可靠性和计算精度。此外，本文进一步探讨了质量比和结构阻尼比对涡激振动特性的影响规律，通过系统参数分析揭示了位移振幅随折合速度的变化趋势。研究结果表明，质量比和结构阻尼比是影响立管涡激振动响应的关键参数，其变化会显著改变振动幅值和锁定区间范围。这一发现为工程实践中立管涡激振动的预测与控制提供了理论依据。

基金项目

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