环境的枯竭和日益增长的能源需求形成了一个尖锐的问题，而海洋覆盖了地球的绝大部分区域并蕴藏着丰富的能源。当海流流过钝体时会在其后脱落旋涡，并产生周期性的升力，发生流致振动(flow induced vibration)现象^[1]。这种振动往往被认为是一种会破坏结构强度并造成疲劳损坏的有害现象，并且它广泛存在于海洋结构物中，例如油气开采平台、系泊系统、海洋锚链、立管甚至是桥梁建筑。然而另一方面，研究也发现在产生振幅的同时，结构也会吸收大部分流体动能。相较于涡轮机以及波浪能发电，利用流致振动获取能量具有显著优势^[2]：①适用范围更广，在低流速海域即可产生较大的振幅；②所获取的能量是完全清洁的可再生能源；③能量获取系统的设备危险系数低，对海洋游鱼等生态造成的影响小；④整个设备可完全工作于海面以下，不会对水面舰船造成影响。因此，对于流致振动响应及其能量获取特性的研究有相当的价值。

流致振动的主要形式包括涡激振动(vortex induced vibration)和驰振(galloping)，Williamson团队对圆柱涡激振动进行了大量研究，他们将其响应分为初始分支、上分支和下分支，并将脱涡模式总结为S和P模式，其中S代表单个涡，P代表脱落一对涡^[3-4]。此外，Sarpkaya^[1]和Bearman等^[5]也对圆柱涡激振动进行了总结与综述，其中包含了更多相关参数和细节。不同于涡激振动，驰振是一种高振幅、低频的升力失稳振动现象，由于其种种特性而逐渐被应用于流致振动能量获取系统^[6]。对于光滑圆柱，驰振并不会发生，而通过PTC(passive turbulence control)方法，如改变圆柱表面粗糙度，则可以使其出现驰振现象，并具备一定的能量获取能力^{[2, 6]}。此外，非圆截面柱体则可能在高折合速度(U_r)下发生驰振(折合速度是基于速度定义的无量纲参数，详见1.1)。方形截面因此是研究和应用驰振的理想模型，这不仅是因为其是工程应用中的常见形状，也因为其具有广为人知的驰振不稳定性。Parkinson和Smith^[7]运用quasi-steady方法对方柱驰振建立了数学模型，他们认为造成驰振发生的原因是方形截面对于瞬时来流的不对称性。同时，他们运用此模型成功预测了驰振不稳定性的起始点。文献[8-9]对方柱流致振动进行了实验和数值研究，在低雷诺数(Re)的数值模拟中未观察到驰振；而在高雷诺数的实验中，当攻角为0°时捕捉到驰振现象。徐枫和方平治等^[10-11]对在空气来流下的方柱涡激振动进行了数值模拟，重点分析了流体特性，振动响应和脱涡模式，然而由于这些研究所选取得折合速度不够广泛并未捕捉到驰振现象，此外，目前针对以水为介质的方柱流致振动以及其相应的能量获取系统的相关研究还不充足。因此本文以方柱为对象展开研究，给出了其在不同雷诺数(折合速度)下的振动响应，频率响应，水动力特性，脱涡模式以及能量收集效率。此外，本文也分析了方柱的振动形式，以期为相关海洋工程设计提供参考。

1 几何模型和数值方法 1.1 物理模型与计算域

