0 引　言

在海洋油气平台中，立管是连接海底管道和平台生产设备的重要设备，工作过程中受到海流、海浪和海风等各类环境载荷的作用。当海洋载荷的激励频率接近于管道自振频率时，海洋立管会诱发耦合共振，加剧立管疲劳破坏，可能造成严重的经济损失。因此，海洋立管流体动力学特性及流致振动等研究对深海油气设备开发有重要研究意义。

目前，对涡激振动（vortex-induced vibration，VIV）现象进行了深入研究，并提出多种振动抑制技术来克服海洋管道的涡激振动^[1]。当海洋立管工作时，海洋波流绕过立管并在管后分离，产生交替脱落的漩涡，使管道承受周期性的流体作用力。当流体作用力的驱动频率（或漩涡脱落频率）接近于管道的固有频率时，海洋立管的横向振幅就会明显增大，出现涡激振动现象。该现象也被称为“频率锁定”现象，即管道振动会迫使漩涡脱落频率固定在立管固有频率附近，而不按漩涡本身的脱落频率脱落^[2]。Iranpour等^[3]对立管涡激振动进行了试验，发现雷诺数在2000～6000范围内立管存在4种不同形式的漩涡脱落形式。Sarpkaya^[4]综述了涡激振动对柔性立管的影响，可利用管道受迫振动实验来分析自激振动特性。Blevins^[5]发现当流体约化速度在6～7.5范围内，立管容易发生频率锁定现象，且载荷频率接近管道固有频率时，振动幅值会明显增大。娄敏等^[6]实验研究了海洋立管内外流对涡激振动的共同影响。

随着计算机软件技术发展，涡激振动现象数值模拟的研究取得很大进展。白长旭等^[7]利用Morison方程求解流体作用力，运用Abaqus软件对海洋立管进行了动力学分析，发现随海流流速增加，立管振动的位移先增大后减小。Mukundan等^[8]系统研究了海洋立管的涡激振动理论，并对立管疲劳损伤进行了评估分析。Hartlen^[9]分析了弹性支承圆柱的横向涡激振动特性，并建立了“尾流振子”经验模型。任大鹏^[10]利用Fluent软件建立了立管涡激振动模型，研究了管道在均匀来流作用下横向振动的周期。赵鹏良^[11]采用Fluent动网格技术，对立管低质量比和小阻尼等情况进行了计算，成功捕捉到“锁定”、“失谐”等现象。上述研究均采用现有软件计算，虽然能够定量研究海洋立管的涡激振动特性，但存在无法较好耦合海洋波浪以及计算精度差等问题。

在海洋波流载荷作用下，为准确计算海洋立管与流体的相互作用过程，在考虑立管与海洋波面相互影响的前提下，本文利用DNS和LES方法求解Navier-Stokes方程，利用IB法计算流体对海洋立柱的作用力，采用CLSVOF方法描述了海洋波浪的自由界面。与前期工作相比，本文将该方法应用到海洋立柱的数值模拟中，分析立柱的流固耦合特性及立柱对海洋波面的影响，为海洋立柱涡激振动特性的数值模拟提供高精度计算方法。

1 数值方法 1.1 Navier-Stokes方程流体求解器

以海洋立柱直径D和来流速度 $ {U_\infty } $ 进行无量纲化，不可压流体N-S控制方程为：

$ \frac{{\partial {u_{\text{i}}}}}{{\partial t}} + \frac{{\partial \left( {{u_{\text{i}}}{u_{\text{j}}}} \right)}}{{\partial {x_{\text{j}}}}} = - \frac{{\partial p}}{{\partial {x_{\text{j}}}}} + \frac{1}{{Re}}\frac{{{\partial ^2}{u_{\text{i}}}}}{{\partial {x_{\text{j}}}{x_{\text{j}}}}} + {f_{\text{i}}} ，$

(1)

$ \frac{{\partial {u_{\text{i}}}}}{{\partial {x_{\text{i}}}}} = 0 。$

(2)

式中： $ {x}_{i}\left(i=1，2，3\right) $ 代表3个坐标方向 $ x $ ， $ y $ 和 $ z $ ； $ {u_i} $ 表示速度 $ u $ ， $ v $ 和 $ w $ ；t为以 $ {D \mathord{\left/ {\vphantom {D {{U_\infty }}}} \right. } {{U_\infty }}} $ 作无量纲化的时间；p为以 $ \rho U_\infty ^2 $ 作无量纲化的压力， $ \rho $ 为流体密度；Re为雷诺数；f_i为立柱与流体的相互作用力。当流体雷诺数较低时，采用直接求解方法（DNS）进行数值求解；当流体雷诺数较高时，利用大涡模拟方法（LES）求解N-S方程，流体小尺度运动由SGS（samgorinsky eddy viscosity）模型进行等效分析，对应控制方程为：

