0 引　言

为了应对全球气候变化，响应各国对能源开发和环境保护的多重需求，开发利用环保可再生能源成为瞩目热点，其中海洋资源以其“资源丰富、可持续性开发利用”等优点受到越来越多的关注^[1]。

深水油气资源的勘探与开采离不开海洋结构物的支持与保护^[2]，流体与结构相互作用产生交替的漩涡脱落，会在结构表面产生周期性的脉动力导致其发生自激振动。涡激振动和驰振是最常见的自激振动现象，当涡脱频率锁定在结构固有频率时会产生涡激振动^[3]，涡激振动在“共振区间”振动幅值较大，但具有一定区间限制性；而驰振不存在“共振区间”，随着流速的增大，振幅会持续增加，较适合于能量收集^[4]。因此国内外学者对流致振动开展大量研究，在参数敏感性分析中发现雷诺数Re、约化速度U_r、结构形状等是影响流致振动类型的重要因素^{[5 − 6]}，其中研究最多的是发生涡激振动的圆形柱体^{[7 − 8]}，而具有棱角的结构形状可能会发生不同振动模式。方形截面结构常见于桥梁建筑、海上平台、系泊系统等许多海洋工程应用中，而流致振动是诱发海洋工程结构物疲劳损伤等工程问题的重要因素，但随着对流致振动现象认识的加深^[9]，学者们认识到流致振动所获取的能量是完全清洁的可再生能源，适用范围较广，具有显著优势和利用前景^[10]，因此对结构开展流致振动研究具有重要工程意义。

早期的驰振研究集中在空气动力学领域，Hartog^[11]首次提出了覆冰输电导线结构驰振的起振判断依据，随后Parkinson和Smith^[12]建立了准定常理论预测方柱的驰振响应，White^[13]较早的利用有限元方法对结构驰振进行研究。近年来，关于方柱的驰振响应成为研究热点，邓见等^[14]在低雷诺数下采用任意拉格朗日-欧拉（ALE）法对方柱流致振动进行数值研究，Su^[15]对低雷诺数Re=100下的方柱自由振动进行了数值模拟，在弹性支撑情况下捕捉到多种振动模式。Barrero^[16]在低质量比方柱的流致振动实验中，观察到了涡激振动与驰振共存的现象，Manzoor^[17]发现质量比会影响方柱的涡激振动向驰振的转变方式，Li等^[18]通过线性稳定性分析和直接数值模拟方法得出，在低雷诺数下方柱横向驰振的本质为一种单自由度振动，练继建^[19]、胡德江^[20]等对方柱流致振动进行试验、数值模拟研究，以探究影响振动的因素。

综上所述，非圆截面柱体流致振动比圆柱涡激振动复杂得多，对于非圆截面柱体的流致振动研究有待进一步深入。在此背景下，本文利用用户自定义函数（UDF）对Fluent软件进行二次开发，通过数值模拟研究不同来流角度下的方柱流致振动，通过分析振动响应以及尾流结构等特性，探究截面几何特征、来流角度、流速等参数对方柱流致振动特性的影响规律，以期为相关海洋领域设计提供参考。

1 数值方法 1.1 控制方程

为了解释流体动力学中的流致振动，采用无量纲的N-S方程和粘性不可压缩流动的连续性方程：

$ \frac{{\partial {u_i}}}{{\partial t}} + \left( {{u_j} - {\omega _i}} \right)\frac{{\partial {u_i}}}{{\partial {x_j}}} - \frac{1}{{{Re} }}\frac{{{\partial ^2}{u_i}}}{{\partial {x_j}\partial {x_j}}} + \frac{{\partial p}}{{\partial {x_i}}} = 0，$

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0。$

(1)

式中：u_i为流体速度；ω_j为第j个方向网格速度张量；p为压力；x_i为笛卡尔坐标分量；雷诺数Re=DU_∞/v，v为流体运动粘度；U_∞为自由流速度。