如图 1所示, 流致振动能量获取系统可以被简化为质量—弹簧—阻尼模型, 其主要由方型柱体、弹簧、传动机构和发电机组成。方柱由弹簧支撑并在发生振动时驱动发电机工作。表 1列出能量获取系统的主要参数, D代表方柱的特征长度, 即边长; K代表弹簧单位刚度; M^*代表方型振子的质量比, 定义为方柱在空气中的质量(M)除以其排水质量; ζ代表阻尼系数, 定义为, 式中C为系统的总阻尼, 它包括耗散阻尼和有效阻尼, 由于难以确定具体的耗散阻尼数值, 因此在本文中忽略了其影响; 表中U_r代表折合速度, 定义为U_r=U/f_nD, 式中U为来流速度, 而f_n为系统的自然(固有)频率, 它可以表示为。值得一提的是, 有关固有频率f_n的设定还存在一定争议, 不过f_n只与后处理有关, 而不会影响计算结果, 因而在本文后续的计算中也使用了此种处理方法。另外, 由于在大多数工程应用中, 海洋来流多处于湍流状态。因此, 本文将折合速度U_r选取为3~20, 其相应的雷诺数范围为Re=24 000~160 000。可以发现雷诺数与折合速度的比值为常数(Re/U_r=8 000), 这也说明了系统自然频率(f_n)一直保持不变。此外, 整个系统仅被允许在y方向上做单自由度运动, 这是因为对流致振动而言, x方向的振幅往往远小于横向振幅, 且单自由度的运动更方便于能量传递与捕获。

图 2和图 3分别是数值模拟的计算域以及相应的网格示意图。计算域长45D, 宽30D, 方柱的中心距入口边界和上下边界均为15D。入口条件设为均匀流速度入口, 出口边界为完全出流条件, 上下边界采用滑移壁面。为了提高模拟的准确性以及计算效率, 计算域被划分为5个部分并独立生成网格, 其中运动域(1)中的网格经过了进一步加密, 并根据所选湍流模型(SST k-ω)而使用函数(y⁺≈1)来确定方柱边界层的第一层网格高度。当流致振动发生时, 运动域(1)会以铺层方式跟随方柱一起振荡以保证其中的网格不会更新及变化。这种分块动网格的方法可以有效避免网格的拉伸与畸变, 相比重叠网格减少了信息交换量, 保证计算精度的同时一定程度降低计算机消耗。

1.2 流体控制方程

流体域的控制方程采用二维unsteady Reynolds averaged Navier-Stokes(URANS)方程, 其表达式如下

(1)

式中:x_i代表位置向量, 而u_i代表相对应的平均速度向量; ρ, p和ν分别为流体密度, 平均压力和运动粘性系数; 基于Boussinesq假设, 可以将(1)式动量方程中的雷诺应力项写为

(2)

式中, v_t, k和δ_ij分别代表了湍动黏性系数、湍动能和克罗内克符号。

本文采用SST k-ω湍流模型对上述方程进行封闭。湍动黏性系数v_t可以通过(3)式求得, 其中a₁为常数等于0.31, 而ω, S, F₂为比耗散率、应变率和第二混合函数。

(3)

SST k-ω是一种两方程模型, 它包含了有关湍动能(k)和比耗散率(ω)的控制方程, 具体可以写为

(4)

(5)

式中:F₁为第一混合函数; P_k为动能项; α, σ_k, σ_w, β_k, β为常数, 其相应的具体表达和数值详见Menter^[12]。

一般当雷诺数大于300时, 流动会进入湍流阶段, 而此时认为流体流过钝体后所脱涡的旋涡有强烈的三维效应。然而, 不同于绕流情形, 当钝体发生振荡时, 其展向的相关性会大幅增加, 这使得流场表现出二维的脱涡模式, 并可以通过二维模拟预测流场的宏观变量^[13-14]。因此二维RANS方程对于预测涡激振动是合理的, 其具有宏观统计优势以及可以节约计算资源, 这种方法也被用在许多已出版文献中^[15-16]。

1.3 结构振动方程

用经典的质量—弹簧—阻尼振子模型描述方柱振动的动力学特性。在单自由度y方向上的运动由二阶线性方程建模如下

(6)

式中：F_y代表了方柱在垂直来流方向所受的升力; A_y表示系统在y方向上的位移。

本文用C语言编写了Newmark-β算法, 并将其嵌入到了用户定义函数(UDF)中以求解振动方程(6)。该数值方法形式如下

(7)