$ \frac{{\partial {{\bar u}_{{i}}}}}{{\partial t}} + \frac{{\partial \left( {{{\bar u}_{{i}}}{{\bar u}_{{j}}}} \right)}}{{\partial {x_{{j}}}}} = - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_{{i}}}}} + {v_{{t}}}\frac{{{\partial ^2}{{\bar u}_{{i}}}}}{{\partial {x_{{j}}}\partial {x_{{j}}}}} - \frac{{\partial \tau _{{{ij}}}^d}}{{\partial {x_{{j}}}}} + {f_{{i}}} ，$

(3)

$ \frac{{\partial {{\bar u}_{\text{i}}}}}{{\partial {x_{\text{i}}}}} = 0。$

(4)

式中：横线表示过滤后的变量， $ {v_t} $ 为湍流涡流粘性； $ \tau _{i,j}^d $ 为亚格子应力， $ \tau _{i,j}^d = - 2{C_s}{\bar \Delta ^2}\left| {\bar S} \right|{\bar S_{ij}} $ ， $ {C_S} $ 为动态参数； $ \bar \Delta $ 为湍流尺度过滤， $ \bar \Delta = {\left( {\Delta x\Delta y\Delta z} \right)^{1/3}} $ ， $ \Delta x $ ， $ \Delta y $ 和 $ \Delta z $ 为3个方向上的网格尺寸； $ {\overline S _{ij}} $ 为应变率张量，幅值 $ \left| {\overline S } \right| = \sqrt {2{{\overline S }_{ij}}{{\overline S }_{ij}}} $ ， $ {\overline S _{ij}} = \dfrac{1}{2}\left( {\dfrac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \dfrac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) $ 。

N-S方程各变量分布在三维正交网格上，控制方程采用中心差分格式离散，由Runge-Kutta格式实现时间推进，其中压力泊松方程通过调用PETS库函数进行求解。本文所使用的N-S求解器在Linux系统下使用Fortran语言编写，采用MPI命令实现程序并行运算，具体算法详见文献[12]。本文计算所需硬件设备为国家超级计算长沙中心“天河”计算机，二维和三维算例分别采用64个和256个CPU进行并行计算。

1.2 CLSVOF方法

将空气和海洋流体看作具有不同密度和粘度的单个系统，其N-S控制方程为:

$ \frac{{\partial {\boldsymbol{u}}}}{{\partial {\text{t}}}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}} = \frac{1}{\rho }\left( { - \nabla p + \nabla \cdot \left( {\mu {\mathbf{S}}} \right) + \rho g + \sigma {\mathbf{\kappa }}\delta \left( {{x_s}} \right) - \nabla \cdot {{\mathbf{\tau }}_{sgs}}} \right) ，$

(5)

$ \nabla \cdot {\mathbf{u}} = 0。$

(6)

式中：u为速度， $ \rho $ 为密度， $ p $ 为压力， $ \mu $ 为动态粘度， $ S $ 为剪切张量， $ g $ 为重力加速度，σ为表面张力， $ {\mathbf{\kappa }} $ 为界面曲率， $ \delta \left( {{x_s}} \right) $ 为在多相流界面位置 $ {x_s} $ 上的平滑函数， $ {\tau _{{\text{sgs}}}} $ 为网格尺度应力，由LES方法Smagorinsky模型来动态决定。其中，密度 $ \rho $ 和动态粘度 $ \mu $ 是由取决于体积分数F来决定，表达式为：

$ \rho = {\rho _{\text{a}}}\left( {1 - F} \right) + {\rho _{\rm{w}}}F ，$

(7)

$ \mu = {\mu _{\text{a}}}\left( {1 - F} \right) + {\mu _{\rm{w}}}F 。$

(8)

式中：下标“a”和“w”分别代表空气和海水。在空气中F值为0；在水中F值为1；在空气-海水界面的网格单元中F值为0 ～ 1。本文采用CLSVOF方法捕获空气-海水界面，迭代流程如图2所示。Level-set函数 $ \phi $ 和VOF函数F每一步通过求解下列对流方程进行更新：

