0 引　言

为应对全球气候变化并减少碳排放，开发和利用可再生清洁能源已成为全球关注的焦点。海洋资源具有“资源丰富、可持续开发利用”的显著优势，逐渐成为助力我国实现“双碳”目标的重要蓝色路径。当海流流经海洋结构物时会在结构物表面交替脱落漩涡从而诱发振动，这一现象称为流致振动^[1]，流致振动常发生在海洋系泊缆、海洋锚链、立管或桥梁建筑等结构^[2]，是海洋工程领域中较为常见的现象。

在流致振动现象中，流体的动能与结构振动的机械能会持续进行能量交换，海洋结构物流致振动过程蕴藏着丰富的能量，这些都是海洋能收集的潜在方向。自从卡门涡街被发现后，流致振动特性成为国内外的研究热点，涡激振动是目前流致振动研究中开展最为完善的研究方向^[3]。Khalak等^[4]开展涡激振动试验，将涡激振动的响应分支分为初始分支、上分支和下分支，学者们也关注到质量比、阻尼等参数敏感性对振动响应特性的影响。Bernitsas等^{[5 − 6]}提出一种利用流致振动的新型海流发电装置（Vortex Induced Vibration Aquatic Clean Energy，VIVACE），其原理是利用柱体结构的流致振动带动发电机进行发电，从而将海流能转化为电能。为了更深入地研究流致振动能量收集原理，研究学者们开展了大量的流致振动基础研究^[7]。

不同于圆柱的涡激振动响应，非圆截面柱体动力学响应特性较为复杂，可能会出现一种高振幅低频率的驰振响应现象，也可能会存在不同的振动模式。张乾坤^[8]等基于统计能量法研究了加强形式对翼型结构流致振动的影响。早期的驰振研究主要集中在空气动力学领域。Den等^[9]首次提出了判断覆冰输电导线结构驰振起振的依据；随后，Parkinson等^[10]建立了用于预测方柱驰振响应的准定常理论；White等^[11]较早地采用有限元方法对结构驰振进行研究。由于三棱柱结构几何存在尖角，其有利于流体的分离从而触发边界层转捩，是诱发驰振的理想结构。近年来，学者们针对非圆截面柱体开展大量的流致振动研究，Ding等^[12]在较大雷诺数范围内开展数值模拟研究，分析非圆截面柱体的振动响应转变特性。Alonso等^[13−14]进行了大量的三棱柱流致振动风洞试验，确定了三棱柱的驰振失稳关系；Kumar等^[15]利用数值方法模拟研究雷诺数在10≤Re≤250范围内正三棱柱的二维流动，田依农等^[16]研究了截面长细比对三棱柱流致振动的影响，练继建等^[17]对非圆截面柱体流致振动进行试验，对刚度、阻尼比、雷诺数等参数进行研究。

目前针对圆柱截面的涡激振动研究尤为完善，非圆截面柱体的流致振动主要集中在振动模式、振动响应转变特性以及物理参数对非圆截面流致振动的影响研究，而非圆截面柱体流体力、尾流脱落形式与振动响应模式之间的耦合机理尚不清晰，且由于流体流动摩擦及转化中存在的机械摩擦造成能量损失，因此如何提高流致振动中能量转换效率有待进一步研究。因此，本文针对不同来流角度下的三棱柱开展数值模拟研究，通过自定义Fluent软件中的UFD函数进行二次开发，系统分析了三棱柱在不同来流角度下的动力响应、流体力、相位特性及尾流模态的变化，揭示了三棱柱流致振动的影响规律，为流致振动能量收集设计提供了参考依据。

1 数值计算方法 1.1 数值理论

本文数值模拟采用基于二维数值计算方式，考虑弹性支撑结构单自由度振动，其振动控制方程为二阶常微分方程如下：

$ m\ddot{y}(t)+C_{\mathrm{total}}\dot{y}(t)+Ky(t)=F_{\mathrm{fluid}}(t)。$