本文采用剪切压力传输SST k-ω湍流模型，考虑湍流剪应力的输运效应，能对由压力梯度引起的分离现象模拟更加精准，表达式为：

$ \begin{gathered} \frac{{\partial (\rho k)}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}({\Gamma _k}\frac{{\partial k}}{{\partial {x_j}}}) + {G_k} - {Y_k} + {S_k}，\\ \frac{{\partial (\rho \omega )}}{{\partial t}} + \frac{{\partial (\rho \omega {u_i})}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}({\Gamma _\omega }\frac{{\partial \omega }}{{\partial {x_j}}}) + {G_\omega } - {Y_\omega } + {D_\omega } + {S_\omega } 。\end{gathered} $

(2)

式中：G_k为层流速度梯度产生的湍流动能；G_ω均为比耗散率；Γ_k、Γ_ω均为方程扩散率；Y_k、Y_ω均为扩散产生的湍流项；D_ω、S_k、S_ω均为正交发散项。

斯特劳哈尔数S_t为一个表征漩涡脱落的无量纲数，其与涡泄频率、结构直径和雷诺数有关，表达式为：

$ {S_t} = \frac{{{f_s}D}}{U}。$

(3)

式中：f_s为漩涡脱落频率，Hz；D为柱体特征长度，m；U为流速，m/s。

阻力系数${\overline C _{{d}}}$和升力系数$ C_{_{\text{l}}}^{'} $为表征顺流向和横向流体力大小的参数。

$ C_{_{\text{l}}}^{'} = \frac{{2{F_{{L}}}}}{{\rho {U^2}D}}，$

$ {\overline C _{{d}}} = \frac{{2{F_{{D}}}}}{{\rho {U^2}D}}。$

(4)

式中：F_L、F_D分别为结构单位长度上的升力、阻力，N；ρ为流体密度，kg/m³。

1.2 控制方程

方柱的横流向动力响应由二阶常微分方程无量纲形式控制：

$ {\ddot y_{\left( t \right)}} + 4\pi {F_n}\zeta {\dot y_{\left( t \right)}} + {\left( {2{\text{π}} {F_n}} \right)^2}{y_{\left( t \right)}} = {{C_l^{'}} \mathord{\left/ {\vphantom {{C_l^{'}} {2{m^*}}}} \right. } {2{m^*}}} 。$

(5)

式中：$ {\ddot y_{\left( t \right)}} $、$ {\dot y_{\left( t \right)}} $、$ {y_{\left( t \right)}} $为结构的瞬时加速度、速度、位移；F_n为结构的归一化固有频率，F_n=f_nD/U；m^*为结构的质量比，m^*=m/ρD²，其中m为方柱单位长度的质量；P为当前表面压力，Pa；P_∞为自由表面压力，Pa。

采用四阶Runge-Kutta法对控制方程进行求解，计算出每个时间步下结构对应的运动状态，基于此控制柱体运动和网格的更新。