$ \frac{{\partial \phi }}{{\partial {{t}}}} + \nabla \cdot \left( {u\phi } \right) = 0 ，$

(9)

$ \frac{{\partial F}}{{\partial {{t}}}} + \nabla \cdot \left( {uF} \right) = 0 。$

(10)

在每个步骤中，利用Level-set函数对多相流界面进行重构，采用分段线性界面计算方法（piecewise-linear interface calculation，PLIC）进行计算，使F在重构界面上的通量得到精确解。当Level-set函数和VOF函数更新之后，F的函数值用来修正 $ \phi $ 值，以确保各相流体的质量守恒。

1.3 浸入边界方法

采用Level-set方法确定海洋管道界面，管道截面为半径为R₀的圆形，对应Level-set函数为：

$ \varphi \left( x \right) = \sqrt {{{\left( {{{x - }}{{{x}}_0}} \right)}^2} + {{\left( {{{y - }}{{{y}}_0}} \right)}^2}} - {R_0} 。$

(11)

由于Level-set函数 $ \varphi \left( {x,t} \right) $ 为距离函数，管道边界的函数值 $ \varphi \left( {x,t} \right){\text{ = 0}} $ ，管道内部任意点的距离函数 $ \varphi \left( {x,t} \right) $ 为负值，流体域内的距离函数值 $ \varphi \left( {x,t} \right) $ 为正值。在管道边界附近，可通过判断边界附近网格点距离函数 $ \varphi \left( {x,t} \right) $ 的正负值来捕捉边界位置。例如，当网格点 $ \left( {i,j} \right) $ 的距离函数为负值（流体域内的点），且任意相邻网格点的距离函数为正值（固体域内的点），则判断点 $ \left( {i,j} \right) $ 为最靠近管道边界的点，也称为浸入边界法中的加力点。

在本文采用IB方法中，海洋波流对管道的作用力由体积力代替，体积力的大小由流场和边界条件决定。添加了体积力 $ f $ 后流体的动量方程为：

$ \frac{{\partial u}}{{\partial t}} + \nabla \cdot \left( {uu} \right) = - \nabla p + v{\nabla ^2}u + f ，$

(12)

动量方程的时间离散形式为：

$ \frac{{{u^{{{n + 1}}}} - {u^{{n}}}}}{{\Delta t}} = \nabla \cdot \left( {uu} \right) - \nabla p + v{\nabla ^2}u + {f^{{{n + 1}}}} 　，　　 $

(13)

在管道边界上，计算体积力 $ {f_{n + 1}} $ 需要满足边界速度条件，即 ${u^{{{n + 1}}}} = {V^{{{n + 1}}}}$ ，即

$ {f^{{{n + 1}}}} = \frac{{{V^{{{n + 1}}}} - {u^{{n}}}}}{{\Delta t}} - \left( {\nabla \cdot \left( {uu} \right) - \nabla p + v{\nabla ^2}u} \right)。$

(14)

式(14)通过在加力点上施加速度边界条件，计算管道边界体积力来实现加力过程。在计算体积力过程中，加力点的位置通常不与边界位置重合，所以可根据周围固体点速度和流体点速度共同插值得到加力点速度。以图3为例，首先通过Level-set函数 $ \varphi \left( {x,t} \right) $ 得到管道界面附近的加力点0，利用Level-set函数 $ \varphi \left( {x,t} \right) $ 的外法向矢量 $ \vec n = {\left. {\left( {{{\nabla \varphi } \mathord{\left/ {\vphantom {{\nabla \varphi } {\left| {\nabla \varphi } \right|}}} \right. } {\left| {\nabla \varphi } \right|}}} \right)} \right|_{\varphi = 0}} $ 寻找边界固体点1，点1与点0之间的连线垂直于管道边界。设加力点0的坐标为 $ \left( {{x_0},{y_0}} \right) $ ，则固体点1的坐标为：

$ \left( {{x_1},{y_1}} \right) = \left( {{x_0} - {n_{\text{x}}}{\varphi _0},{y_0} - {n_{\text{y}}}{\varphi _0}} \right)。$

(15)

式中： $ {n_{\text{x}}} $ 和 $ {n_{\text{y}}} $ 分别为法向矢量沿 $ x $ 方向和 $ y $ 方向的分量； $ {\varphi _0} $ 为加力点的Level-set函数值。加力点相邻流体点（流体点2和3）的位置 $ \left( {{x_2},{y_2}} \right) $ 和 $ \left( {{x_3},{y_3}} \right) $ 公式如下：