(1)

式中：$ {{\ddot y}_{( t )}} $、$ {{\dot y}_{( t )}} $、$ {y_{( t )}} $为结构的瞬时加速度、速度、位移；m为振动结构的质量，kg；C_total为结构的阻尼系数，N·s/m；K为振动系统刚度，N/m；F_fluid(t)为结构受到流体力，N。

结构在升力频率和结构主导频率相近时会发生共振，此区间适合于能量收集，根据Bearman^[18]提出的简谐理论近似等效结构的运动状态，即：

$ f=f\mathrm{_{fluid}}=\frac{1}{{T}}，$

(2)

$ {y_{(t)}} = A \cdot \sin (2 {\text{π}} ft) 。$

(3)

式中：A为振子的振幅，m；f为系统的振动主导频率，Hz；f_fluid为升力频率，Hz。

一个振动周期内流体动能转换的电能计算式为：

$ {P}_{UL}=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}{F}_\mathrm{fluid}(t)\dot y(t)\mathrm{d}t}。$

(4)

将式(2)～式(4)代入式(1)中，最终得到功率的表达式为：

$ \begin{split} & A\cdot f=\frac{C_y\left[\tan\left(2\text{π}ft\right)\cos\phi+\sin\phi\right]D\cdot f_nU_r^2}{8\text{π}^2\zeta_{\mathrm{total}}\left(m^*+C_a\right)}= \\ & \quad\quad\quad\frac{U^2\rho DL\cdot C_y\sin\phi}{4\text{π}C\mathrm{_{total}}}，\\ & P_{UL}=8\text{π}^3m\zeta f_nA^2f^2。\end{split} $

(5)

采用四阶Runge-Kutta方法对式(1)进行离散求解，计算出每个时间步下结构对应的运动状态，并将计算得到的位移和速度作为下一步计算的初始条件，基于此控制柱体运动和网格的更新。求解过程如下：

$ \ddot {y_{\left( t \right)}} =\frac{{{F_{y(t)}}}}{m} -2\zeta{\omega _n} \dot {y_{\left( t \right)}}-{\omega _{n}^2} {y_{\left( t \right)}}, $

(6)

$ {\dot y_{(t + \Delta t)}} = {\dot y_{\left( t \right)}} + \frac{{\Delta t}}{6}({K_1} + 2{K_2} + 2{K_3} + {K_4})，$

(7)

$ {y_{(t + \Delta t)}} = {y_{\left( t \right)}} + {\dot y_{\left( t \right)}}\Delta t + \frac{{{{\left( {\Delta t} \right)}^2}}}{6}({K_1} + {K_2} + {K_3})。$

(8)

式中：Δt为计算时间步长。其中，K₁、K₂、K₃、K₄表达式为：

$ {K_1} = \frac{{{F_{y(t)}}}}{m} - 2\zeta {\omega _n}{\dot y_{(t)}} - \omega _n^2{y_{(t)}}，$

(9)

$ {K_2} = \frac{{{F_{y(t)}}}}{m} - 2\zeta {\omega _n}\left( {{{\dot y}_{(t)}} + \frac{{{{\Delta }}t}}{2}{K_1}} \right) - \omega _n^2\left( {{y_{(t)}} + \frac{{{{\Delta }}t}}{2}{{\dot y}_{(t)}}} \right)，$

(10)

$ {{K_3} = \displaystyle\frac{{{F_{y\left( t \right)}}}}{m} - 2\zeta {\omega _n}\left( {{{\dot y}_{\left( t \right)}} + \frac{{\Delta t}}{2}{K_2}} \right){\text{ }} - \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) }，$

(11)

$ {{K_4} = \displaystyle\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)} 。$

(12)

1.2 数值分析方法