$ \left\{\begin{split} & {\ddot y_{\left( t \right)}} = \frac{{{F_{{\text{L}}\left( t \right)}}}}{m} - 2\zeta {\omega _n}{\dot y_{\left( t \right)}} - \omega _n^2{y_{\left( t \right)}}，\\ & {\dot y_{\left( {t + \Delta t} \right)}} = {\dot y_{\left( t \right)}} + \frac{{\Delta t}}{6}\left( {{\beta _{\text{1}}} + 2{\beta _2} + 2{\beta _3} + {\beta _{\text{4}}}} \right)，\\ & {y_{\left( {t + \Delta t} \right)}} = {y_{\left( t \right)}} + {\dot y_{\left( t \right)}}\Delta t + \frac{{{{\left( {\Delta t} \right)}^2}}}{6}\left( {{\beta _1} + {\beta _2} + {\beta _3}} \right)。\end{split} \right.$

(6)

式中，Δt为计算时间步长，中间函数的表达式为:

$ {\left\{\begin{split} & {\beta _1} = \frac{{{F_{y\left( t \right)}}}}{m} - 2\zeta {\omega _n}{\dot y_{\left( t \right)}} - \omega _n^2{y_{\left( t \right)}} ，\\ & {\beta _2}=\frac{{{F_{y\left( t \right)}}}}{m} - 2\zeta {\omega _n}\left( {{{\dot y}_{\left( t \right)}} + \frac{{\Delta t}}{2}{K_1}} \right) - \omega _n^2\left( {{y_{\left( t \right)}} + \frac{{\Delta t}}{2}{{\dot y}_{\left( t \right)}}} \right)，\\ & {\beta _3} =\frac{{{F_{y\left( t \right)}}}}{m} -2\zeta {\omega _n}\left( {{{\dot y}_{\left( t \right)}}+ \frac{{\Delta t}}{2}{K_2}}\right) - \omega _n^2\left( {{y_{\left( t \right)}} +\frac{{\Delta t}}{2}{{\dot y}_{\left( t \right)}} +\frac{{{{\left( {\Delta t}\right)}^2}}}{4}{K_1}} \right)，\\ & {\beta _4} =\frac{{{F_{y\left( t \right)}}}}{m} -2\zeta {\omega _n}\left( {{{\dot y}_{\left( t \right)}} + \Delta t{K_3}} \right) -\omega _n^2\left( {{y_{\left( t \right)}}+\Delta t{{\dot y}_{\left( t \right)}} + \frac{{{{\left( {\Delta t} \right)}^2}}}{2}{K_2}} \right)。\end{split}\right.} $

(7)

2 数值模型 2.1 模型参数

为避免结构在振动过程中出现漩涡碰壁现象从而影响计算的准确性，本文的整体计算域大小设置为50D×40D。其中，D为结构的特征长度，计算域入口距离振动结构20D，出口距离振动结构30D，上下边界对称距离结构20D，本文流体域的阻塞比为2.5%，可忽略计算流体域选取对计算结果的影响。为了避免负网格造成的计算发散，选择了嵌套网格和动网格相结合的研究方法。计算域入口采用均匀流，侧向边界施加对称性条件，出口设置为完全出流条件；求解器采用基于压力的压力耦合求解器，算法为Couple算法；考虑计算结果的精度，采用二阶迎风格式进行空间离散、一阶隐式格式进行时间离散，柱体运动方程的求解用四阶Runge-Kutta并通过UDF实现。

考虑参数y⁺≈1的影响，在柱面附近的第一级网格间距设置为离柱面0.006D，在法线方向上，网格间距以1.05的扩展速率增加，本文计算域网格质量均可保证在0.9以上，图1为本文验证网格独立性所采用的4套网格。本文整体网格分为2种类型，即可以定义为背景网格（灰色网格）和黑色网格（前置网格），嵌套网格可避免传统动网格技术面临的网格动态更新时易出现的负体积问题，始终保持较高网格质量。

以Re=2.2×10⁴下的方柱绕流为例，对4种疏密程度不同的网格进行网格收敛性分析。表1给出4套网格密度计算得到的平均阻力系数、均方根升力系数和无量纲涡脱频率斯特劳哈尔数及各结果的相对误差情况。由表可知，当网格数量达到13万以上后计算误差较小，结果趋于稳定。综合计算成本和效率，该计算工况选取178115网格密度进行模拟。

2.2 模型验证

为验证本文选取计算方法、网格划分方式的正确性，选取方柱在Re = 2.2×10⁴工况的数值模拟结果与相关参考文献的研究结果进行对比，观察表2可知，数值结果与各参考文献中的参考数值较为接近，证明本文后续开展的数值模拟具有可靠性。

图2给出了Re = 2.2×10⁴时，方柱固定绕流下的尾流结构图，在方柱中观察到P+S漩涡脱落模式，即一个振动周期内有一对漩涡和一个单涡脱落。

方形截面为目前研究和实现驰振响应的理想模型，为了探究方柱的流致振动特性，本文主要对5种不同来流角度（0°、15°、22.5°、30°、45°）下的方柱流致振动进行数值模拟，如图3所示。

3 计算结果 3.1 频率与振幅响应分析

图4为不同来流角度θ下方柱的振动振幅、频率比随约化速度的变化趋势，明显观察到不同来流角度下的方柱振动响应存在区别，具体表现为：

1）θ=0°时，方柱的振动响应规律与其他工况下的存在明显区别。方柱振动响应分为：涡激振动初始分支（2.08 ≤U_r ≤5.21）、涡激振动上分支（5.21<U_r ≤12.5）和驰振分支（U_r >12.5）这3个阶段，振动过程中振幅随约化速度的增加而呈现增大的趋势，频率比随约化速度的增加而呈现出先增加后减小的趋势，但频率比数值始终小于1。该工况下，方柱振动并未出现圆柱涡激振动的“锁定区间”以及“自限性”振动特性，而更多表现出低频高幅的驰振基本特性。

2）θ≠0°时（15°、22.5°、30°），方柱在这3种来流角度下，振动响应特性呈现出相似的变化规律，仅在数值上存在差异。振动分支为涡激振动初始分支（2.08 ≤ U_r ≤ 5.21）、涡激振动上分支（5.21 < U_r ≤ 11.46）以及高频高幅分支（U_r >11.46）。其中，频率比随约化速度的增加而呈现持续升高的趋势，Names等^[23]认为高频高振幅响应分支是由涡激振动与驰振共同作用形成的共振分支，上述3种工况下方柱振动受涡激振动与驰振相互影响，振动较为复杂。

3）θ = 45°时，方柱的振动响应规律与圆柱振动响应类似，方柱振动响应主要分为：涡激振动初始分支（2.08 ≤U_r ≤5.21）、涡激振动上分支（5.21 <U_r ≤11.46）和涡激振动下分支（U_r >11.46）这3个阶段。其中，频率比随约化速度的增加而呈现上升趋势，振幅随约化速度的增加而呈现先增大保持稳定再下降的趋势。在θ=45°下，方柱振幅出现类似圆柱涡激振动的“锁定区间”，体现出圆柱涡激振动“自限性”的振动特性。由此表明，θ=45°工况下的方柱流致振动响应为涡激振动响应。

3.2 相位分析

相位图能直观的描述结构振动与流场能量交换之间的关系。若升力与振动位移同步变化，则流体力做正功，流场促进柱体振动，反之则抑制振动。按照笛卡尔坐标系将相位图分为4个象限，定义当相位图位于一、三象限时，此时能量由流体流向结构，促进振动；而当相位图位于二、四象限时，相位图呈反相分布，抑制振动。

图5为方柱在来流角度θ =0°时，横流向位移和升力系数之间的相位关系图。在涡激振动分支时，如图5（a），相位图位于一、三象限，此阶段下柱体所受流体力与振动位移呈现出同相的状态，随着约化速度的增加，相位开始发生转变，如图5（b）；当响应转换为驰振时，图5（c）中相位图位于二、四象限的占比增大，此时由于能量交换不稳定，但未完全“反相”；进入驰振分支后，图5（d）中相位图出现“回归”一、三象限的趋势，此阶段振幅持续增大，能量从流场流向结构。

图6为方柱在来流角度θ=22.5°下的相位关系图，相位图整体向下偏移。在涡激振动初始分支，如图6（a）所示相位图主要分布在一、三现象；当进入高频高幅分支后，如图6（c）和图6（d）所示，相位图位于二、三象限，且随着流速的增大，第三象限占比较大，这也导致振动幅值没有出现大幅下降的现象。

图7给出了方柱在来流角度θ = 45°时的相位关系图，整体呈现从“同相”到“反相”的变化过程，在涡激振动上分支时，如图7（a）所示，此时流场能量流向结构；进入上分支后，此时能量交换较为平衡，见图7（b）和图7（c）；当进入下分支后，如图7（d）所示，此阶段处于“反相”过程，能量由结构流向流场，导致振幅减小。

3.3 流体力及流场分析

为探究方柱的升阻力变化趋势，图8为相同数值模拟条件下的孤立圆柱的受力特性进行对比，其阻力系数与升力系数随约化速度的变化过程总体相似，表现为先增后减的两段式进程。图9为不同来流角度下方柱的升阻力系数变化趋势，观察可发现，所有工况下的阻力系数变化趋势相似，整个振动过程中来流角度为45°时阻力系数最大，30°、22.5°、15°依次减小，是因为不同来流角度下的方柱流向投影长度存在差异，结构的压差阻力随投影长度的增加而增大，与孤立圆柱相比，方柱阻力系数没有出现大幅下降的现象。

方柱的升力系数变化趋势也相似，其中来流角度0°时方柱在U_r>8.33后，升力系数持续下降后有增加的趋势，但与孤立圆柱相比下降趋势减缓，是由于方柱进入驰振分支后，此时振动具有“自激性”，不受漩涡脱落和固有频率的影响；在来流角度45°下，7.29≤U_r≤11.46范围内，升力系数出现剧烈波动，是由于此时方柱处于“共振”区间，结构振动受到漩涡脱落的影响，此时结构升力波动性较大，振动不稳定。

本文对来流角度0°、15°下的方柱尾流结构图进行分析，根据不同的分支来观察漩涡类型，进一步分析柱体的振动响应与尾涡类型的关系。

来流角度0°下方柱发生驰振响应。如图10（a）为初始分支时的方柱尾涡结构，在一个完整的振动周期内观察到“P+S”的漩涡脱落模式。图10（b）为上分支时的方柱尾涡结构，此时表现为“2P+2S”的漩涡脱落模式，在一个周期内，方柱上脱落6个漩涡；图10（c）为发生驰振的方柱尾涡结构，此时漩涡脱落模式为“4P+2S”，一个周期内有10个漩涡脱落，此阶段尾流涡街较宽且整体出现偏移，这也导致了方柱呈低频振动。

来流角度15°下方柱在涡激振动与驰振共同作用发生振动响应。图11（a）所示为初始分支时的方柱尾涡结构，观察到“2S”的漩涡脱落模式，图11（b）所示为上分支时的方柱尾涡结构，此时表现为“P+S”的漩涡脱落模式，尾涡密度较小；图11（c）为涡激振动与驰振共同作用下的方柱尾涡结构，此时涡脱模式为“2S”，呈现出明显的卡门涡街规律性，但涡街较宽且尾涡强度较强，这也导致了方柱的大幅振动。

4 结　语

1）来流角度是影响方柱流致振动响应和振动分支模式的重要因素。当来流角度0°时，方柱振动分支为初始分支、上分支和驰振分支，发生驰振响应；当来流角度为15°、22.5°、30°时，振动分支为初始分支、上分支和高频高幅分支，方柱产生涡激振动与驰振相互作用的振动响应；当来流角度为45°时，振动分支为初始分支、上分支和下分支，方柱发生涡激振动响应。

2）结构振动位移与升力之间的相位关系影响结构的振动。当处于“同相”时，能量由流场流向柱体结构，促进振动的发生，当处于“反相”时，能量由柱体结构流向流场，抑制柱体的振动。

3）尾流结构对振动的影响显著。尾流涡街发生偏移且周期内漩涡脱落丰富导致柱体发生低频高幅的驰振，在涡激振动与驰振共同作用下，较强的尾涡强度也能维持结构的大幅振动。

4）从能量俘获角度，可通过控制非圆柱体的来流角度实现大范围流速下的能量汲取。控制来流角度使方柱发生驰振响应，从而实现能量的汲取。