$ \left\{ {\begin{array}{*{20}{c}} {\left[ {{x_2},{y_2}} \right] = \left[ {{x_0},{y_0} + {\text{sign}}\left( {{n_{\text{y}}}} \right) \cdot \Delta y} \right]}，\\ {\left[ {{x_3},{y_3}} \right] = \left[ {{x_0} + {\text{sign}}\left( {{n_{\text{x}}}} \right) \cdot \Delta x,{y_0}} \right]} 。\end{array}} \right. $

(16)

在管道壁面附近时，速度场近似满足线性分布，所以近壁流场采用多项式插值方法来确定加力点速度。设加力点0速度 $ {u_{{\text{0i}}}}\left( {i = 1,2,3} \right) $ ，由周围固体点和流体点速度线性插值得到：

$ {u_{{\text{0i}}}} = {b_{\text{1}}} + {b_{\text{2}}}{x_{\text{0}}} + {b_{\text{3}}}{y_{\text{0}}} 。$

(17)

式中： $ {b_1} $ $ {\mathrm{b}}_{1} $ ， $ {b_2} $ $ {\mathrm{b}}_{2} $ 和 $ {b_3} $ $ {\mathrm{b}}_{3} $ 为插值系数，表达式为：

$ \left[ {\begin{array}{*{20}{c}} {{b_{\text{1}}}} \\ {{b_{\text{2}}}} \\ {{b_{\text{3}}}} \end{array}} \right] = {{\mathbf{A}}^{{{ - 1}}}}\left[ {\begin{array}{*{20}{c}} {{u_{{\text{1i}}}}} \\ {{u_{{\text{2i}}}}} \\ {{u_{{\text{3i}}}}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\text{1}}&{{x_{\text{1}}}}&{{y_{\text{1}}}} \\ {\text{1}}&{{x_{\text{2}}}}&{{y_{\text{2}}}} \\ {\text{1}}&{{x_{\text{3}}}}&{{y_{\text{3}}}} \end{array}} \right]^{{{ - 1}}}}\left[ {\begin{array}{*{20}{c}} {{u_{{\text{1i}}}}} \\ {{u_{{\text{2i}}}}} \\ {{u_{{\text{3i}}}}} \end{array}} \right]\left( {i = 1,2,3} \right)。$

(18)

通过求解系数矩阵 $ {\mathbf{A}} $ ，得到加力点0上的速度信息 $ \left({u}_{0\text{1}}，{u}_{0\text{2}}，{u}_{0\text{3}}\right) $ ，即加力点0的速度 $ \left( {u,v,w} \right) $ ，然后根据式(14)求解对应的体积力。此外，由第n步时间推进到第n+1步时，还可以实时监测管道和流体相对运动，实现重复迭代，增强数值计算稳定性。

2 海洋管道绕流研究

在海洋管道静止绕流研究中，定义雷诺数Re为 $ {{{U_\infty }D} \mathord{\left/ {\vphantom {{{U_\infty }D} \nu }} \right. } \nu } $ ， $ {U_\infty } $ 为海水运动速度， $ D $ 为管道直径。设置计算域为100D×70D，管道周围网格尺寸为0.02D，总网格数为512×384，无量纲时间步长为0.001，边界条件设置进出口边界条件。当雷诺数分别为40和100时，计算结果如图4所示。当Re=40，静止管道后缘会形成2个不脱落的对称涡；当Re=100时，在管道尾部附着的尾涡交替脱落，并对管道的作用力出现周期性振荡，尾涡为卡门涡街结构。