本文采用剪切应力传输SSTk-ω湍流模型，充分考虑湍流剪应力的输运效应，以提高对压力梯度诱导分离现象的模拟精度。表达式如下：

$ \left\{\begin{aligned} &\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} ,\\ &\begin{split} \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{split}\end{aligned} \right. $

(13)

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

在振动特性研究中，常用约化速度代替为流速，其表达式为：

$ {U_r} = \frac{U}{{{f_n}D}} 。$

(14)

式中：U为来流速度，m/s；f_n为固有频率，Hz。

流体流经结构会在其两侧产生旋涡泄放现象，升阻力系数表达式分别为：

$ C_{_{\text{l}}}^{{'}} = \frac{{2{F_{{L}}}}}{{\rho {U^2}D}}\\{\overline C _{{d}}} = \frac{{2{F_{{D}}}}}{{\rho {U^2}D}} 。$

(15)

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

利用用户自定义函数对Fluent进行二次开发，采用Overset重叠网格在Fluent软件中建立二维流固耦合模型，结合SST k-ω湍流模型求解运动方程获取流场信息，通过四阶Runge-Kutta方法求解柱体振动方程，控制方程采用Coupled算法进行处理，借助UDF及DEFINE-CG-MOTION宏命令实现网格运动，从而实现三棱柱振动特性的数值模拟研究。

2 数值模型 2.1 数值模型建立

参考时忠民等^[19]对于流场计算域的研究，流体域中的宽度和尾流区长度需满足计算精度要求，本文的整体计算域大小设置为50D×40D，其中D为三棱柱的特征长度，计算域入口距离振动结构为20D，出口距离振动结构为30D，上下边界对称距离结构为20D，如图1所示。考虑到阻塞比B（B=d/H，H为计算域宽度，d为结构的投影长度），为避免计算过程中由于阻塞比较大产生的边界效应影响，本文流体域的阻塞比为2.5%，因此可以忽略计算流体域选取对计算结果的影响。

本文计算域入口设置为均匀流条件，即u=U，v=0；出口采用完全出流条件，侧向边界施加对称性条件。为确保边界层具有足够的分辨率，柱面附近的第一级网格间距设置为距柱面0.006D，在法线方向上网格间距按1.05的扩展速率逐渐增加。整体网格如图2所示，采用背景网格（浅色网格）和前置网格（深色网格）2种类型，结合嵌套网格与动网格方法，确保网格在运动过程中始终保持较高质量。本文主要对5种不同来流角度（0°、30°、45°、60°、90°）下的三棱柱流致振动进行数值模拟，如图3所示。通过分析振幅、振动频率等振动特性以及流场特性，以探究三棱柱的振动机理。

2.2 模型验证

对4种不同疏密程度的网格进行网格收敛性分析，以雷诺数Re=2.2×10⁴为条件，以针对三棱柱绕流进行研究，以验证本文计算方法与网格划分方式的合理性。表1为不同网格密度下计算所得的均方根升力系数、平均阻力系数、无量纲涡脱频率（斯特劳哈尔数St）及其相对误差根据对比结果，在兼顾计算时间与计算效率的基础上，本文选取网格密度为178115的方案进行模拟分析。将相关参考文献的研究结果与Re=200工况的数值模拟结果进行对比，结果如表2所示。本文模拟结果与参考文献中的数据具有良好的一致性，进一步证明了本文研究方法的有效性和可靠性。

3 三棱柱流致振动数值研究 3.1 频率与振幅响应分析

图4为不同来流角度下三棱柱的频率、振幅随约化速度U_r的变化趋势，对应结果存在显著的区别，具体表现在以下方面：

1）来流角度θ=0^o时，三棱柱的振动响应规律与其他来流角度响应存在明显区别。三棱柱振动响应主要分为：涡激振动初始分支（2.04≤U_r≤5.10）、涡激振动-驰振转变分支（5.10<U_r≤10.21）和驰振分支（U_r>10.21）3个阶段。其中，当约化速度增大时，振幅随之增大，频率比（f_s/f_n）呈现出先增加后减小的趋势。在驰振分支，由于振动频率不受固有频率和漩涡脱落的影响，振动频率呈现稳定下降趋势且数值始终小于1，振幅持续增长未出现下降趋势，该工况下三棱柱表现出低频高幅的驰振特性。

2）来流角度θ=30^o、90^o时，三棱柱的振动响应特性仅在频率比数值上存在差异。2种来流角度下的振动响应主要分为涡激振动初始分支（2.04≤U_r≤5.10）、过渡分支（5.10<U_r≤14.29）以及高振幅分支（U_r>14.29）3个阶段。其中，振幅比持续增大，频率比随约化速度的增加先保持稳定而后迅速增大，在过渡分支时，频率比保持在0.5附近波动，呈现“弱锁定”状态；值得注意的是，高振幅分支并不是低频高振幅的驰振分支，而是高频高振幅响应分支，此阶段存在多阶振动，振动成分较为复杂。

3）θ = 45^o时，振动响应分为涡激振动初始分支（2.04≤U_r≤6.12）、过渡分支(6.12<U_r≤12.25)以及高振幅分支（U_r>12.25）3个阶段，该工况下振幅和频率保持上升趋势。

4）θ = 60^o时，整个振动过程中观察不到明显的振动幅值，频率比随约化速度的增加而持续增大。该工况下振动过程中振幅比出现最大值而后减小，可认为此工况下的三棱柱流致振动响类似涡激振动响应，该工况下振幅较小且存在“自限性”。

3.2 相位分析

柱体流致振动中升力和振幅的相位差变化反映了流体力做功的正负，即能量由流场输入到柱体振动体系，或相反。以原点为中心，采用笛卡尔坐标系将相位图划分为4个象限。当相位点位于第一象限和第三象限时，结构振动位移与升力系数呈同相分布，此时能量由流体传递至结构；而当相位点位于第二象限和第四象限时，则呈反相分布，与第一象限和第三象限的情况相反。

由以上分析可知，来流角度θ=0°时三棱柱发生驰振响应，图5为三棱柱处于不同振动分支下的无量纲时程曲线以及相位图。当振动处于涡激振动初始分支时，见图5(a)，三棱柱呈现不规律振动并且振动幅值较小，随着约化速度的增加，三棱柱进入涡激振动-驰振转变分支，此阶段位移时程曲线始终与升力系数时程曲线保持一致；但振动存在不规则波动，如图5(b)、图5(c)所示；当进入驰振分支后，如图5(d)所示，位移与升力依然保持一致并呈现出稳定的周期性振动，在整个振动过程中，相位图始终位于一、三象限，柱体所受流体力始终与振动位移“同相”，可认为结构始终从流场中获取能量，这也导致了三棱柱发生高幅值的驰振响应。

由上述图4分析可知，三棱柱在来流角度θ=30°、90°下的响应趋势相似，因此，图6、图7分别给出来流角度30°、60°时三棱柱的无量纲时程曲线和相位图。在整个振动过程中，位移曲线和升力系数曲线均向上偏移，在涡激振动初始分支时，时程曲线呈现规则的多周期振动，相位图分布在一、三象限，如图6(a)所示；在过渡分支则表现为不规则的多周期振动，此时相位图在4个象限均有分布如图6(b)、图6(c)；当进入高振幅分支后，见图6(d)，此阶段一个位移周期内有多个升力系数峰值，相位图偏移在第一象限，位移与升力之间存在相位差，这导致该工况下振动幅值小于来流角度θ=0°下的振幅。

3.3 流致振动能量及流场分析

由式(14)可知，振幅和振动主频率是描述系统振动的重要参数，通过优化参数的选取，有利于提高系统的能量利用率。由于流体流动摩擦及转化中存在的机械摩擦损耗能量，可以认为流体对振子做功功率是流致振动能量收集的上限可利用功率，可作为衡量流致振动发电潜能的有效依据。不同来流角度流体对振子做功功率如图8所示，分析可知，在多数约化速度范围内0°工况功率较高，45°工况在高约化速度时达到功率最大值。

本文对来流角度为 0° 和 30° 时三棱柱的尾流结构进行了分析，通过观察不同分支的涡旋类型，探讨了柱体振动响应与尾涡类型之间的关系。

来流角度0°下的三棱柱发生驰振响应，根据不同的分支来观察漩涡类型。图9为初始分支时的三棱柱尾涡结构，在一个完整的振动周期内观察到“2P”漩涡脱落模式，在柱体上、下方分别脱落一对漩涡；如图9(b)所示为转变分支时的尾涡结构，此时观察到“2P+4S”漩涡脱落模式；当进入驰振分支后，尾涡脱落呈现更复杂的模式，如图9(c)所示，此时观察到“2T+2S”漩涡脱落模式，“2T”模式是指每半个周期内产生3个漩涡，2个符号相同的漩涡朝一个方向旋转，另一个符号相反的漩涡朝另一个方向旋转，本文观测到的三棱柱尾流类型与Ding等^[15]研究结果一致。在漩涡脱落过程中，三棱柱后沿尖角对于漩涡脱落有着显著影响，由于后延尖角的存在，使得漩涡能从分离剪切层中快速切断，没有形成长涡街，整个尾涡区斜向偏移，从而使三棱柱呈低频振动。

来流角度30°下，三棱柱在初始分支时的“P+2S”，如图10(a)所示，图10(b)所示为过渡分支时的尾涡结构，此时呈现“2S”；图10(c)所示为高频高振幅分支下的尾涡结构，此时涡脱模式为“2S”，呈现出明显的卡门涡街规律性，漩涡脱落分离点固定于来流棱角处，且尾涡强度较强，这也导致了高振幅分支的大幅振动。

4 结　语

本文对不同来流角度下三棱柱的流致振动响应特性进行了数值模拟研究，系统分析了来流角度对柱体流致振动响应的影响，揭示了截面几何特征与来流角度对柱体振动响应、流体力及相位变化、尾流结构的影响规律，深入探索了三棱柱的振动激励机制，为后续非圆截面流致振动研究提供重要参考。本文主要结论如下：

1）不同截面形状柱体的流致振动激励机制存在差异。圆柱的流致振动为典型的涡激振动，表现出明显的“自限性”特征，三棱柱由于存在几何棱角，导致结构与分离剪切层作用更为复杂，当来流角度0°时三棱柱由于攻角变化导致升力失稳是其诱发驰振的主要原因，涡激振动是三棱柱诱发驰振响应的前提。

2）来流角度是影响三棱柱流致振动响应和振动分支模式的重要因素。当来流角度为0°时，三棱柱发生驰振响应；来流角度为30°、45°、90°时三棱柱发生高频高振幅的振动响应，来流角度为60°时出现类似涡激振动的振动响应。

3）三棱柱的几何棱角在漩涡脱落中起了重要作用，从而使三棱柱呈多种振动模式：来流角度0°时，三棱柱后延尖角的存在，使漩涡能从分离剪切层中快速切断，没有形成长涡街，整个尾涡区斜向偏移，从而导致三棱柱呈低频振动；来流角度30°时，漩涡脱落分离点固定于来流棱角处，且尾涡强度较强，从而导致高振幅分支的大幅振动。

传统的CFD方法依赖于大量计算资源和高精度数值求解，常面临计算耗时长、复杂流场难以快速预测等挑战。AI为CFD提供了新解决方案。AI在提高CFD仿真效率、提高建模精度及实现仿真的智能化和自动化体现出了巨大潜力，从而加速CFD在海洋工程、航空航天、能源、环境等领域的应用与发展。