式中：分别是t时刻的横向位移, 速度和加速度; 参数α和β选取为0.5和0.25可以保证积分精度和稳定性。在选取了合适的参数后可以计算如下积分常数

(8)

将(7)式和(8)式相结合后, 可以得到在t+Δt时刻的速度及加速度

(9)

式中仍有3个未知量(和)需要求解, 将其引入到t+Δt时刻的质量—弹簧—阻尼方程中进行约化求解, 最终可以得到如下方程

(10)

通过此方法可以获得系统的振动响应并结合动网格技术完成网格更新。

1.4 能量收集方程

在方柱振动的一个周期T中, 其能量获取功率可以表示为

(11)

基于所采用的求解方法, 可以在每一个时间步利用UDF程序提取方柱所受到的流体力和速度, 并通过基于C语言的二阶精度求积方法对方程(11)进行直接求解。事实上, 在振动过程中刚度和质量项并不会对能量转换起到作用, 唯一的非零项仅为阻尼项。流体扫过方柱时的能量可以根据伯努利方程表示为

(12)

式中，L为柱体长度, 因此可以得到整个系统的能量转化效率, 它被定义为:η=100%×P_energy/P_fluid。

2 数值模型验证

数值模型会受到网格数量以及时间步长的影响, 因此, 在开始计算前需要考虑两者的独立性。本文选择了4种不同密度网格模拟雷诺数为95 000的方柱绕流, 并将结果列在了表 2之中, 其中升力系数定义为C₁=2F_y/ρU²DL, 而阻力系数定义为C_d=2F_x/ρU²DL, 表中升力系数的数值是通过计算其均方根并乘以2^0.5得到。可以看出, 当网格数量超过65 000时, C₁和C_d都保持了良好的一致性。综合计算效率和准确性, 我们选择Mesh3作为本文的计算网格。为了选取合适的时间步长, 文章计算了在折合速度U_r=12时, 不同时间步下的方柱流致振动响应。图 4绘制了无量纲振幅(D^*)的时程曲线, D^*被定义为实际振幅(A_y)除以方柱特征长度(D)。经过对比可以发现当时间步小于0.003 s时, 振幅几乎不受时间步长的影响。可以认为, Δt=0.003 s可以确保计算准确性并节省计算资源。

为了保证数值方法的可靠性, 首先模拟了单个圆柱体的振动响应, 并将计算结果与已发表的数据^{[3, 16]}进行了比较。本次计算时所选择的主要参数与文献[3]中的实验一致。

如图 5a)所示, 本文基于快速傅里叶变换(FFT), 绘制了频率响应的曲线。图中2条虚线A和B分别表示固有频率和脱涡频率。可以看出, 数值模拟结果可以清楚地描述这2种频率对流致振动的影响, 在锁定区域内脱涡频率(B)和振动频率(f_osc)均锁定在固有频率(A)附近。图 5b)给出了无量纲振幅响应的结果, 图中数值模拟曲线呈现出明显的初始分支, 上分支和下分支的变化, 这与实验结果及其他已出版的CFD结果类似。因此, 可以认为, 本文所采用的数值方法可以较好地模拟出流致振动的重要特征。

3 结果与分析

图 6给出了方柱的无量纲振幅随折合速度的变化曲线, 为了方便比较, 图中同时给出了单个圆柱的流致振动振幅响应。总体而言, 方柱的振动明显更为剧烈, 其振幅随折合速度的增加而不断增长, 在U_r=20时达到最大值4.153D。另一方面, 与圆柱不同, 方柱的振幅变化并没有出现明显的初始分支, 上分支和下分支, 这也与Cui^[15]的结果相一致。

