细长柔性圆柱体的涡激振动(vortex-induced vibration, VIV)问题常见于土木工程和海洋工程，如受风力作用的悬索和受海流作用的海洋立管等。在涡激力的长期作用下，细长柔性结构易发生疲劳损坏，造成经济损失，因此受到研究者的关注。从现有文献来看，二维圆柱涡激振动研究较为充分，并取得了系统性的成果。但三维细长柔性立管涡激振动的研究不足，一些机理尚不明确。

Vandiver^[1]实验研究了均匀来流条件下细长柔性圆柱的涡激振动，提出了通过横向振幅估算阻力系数的公式。Brika等^[2]通过验研究了均匀流中柔性悬线涡激振动的一阶模态，结果表明柔性悬线的涡激振动具有迟滞性，响应分支的跳跃伴随着振动位移和泻涡相位差的突变，不同响应分支的临界折合速度与Williamson等^[3]的结论基本一致。Huera-Huarte等^[4-5]开展了部分浸入水中的立管物理模型实验，顺流向的响应分支除初始分支、上端分支和下端分支外，在大折合流速区间内，再次出现比较大的振动响应；泻涡模式呈现三维性，支座附近的尾涡模式均为2S，跨中的尾涡模式则比较复杂。Chaplin等^[6]得到了与文献[4-5]相似的结论。在数值模拟方面，Huang等^[7]通过雷诺平均模型(RANS)对长细比超过3 000的柔性立管涡激振动进行了数值模拟，立管的流向弯曲对两个方向的振动响应有显著影响，靠近支座处出现了2C泻涡模式。Newman等^[8]通过直接数值模拟(DNS)研究了低雷诺数下均匀来流中长细比为4π的柔性立管涡激振动，得出两种振动模式，即行波模式和驻波模式。由于维持行波模式所需的能量较小，驻波模式最终会转换为行波模式。此外，Newman等^[8]发现行波模式下的尾涡为倾斜的，驻波模式下为平行交叉泻涡。Bourguet等^[9]通过DNS方法研究了不同雷诺数下剪切流中柔性立管的涡激振动，指出对于多模态涡激振动，振动响应在某一时刻下表现为其中的某一阶模态。Lucor等^[10]对不同剪切流条件下柔性圆柱涡激振动进行了数值模拟，泻涡模式与刚性圆柱相近，并与作用在圆柱上的流体力有直接关系。

海洋工程中也会遇到来流与立管倾斜的情况，如悬链线立管触底段。一般来说，倾斜角度α定义为来流方向与立管截平面的夹角，α=0°代表来流方向正交于立管轴线。

倾斜来流条件下固定圆柱绕流问题的研究较为充分^[11-15]。大倾角时，圆柱下游形成的展向涡排倾斜于圆柱轴线。若干研究已经验证了通过正交于圆柱轴线的流速分量来预测倾斜圆柱受力特征的可行性，这种方法称为"余弦法则"或者"独立性原理"。尽管各研究给出的独立性原理有效范围不一致，但通常认为，在α＜40°时独立性原理能提供精确的预测结果。更大倾较角情况下，轴向流速分量作用明显，独立性原理的预测值出现偏差。

相比于倾斜来流下圆柱绕流的研究，倾斜来流下细长柔性立管涡激振动的研究则少之又少。Bourguet等^[16]对倾斜来流下细长柔性立管涡激振动进行了直接数值模拟，并考虑了不同顶张力的影响。大顶张力时，立管的顺流向弯曲很小，横流向振动呈现2阶模态，独立性原理适用于预测圆柱振动、流体力等；顶张力小时，立管的顺流向弯曲较大，独立性原理出现了偏差。

本文旨在对大倾斜角(α=60°)来流条件下的细长柔性立管涡激振动进行直接数值模拟，深入分析立管的振动响应、振动能量谱、相位差、振动轨迹、流体力、能量传递参数等，以及尾涡结构与振动响应之间的关系。

1 浸入边界法和参数设置 1.1 数值方法

本文应用基于浸入边界法(immersed boundary method，IBM)的三维水动力并行计算程序CgLES_IBM，并结合隐式结构动力计算程序X-code，对倾斜来流作用的细长柔性立管涡激振动开展数值模拟研究。

流体运动的控制方程如下

$ \left\{ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} = - \nabla \mathit{\boldsymbol{uu}} - \nabla p + \nu {\nabla ^2}\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{f}}_b} + {\mathit{\boldsymbol{f}}_i}\\ \nabla \mathit{\boldsymbol{u}} = 0 \end{array} \right. $

(1)