表1列出了Re＝40时管道平均阻力系数 $ {C_d} $ 及尾迹涡长度。当Re＝100时，管道平均阻力系数 $ {C_d} $ 、升力系数的脉动幅值 $ {C'_{\text{L}}} $ 以及Strouhal数如表2所示。其中， $ {\text{St}} = {{fD} \mathord{\left/ {\vphantom {{fD} {{U_\infty }}}} \right. } {{U_\infty }}} $ ，f 为尾涡脱落频率。由表1和表2可知，计算结果与其他研究者的数值结果非常接近，说明本文算法能准确模拟静止海洋管道与流体之间的相互作用过程。

3 海洋管道涡激振动研究

作为典型的流固耦合问题，海洋管道涡激振动的实验和数值模拟研究较多^[20]。本文利用CLSVOF-IB方法研究模拟海洋管道的涡激振动特性，如图5所示。海洋管道简化为一个弹性固定圆柱绕流模型，圆柱直径为D。该模型只在横向振动，侧向位移为Y₀，无量纲化的控制方程为：

$ \ddot y + 2\zeta \left( {\frac{2}{{{U_{red}}}}} \right)\dot y + {\left( {\frac{{2{\text{π}} }}{{{U_{red}}}}} \right)^2}y = \frac{1}{{2n}}{C_L}。$

(19)

式中：y为无量纲化位移Y₀/D； $ \zeta $ 为阻尼比， $ \zeta = {0.5c} / {\sqrt {km} } $ ，c为阻尼系数，k为刚度；C_L为升力系数， $ {C_{\text{L}}} = {{{F_y}}/ {\left( {0.5\rho DL{U^2}} \right)}} $ ，F_y为流体作用力；n为质量比， $ n = {m/ {\left( {\rho {D^2}} \right)}} $ ，m为圆柱质量；U_red为约化速度， $ {U_{{\text{red}}}} = {{{U_\infty }}/ {\left( {{f_n}D} \right)}} $ ， $ {U_\infty } $ 为流体速度； $ {f_n} $ 为共振频率， $ {f_{\text{n}}} = {{\sqrt {{k / m}} }/ {\left( {2\pi } \right)}} $ 。参照文献[21-23]算例，选择质量比n为117.1，阻尼比 $ \zeta $ 为0.0012。

在海洋管道涡激振动算例中，计算域、网格尺寸分布、时间步长和进出口边界条件设置均与管道绕流相同 $ 0.005\mathrm{D}/{\mathrm{U}}_{\mathrm{\infty }} $ 。根据文献[23]，选择5组算例参数，如表3所示。通过选择不同雷诺数和约化速度U_red，得到圆柱不同的涡脱落频率。当涡脱落频率与系统共振频率接近时，圆柱振动幅值明显增加，对应频率也称为“锁频频率”；否则当2个频率相差较多时，圆柱振动幅值较小。

在稳定状态下，涡激振动圆柱的最大振幅和振动频率如图6所示。对于算例A、算例C和算例E，涡激振动圆柱振动位移随时间变化情况如图7所示。结果表明，当系统固有频率在锁频频率之内，圆柱侧向位移明显增大，无量纲化位移达到0.4；当系统固有频率在锁频频率之外，圆柱侧向位移较小，最大位移约为0.003。总体上，当系统固有频率远离尾涡脱落频率，圆柱振幅较小，阻力系数振荡较小；但当系统处于“锁频频率”状态（算例C和算例D），圆柱振动幅值明显的增加，该结果与文献[23]结论相同。

4 海洋管道对波面影响分析

为验证CLSVOF-IB数值方法对海洋波面的捕捉程度，设置管道与海洋波面相互作用的算例。以入口速度U和管道直径D为约化参数，设置雷诺数Re为2.7×10⁴，Fr参数分别为0.8，1.2和1.6，其中Froude数定义为 $ \sqrt {{{{U^2}}/{\left( {gD} \right)}}} $ ，然后分析不同Fr数时管道对流体界面的影响。如图8所示，设置计算域为60D×30D×10D，对应方向上的网格数为512×384×256，同样使用局部网格加密技术，在管道和流体/空气波面（y/D=0）附近使用加密网格，网格分辨率为0.02D×0.02D×0.02D。

当Fr数为1.2时，管道绕流对应自由波面变化的俯视图如图9所示，当Fr数分别为0.8，1.2和1.6时，管道对海洋波面影响的俯视图如图10所示。结果表明CLSVOF-IB数值算法能够较好捕捉海洋管道与流体自由界面的相互作用，计算多相流界面时有较高精度。

5 结　语

针对海洋管道在外界激振力引发的流致振动问题，本文结合CLSVOF方法和IB方法优势，建立模拟海洋管道与流体相互作用的流固耦合模型，采用Level-set距离函数的特性确定海洋管道边界，利用浸入边界法计算流体介质对管道的作用力，采用CLSVOF方法模拟海洋的自由波面。通过改变流体Fr数，分析了管道与海洋自由波面的相互作用过程。结果表明CLSVOF-IB方法能够准确有效分析海洋管道的涡激振动特性，并且能较好捕捉海洋管道与自由波面的相互作用过程。