图 7绘制了方柱和圆柱的频率响应, 图中的振动频率(f_osc)是无量纲振幅谱的峰值频率, 在f_osc随时间变化且有多个峰值的情况下, 则选择最主要的峰值频率。对于单个圆柱的流致振动而言, 在锁定区之外, 振动频率和涡流脱落频率遵循Strouhal定律, 即f_osc是U_r的线性函数。当共振发生时, 脱涡频率和振动频率会锁定在一起, 并靠近固有(自然)频率。而方柱的频率响应明显不同, 其振动几乎不受到脱涡频率的控制且主频一直低于固有频率, 这也导致了大幅振动的发生。结合图 6和图 7可以发现, 当雷诺数大于120 000时, 方柱的振幅突然大幅增加且一直保持增长趋势, 与此同时其频率却开始下降并不再随折合速度或雷诺数的变化而变化。这种现象表征方柱此时从涡激振动转为驰振。

对于涡激振动, 其共振区域发生在折合速度为U_v=1/St附近, 其中St为Strouhal数。对于本文所选取雷诺数范围, St可大约等于0.13。另一方面, 根据quasi-steady理论, 驰振不稳定性则发生在折合速度大于U_g=8πM^*ζ/A₁时, 其中A₁是一个与雷诺数, 截面形状, 湍流强度有关的函数。参数A₁可以通过测量绕流下的升力系数C_L0与来流攻角α在α=0°处的斜率获得(A₁=∂_L0/∂α|_α=0°)。在本文中, A₁的取值可近似为2.69, 与文献[7]一致。U_g和U_v的关系决定了方形柱体的振动形式:对于U_g < U_v, U_g=U_v或U_g略大于U_v的情况, 方柱会发生耦合的涡激振动与驰振, 且两者有强相互作用; 相反, 当U_gU_v时, 涡激振动与驰振的相互影响可以忽略, 而两者也会呈现分开的振动形式。针对本文所选取的阻尼系数, 质量比和雷诺数, 可以发现U_g/U_v的值远小于1, 这意味着方柱的流致振动存在强烈的涡激振动与驰振相互作用。可以认为驰振持续影响了涡激振动的振幅和锁定区域, 同时涡激振动也改变了驰振的振动形式。因此, 如图 6和7所示, 方柱并未如圆柱一样呈现出典型的锁定区以及分支特性, 而且直到U_r>18后, 驰振才脱离涡激振动影响展现出典型的高振幅低频特性。以上现象也表明, 本文所用数值方法可以成功捕捉涡激振动与驰振的相互作用。

图 8绘制了方柱在不同折合速度下的无量纲振幅时程曲线, 为了观察方柱位移与其所受流体力的关系, 图中同时给出了升力系数的时程曲线。在低雷诺数下(见图 8a)), 方柱呈现出不规律振动且振幅较小, 而当折合速度大于10且未发生驰振时, 如图 8b)所示, 方柱的振动开始呈现出明显的周期性波动。在驰振发生前(见图 8a)和8b)), 可以发现方柱的位移曲线与升力曲线保持了良好的一致性而并没有发生“相变”现象。然而, 当方柱在由涡激振动转换为驰振的过程中(如图 8c)所示), 其振动再次出现明显的不稳定性, 与此同时, 在振动位移波动一个周期的过程中, 升力系数曲线则发生多个周期震荡。在方柱进入完全驰振状态(Re>144 000, U_r>18)后, 其振动又恢复为规律波动(见图 8d))。

对位移曲线和升力系数曲线进行FFT变换后得到了图 9的频谱图。在图 9a)中可以发现无论是位移谱还是升力谱都有多个峰值频率, 且两者保持了较好的一致性, 其主频处在接近固有频率的位置。当驰振发生时, 位移谱频率主频降低, 而升力谱频率更为明显的受到脱涡频率的控制, 在图 9c)中, 方柱进入完全驰振状态时, 可以看出在位移谱中低频所占的比重增大, 这也是驰振充分发展的一个重要信号。