式中：u为速度；t为时间；p为压强；ν为运动粘滞系数；f_b为附加体积力矢量，代表流固耦合边界条件；f_i为重力矢量。

对以上控制方程采用具有二阶精度的Adams-Bashforth格式进行时间离散，可得控制方程的守恒形式如下

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{u}}^{n + 1}} = {\mathit{\boldsymbol{u}}^n} + {\delta _t}\left( {\frac{3}{2}{\mathit{\boldsymbol{h}}^n} - \frac{1}{2}{\mathit{\boldsymbol{h}}^{n - 1}} - \frac{3}{2}\nabla {p^n} + } \right.\\ \;\;\;\;\;\;\;\;\left. {\frac{1}{2}\nabla {p^{n - 1}}} \right) + {\mathit{\boldsymbol{f}}^{n + \frac{1}{2}}}{\delta _t}\\ \nabla \cdot {\mathit{\boldsymbol{u}}^{n + 1}} = 0 \end{array} \right. $

(2)

式中：δ_t为时间步长，h=▽·(-uu+ν(▽u+▽u^T))由对流项与扩散项组成，附加体积力表示为

$ \begin{array}{l} {\mathit{\boldsymbol{f}}^{n + \frac{1}{2}}}{\delta _t} = D\left( {{\mathit{\boldsymbol{V}}^{n + 1}} - I\left( {{\mathit{\boldsymbol{u}}^n} + {\delta _t}} \right.\left( {\frac{3}{2}{\mathit{\boldsymbol{h}}^n} - \frac{1}{2}{\mathit{\boldsymbol{h}}^{n - 1}}} \right.} \right. - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {\left. {\frac{3}{2}\nabla {p^n} + \frac{1}{2}\nabla {p^{n - 1}}} \right)} \right)} \right) \end{array} $

(3)

式中：Vⁿ⁺¹为物面边界的速度，I(ϕ, X_i)和D(Φ, x)为插值函数，表示为

$ \mathit{\Phi }\left( {{\mathit{\boldsymbol{X}}_i}} \right) = I\left( {\phi, {\mathit{\boldsymbol{X}}_i}} \right) = \sum\limits_{x \in {g_h}} {\phi \left( x \right)} {\delta _h}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{X}}_i}} \right) $

(4)

$ \phi \left( \mathit{\boldsymbol{x}} \right) = D\left( {\mathit{\Phi }, \mathit{\boldsymbol{x}}} \right) = \sum\limits_{i = 1}^{{N_{{\rm{IBP}}}}} {\mathit{\Phi }\left( {{\mathit{\boldsymbol{X}}_i}} \right)} {\delta _h}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{X}}_i}} \right)\Delta {V_i} $

(5)

式中：ϕ(x)表示定义在流体正交网格结点x上的物理量，Φ(X_i)表示定义在可移动物面边界点X_i上的物理量。I(ϕ, X_i)将定义在正交网格上的物理量插值到物面边界点上；与之相反，D(Φ, x)将定义在物面边界点上的物理量投射到正交网格上。g_h表示正交网格点的集合，N_IBP表示插值时涉及某个正交网格点的所有的物面边界点的个数，ΔV_i为物面边界点的体积。δ_h为狄拉克函数，表示为

$ {\delta _h}\left( \mathit{\boldsymbol{x}} \right) = {d_h}\left( x \right){d_h}\left( y \right){d_h}\left( z \right) $

(6)

式中：d_h二阶空间精度的插值函数，表示为

$ {d_h}\left( r \right) = \left\{ \begin{array}{l} \frac{1}{{8h}}\left( {3 - 2\frac{{\left| r \right|}}{h} + \sqrt {1 + 4\frac{{\left| r \right|}}{h} - 4{{\left( {\frac{{\left| r \right|}}{h}} \right)}^2}} } \right), \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| r \right| \le h\\ \frac{1}{{8h}}\left( {5 - 2\frac{{\left| r \right|}}{h} - \sqrt { - 7 + 12\frac{{\left| r \right|}}{h} - 4{{\left( {\frac{{\left| r \right|}}{h}} \right)}^2}} } \right), \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;h \le \left| r \right| \le 2h\\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2h \le \left| r \right| \end{array} \right. $

(7)

针对传统浸入边界法施加边界条件精度不高的情况，Ji等^[17]提出了基于嵌入式迭代的浸入边界法，将浸入边界法嵌入到压强泊松方程的迭代求解中，利用压强的中间解比初始值更接近真实值的特点，迭代修正附加体积力，在不显著增加额外计算耗时的前提下，提高整个算法的求解精度。

实际工程中，海洋立管的长细比较大，对结构固有频率的起主要作用的是立管所受的顶张力，结构自身刚度对固有频率影响很小，可忽略。本文将细长柔性立管简化为受顶张力作用且两端铰接的弹性连续索，其动力平衡方程表示为

$ \mathit{\boldsymbol{M\ddot A}} - \mathit{\boldsymbol{TA''}} + \mathit{\boldsymbol{D\dot A}} = \mathit{\boldsymbol{F}} $

(8)

式中：M为集中质量矩阵，T为顶张力矩阵，D为阻尼矩阵，F为荷载向量，A为位移向量, 上标·和'分别表示对时间和空间的导数。

采用三维二结点索单元对立管进行空间离散，每个结点具有3个自由度。结构动力方程采用Newmark-β法进行求解，该方法为无条件稳定，具有二阶时间精度，时间步长的取值可以灵活地与求解流场的时间步长协调一致。

1.2 验证算例

本文对刚性圆柱绕流进行了三维数值模拟，圆柱长细比为L/D=10，雷诺数Re=1 000。计算域在顺流向、横流向和展向的大小为50D×50D×10D，采用768×768×128的网格对计算域进行离散，顺流向和横流向最小网格边长均为D/64，展向采用均匀网格，边长为D/12.8。

表 1给出了不同结果的圆柱背流侧基压-C_pb、平均阻力系数$ {\bar C_D} $、脉动升力系数C'_L、泻涡频率St的对比。可以发现，本文结果与文献[18, 19]的结果吻合较好，验证了本文方法的准确性。

1.3 计算参数

计算域为长方体, 顺流向、横流向和展向的大小为70D×40D×50D。入流边界为Dirichlet型边界(u=U, v=0, w=0), 出流边界为Neumann型边界(∂u/∂x=0, ∂v/∂x=0, ∂w/∂x=0), 横流向侧边界为可滑移边界(∂u/∂y=0, v=0, ∂w/∂y=0), 立管壁面为无滑移边界(u=u_s, v=v_s, w=w_s, 其中u_s、v_s、w_s表示结构物的速度), 展向为周期性边界，如图 1所示。立管两端铰接, 直径为D, 长度为50D, 质量比m^*=6。初始状态下, 立管中心距入流边界20D，距出流边界50D，距侧边界均为20D。本文Re=ρ_fUD/μ=500。其中, U为未受扰动的来流速度, 方向与立管横截面的夹角为α=60°，ρ_f为流体的密度，μ为流体的动力粘滞系数。无量纲顶张力为T=τ/ρ_fD²U²=124，τ为顶张力。为了获得更大的振动幅值，结构阻尼设为D=0。立管布置见图 2。

在x-y平面内, 采用正交网格对计算域进行空间离散，网格数量为576×256。为提高精度, 在立管附近16D×6D的范围内采用均匀正交网格进行加密。加密区内网格尺寸为D/32, 加密区外为渐变网格, 网格尺寸以指数形式增长。展向采用均匀网格, 共划分192层网格。对于立管的有限元模型, 其划分单元的数目与展向的流体网格相协调, 为192个单元。图 3所示为网格布置，右下角为立管周围的网格细节。

为了满足CFL=UΔt/Δx＜0.5的稳定性条件, 取无量纲时间步长ΔtU/D=0.001。

2 算例结果分析 2.1 结构振动响应

图 4所示为立管顺流向和横流向振动位移的均方根值。顺流向位移扣除了平均位移，其均方根振幅约为0.125D，横流向位移的均方根振幅约为0.55D，顺流向振动呈4阶振型，横流向振动呈2阶振型。这与相同模拟条件下的Bourguet等^[16]得到的结果较为吻合，如图 4所示，进一步验证本文方法的正确性。图 5给出了顺流向和横流向的位移的时空分布，从图中可以看出立管振动表现出明显的驻波振动特性，且振动的波腹数与立管模型的阶数对应。

对立管位移的时空分布进行傅里叶分析，得到振动能量谱沿展向的分布。如图 6所示，展向上的振动能量分布与振型对应，顺流向振动的能量谱体现出4阶振型的特征，横流振动的能量谱体现出2阶振型的特征。此外，频谱图呈现单谱模式，即大部分振动能量均集中在一个较窄频段上。如：顺流向主控频率集中在f_x=0.352，横流向主控频率集中在f_y=0.176，两者之比等于2。这与Bourguet等^[16]的结果f_x=0.348和f_y=0.174吻合良好。

2.2 位移相位分析

Dahl等^[20]实验研究了弹性支撑的刚性圆柱顺流向和横流向振动的同步性问题，结果发现顺流向和横流向振动的相位差与尾流模式和流体力等关系密切。Vandiver等^[21-23]进一步实验研究了剪切来流中细长柔性立管的两向振动的同步性问题，指出两向振动的相位差与立管展向上激励区和阻尼区的分布有一定的联系。

本文采用Huera-Huarte等^[4-5]提出的方法计算两向振动的相位差，即

$ {\varphi _{xy}} = \left( {p{\varphi _x} - q{\varphi _y}, \bmod \;360^\circ } \right) $

(8)

式中：φ_x和φ_y分别为顺流向和横流向位移的瞬时相位，可通过希尔伯特变换获得。p和q为整数，通过顺流向和横流向的瞬时振动频率确定, 即f_x/f_y=q/p。本文结果中顺流向和横流向振动频率之比约为2，故取(p, q)=(1, 2)。之前研究通常将相位差φ_xy和振动轨迹的方向联系在一起。若展向上某点的振动轨迹在横流向最大正位移时刻向上游移动，则该点的振动轨迹为逆时针方向，相位差φ_xy在0°~180°；若振动轨迹在横流向最大正位移时刻向下游移动，则该点的振动轨迹为顺时针方向，相位差φ_xy在180°~360°。需要指出的是，只有当两向振动进入同步振动的状态后, 即φ_xy基本不随时间发生变化，两向振动相位差φ_xy才有意义。

图 7(a)为相位差φ_xy的概率密度的展向分布。可以看出，相位差的权重分布体现出明显的驻波特征。立管展向位置z/D∈(0, 12.5)和z/D∈(25, 37.5)范围内的相位差均在φ_xy∈(0°, 180°), 振动轨迹为逆时针方向；立管展向位置z/D∈(12.5, 25)和z/D∈(37.5, 50)的相位差均在φ_xy∈(180°, 360°)，振动轨迹为顺时针方向。而在横流向振动腹点附近，即z/D=12.5和z/D=37.5处，相位差φ_xy发生180°跳跃。

图 7(b)是立管展向不同位置处的振动轨迹图，从下到上z/D依次为5、10、12.5、15、20、25、30、35、37.5、40、45。可以看出，由于跨中处两向振动均为结点位置，此处立管振幅很小，轨迹形状难以定义。以跨中为界，跨中以下部分和以上部分均是一个完整的轨迹变换过程。进一步观察发现，z/D=12.5和z/D=37.5处，φ_xy发生180°的切换，z/D∈(0, 12.5)和z/D∈(25, 37.5)为逆时针的"8"字形，而z/D∈(12.5, 25)和z/D∈(37.5, 50)则为顺时针的"8"字形。其中，z/D=12.5和z/D=37.5处的"8"字形很狭长，这是由于此处是横流向振动的波腹和顺流向振动的波结造成的。

2.3 流体力与流固能量传递 2.3.1 流体力

流体力是表征流固耦合系统的非线性平衡态的重要参数，可以用来认识流固能量传递的时空特性。图 8给出了流体力的时空分布，可以看出升、阻力均表现为驻波和行波的混合模式，且行波模式在阻力中占优，而驻波模式在升力中占优。流体力这种分布特点与立管振动模式和尾涡有着紧密的联系，图 9给出了流体力系数沿立管展向的分布。可以看出，阻力系数平均值在1.2~2.0变化，阻力系数均方根值在0.4~1.3变化，升力系数均方根值在0.25~1.1。横流向振动波腹附近，阻力系数平均值和均方根值较大，而升力系数的均方根值较小。然而，仔细观察发现，流体力分布并不关于横流向波腹对称，这与相位差φ_xy和振动轨迹在波腹处发生切换相关。

2.3.2 流固能量传输

流固之间的能量交换不仅对VIV现象的机理解释有理论价值，而且对工程实际有指导意义。

流体与固体能量传递时均值的计算公式为

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{C_{{\rm{fv}}}}\left( z \right) = \\ \frac{{\frac{2}{T}\int_{{T_s}} {\left( {{{C'}_D}\left( {z, t} \right){{\dot A}_x}\left( {z, t} \right) + {C_L}\left( {z, t} \right){{\dot A}_y}\left( {z, t} \right)} \right){\rm{d}}t} }}{{\sqrt {\frac{2}{T}\int_{{T_s}} {\left( {\dot A_x^2\left( {z, t} \right) + \dot A_y^2\left( {z, t} \right)} \right){\rm{d}}t} } }} \end{array} $

(9)

式中：$ {{{\dot A}_x}} $和$ {{{\dot A}_y}} $代表立管顺流向和横流向的振动速度，T_s为统计时间，取整数倍的振动周期。C_fv为时均能量传递参数，用来量化流体力对立管做功的"功率"。C_fv为正数，意味着流体向结构振动提供能量，激发结构振动，C_fv为负数，表明流体对立管振动起阻尼作用，抑制结构振动。

图 10为C_fv沿立管展向的分布。如图(10)所示，激励区范围约为z/D∈(0, 20)和z/D∈(25, 40)，阻尼区范围约为z/D∈(20, 25)和z/D∈(40, 50)。C_fv的分布受到横流向和顺流向振动的共同影响。z/D=25和z/D=50(z/D=0)处为两向振动的结点，其附近偏下一定范围内能量传递参数为负值，立管涡激振动处于阻尼状态；z/D=12.5和z/D=37.5处为横流向振动的波腹，能量传递参数达到最大，涡激振动处于激发状态。文献[22]的正交来流细长柔性圆柱涡激振动的实验结果表明，激励区往往对应着逆时针的振动轨迹，而阻尼区对应顺时针的振动轨迹。本文中，C_fv的分布并非完全印证文献[22]的结论。在逆时针振动区域，z/D∈(12.5，25)达到最大值(正值)，且沿展向几乎不变，这与Modarres-Sadeghi等^[22]的结论一致。在顺时针振动区域z/D∈(12.5, 25)，C_fv迅速减小，但仅在z/D∈(20, 25)区域内出现了负值，这与Modarres-Sadeghi等^[22]的结论不同。本文认为，这与倾斜来流造成的升力在顺时针振动区域较小有关，如图 9(c)所示。尽管立管横流向振动的速度在顺时针和逆时针振动区域相同，但升力却在顺时针振动区域较小，造成阻尼效应减弱。需要说明的是，由于受到端部支撑条件的影响，立管端部的顺时针和逆时针振动区域与跨中相应区域的能量传递系数有稍微的区别。

2.4 尾涡模式

图 11为某一时刻不同展向截面的尾涡图。可以看出，立管的泻涡模式为2S模式，即一个振动周期内立管两侧分别脱落一个漩涡y。

Bourguet等^[16]对60°倾斜来流条件下固定圆柱绕流进行了数值模拟，认为圆柱后的尾涡为明显的倾斜泻涡模式，如图 12(a)所示。然而细长柔性立管涡激振动的尾涡模式(见图 12(b))既区别于这种模式，又区别于平行泻涡模式。由图 12(b)可知，倾斜来流下立管涡激振动的尾流场表现出倾斜泻涡特征，但与圆柱绕流的倾斜泻涡模式相比，倾斜角度明显偏小，且倾斜的方向相反。结合图 5和图 8，这种倾斜泻涡模式导致了流体力时空变化呈现一定的行波模式。同时，由于立管长细比没有达到界限值，故而立管位移响应表现为驻波模式。受此影响，流体力时空变化也体现出一定的驻波模式。

3 结论

1) 受限于立管的长细比，位移响应呈现出驻波振动特征，顺流向为4阶振型，横流向为2阶振型，顺流向主控频率约为横流向的2倍。

2) 立管涡激振动的频谱呈现为单谱模式，所有的振动能量都集中在一个窄频段上。

3) 两向振动相位差在顺流向振动结点处出现了180°的相位切换。跨中上、下两侧的相位变化呈现周期性，立管展向位置上z/D∈(0, 12.5)和z/D∈(25, 37.5)范围内的振动轨迹为逆时针的"8"字形，其他部分为顺时针的"8"字形。

4) 流体力的时空分布呈现行波和驻波的混合模式。横流向振动波腹附近，阻力系数平均值和均方根值较大，而升力系数的均方根值较小。此外，流体力分布并不关于横流向波腹对称。

5) 研究了时均能量传递系数沿立管展向的分布，发现：振动结点以下一定范围内涡激振动处于阻尼状态。在逆时针振动区域，能量传递系数达到最大值，且沿展向几乎不变；在顺时针振动区域，能量传递系数迅速减小，但仅在部分区域内出现负值，这与倾斜来流造成的升力在顺时针振动区域较小有关。

6) 倾斜来流条件下细长柔性立管涡激振动的尾流模式表现为倾斜泻涡模式，但与倾斜来流条件下圆柱绕流的倾斜泻涡模式相比，倾斜角度较小，且倾斜方向相反。这种倾斜的泻涡模式决定了流体力的时空分布呈现出一定的行波特征。但由于立管位移响应表现为驻波模式，流体力时空变化也体现出一定的驻波模式。