本文通过观察方柱振子周围的瞬时涡量等值线来研究涡脱落模式, 其中涡量定义为∂ω=∂v/∂x-∂u/∂y。图 10给出了方柱在Re=40 000, 88 000, 144 000时的瞬时涡量图, 压力云图以及速度流线图。“下”代表方柱处于位移最低处, “上”代表方柱位于最高处。可以观察到, 在雷诺数Re=40 000时方柱呈现出典型的S脱涡模式, 即在振动周期中正涡和负涡单独脱出, 因而在远场表现出经典的卡门涡街结构。当雷诺数增大时, 在方柱的尾涡中可以观察到P模式的脱涡方式, 如图 11所示, 在Re=88 000时, 方柱在每一个振荡周期会脱落共6个旋涡, 并以2对涡加2个单独涡的2P+2S方式脱落; 进入驰振状态后, 方柱在每个周期以6P+2S的方式脱落共14个旋涡, 其中2个单独涡分别在位移最高点和最低点脱落。综合图 10至12, 可以看到在最高点方柱脱落负涡, 而最低点脱落正涡, 这也说明并未发生“相变”。

图 13给出了方柱流致振动能量获取系统的转换效率以及功率, 可以发现, 在低雷诺数时, 由于方柱的振动幅度以及频率都较低, 因而其俘能功率以及效率均处于较低水平。然而当来流速度增大后, 方柱的能量转化功率P_energy和效率η也随之同步增加, 并在U_r=11时达到效率最大值7.156%, 发生这种变化的原因是结构振动频率逐渐靠近系统固有频率。进一步地, 当方柱由涡激振动向驰振转化时, 不论其功率还是效率都会突然降低。而进入完全驰振状态后, 功率P_energy又会大幅增加, 并且效率也保持相对稳定。这是由于在驰振转化过程中方柱振动频率降低导致功率和效率下降, 而在进入完全驰振后, 整个系统的振幅又会发生大幅增长进而又提高了发电功率。

可以从图 13中看出, 尽管在驰振过程中的能量转化功率P_energy远高于涡激振动区域, 但其效率却明显低于稳定的涡激振动区效率。造成这种现象的原因是能量转化效率η的分母P_fluid与来流速度U呈三次方关系(见公式(12)), 因此功率P_energy在由涡激振动向完全驰振转化的过程中其增长小于P_fluid的增长。此外, 相较于涡激振动, 方柱在驰振过程中的振动呈现出不稳定性, 且其振幅会随着来流速度的增加而不断增长。这意味着方柱流致振动能量获取系统所俘获的流体动能会随来流的增加而持续增长没有极限。因此, 可以看出方型柱体在应用于流致振动能量转化中具有明显优势, 尤其是在高雷诺数下。

4 结论

本文对方柱流致振动能量获取系统进行了数值模拟。通过URANS方程结合SST k-ω湍流模型获得流场信息, 通过Newmark-β算法求解振动方程。文章重点分析了不同雷诺数下的方柱振动响应, 频率响应, 水动力特性, 脱涡模式, 能量收集效率以及振动形式, 主要结论如下:①方柱的振动相较于圆柱更为剧烈且无明显分支, 其振幅随雷诺数增加而增大, 频率几乎不受脱涡频率控制; 方柱在U_r>14时由涡激振动向驰振转变, 并在U_r>18后进入完全驰振状态；②方柱位移曲线与升力曲线之间无“相变跳跃”现象, 在低雷诺数下两者呈不规律振荡而在高雷诺数下位移出现周期性波动; ③在低雷诺数下, 方柱以2S模式脱落旋涡, 而随着雷诺数增加可以观察到nP+2S的混合脱涡模式, 其中单个旋涡分别在系统达到位移最值点脱落；④方柱能量收集系统的效率和功率在低雷诺数时均不高, 随后在U_r=11时效率达到最大值; 在由涡激振动向驰振转化的过程中, 功率和效率均出现突然下降, 而在完全驰振状态功率开始大幅增加且随雷诺数增加没有尽头。后续研究可对方柱能量收集系统的阻尼、质量、刚度、表面粗糙度等参数进行优化以提高能量获取效率。本研究对海能流发电, 海洋平台及立管等工程实践有借鉴作用。