半潜式浮式风机平台绕流力学特性与尾流分析

邓露 ，刘建祥

（1.工程结构损伤诊断湖南省重点实验室，湖南长沙410082；2.湖南大学土木工程学院，湖南长沙410082）

日益增长的能源需求促进了海上风电行业的发展，半潜式浮式风机平台适用水深范围广，是一种开发具有储量丰富、风速大等优点的深水风能资源的重要海洋平台[1].长期以来，柱体绕流一直是海洋工程领域的重点问题[2].流体流经浮式风机平台的柱体结构时，柱体后方周期性脱落的漩涡会导致平台大幅往复运动，这将加剧平台系泊结构的疲劳损伤，甚至导致系泊结构疲劳破坏[3].因此开展半潜式浮式风机平台绕流研究有着重要的工程意义.

在单柱绕流方面，国内外学者针对单柱绕流的力学特性、流场特性及其影响因素等已进行了大量的研究，并取得了丰硕的成果[4-7].在双柱和多柱式结构绕流方面，贾晓荷和刘桦[8]采用大涡模拟方法分析了并列或串列排布的双圆柱的绕流特性，结果表明串列排布时，上游圆柱对下游圆柱受力有影响，并列排布时，两柱的阻力系数均值基本相等.郑宇华和顾杰[9]的研究表明，间距比（双柱中心的间距与柱宽的比值）影响并列双柱开缝下方的漩涡个数、形态等.Kitagawa T 和Ohtab H[10]采用大涡方法模拟了串列双圆柱的绕流特性，结果表明当间距比小于某一临界值时，只有下游圆柱有涡脱现象.李聪洲等[11]研究了串列双柱绕流的力学特性，结果表明存在一个临界间距比使得圆柱的升力和阻力发生跳跃性变化.刘为民等[12]的研究表明，0°和 45°来流时，阵列四立柱柱群中立柱的涡脱频率大于单柱的涡脱频率，涡脱模式表现为2S 模式.刘正浩等[13]的研究表明，来流角对阵列四柱式海洋平台的立柱的升力和阻力系数影响较大，不同来流角下平台底部的流场不同.

双柱和多柱式结构绕流时，柱体尾流之间的相互干扰会导致立柱的力学特性和流场结构更复杂[11].然而，以上针对双柱和多柱式结构绕流的研究，主要分析了流速、来流角和间距比等因素对结构的受力、流场形态和涡脱模式等的影响，而深入分析流场结构以及尾流之间干扰机理的研究极少.因此，本文采用DDES 方法对不同来流角和流速下多柱式半潜式浮式风机平台的绕流进行了三维数值模拟，分析了平台立柱的阻力系数平均值、升力系数频率等力学指标，从流场的相干结构层面研究了尾流的干扰机理，并分析了尾流的空间相关性.研究结果阐释了多柱式结构绕流的尾流相互干扰机理，为深入理解多柱式结构绕流现象提供了理论依据.

1 数值建模

1.1 DDES 方法

Spalart 等[14]提出的 DES 方法（Detached Eddy Simulation）是一种雷诺时均法（Reynolds-averaged Navier-Stokes，RANS）和大涡模拟（Large Eddy Simulation，LES）相结合的方法，其基本处理原则是：在近壁区采用雷诺时均法进行模拟，在远离壁面的区域采用LES 来模拟分离流动.DES 方法通过长度尺度来实现RANS 到LES 的转变，其定义为：

式中：d 为离最近壁面的距离，Cdes取 0.65，Δmax为局部网格的最大尺寸.

然而，DES 方法具有局限性，在某些情况下，当RANS 转变到LES 时，会出现模型应力损耗现象，并导致“网格诱导分离”问题.因此，Spalart 等[15]对DES方法作了改进，提出了含有延迟函数的DDES 方法，解决了“网格诱导分离”问题.DDES 方法通过引入延迟函数对长度尺度的定义进行修正，修正后的定义为：

式中：fd为延迟函数，其表达式为：

式中：rd为延迟因子，其表达式为：

式中：κ 为卡门常数，取 0.41；ui，j为速度梯度；vt为动力黏性系数；v 为分子黏性系数.

DDES 方法是一种改进的RANS 与LES 相结合的方法，其在模拟流体分离流动时十分有效[16]，已得到广泛的应用.该方法具有多个优点，如不需要精细的网格也能得到较好的结果，适合复杂的结构等.因此，本文采用该方法进行研究.

1.2 平台模型

选取美国国家可再生能源实验室设计的半潜式浮式风机平台作为分析对象，平台的立柱（包括边柱和中柱）、连杆和浮箱底座的结构布置如图1 所示（图中括号内为构件编号）.平台的结构尺寸，吃水深度等设计参数详见文献[17].

图1 半潜式浮式风机平台模型图Fig.1 Model of semi-submersible floating wind turbine platform

1.3 计算域与网格划分

流场计算域及模型平面布置如图2 所示.θ 为来流角，即中柱与边柱2 的连线与x 轴负向的夹角.流域入口边界距离边柱2 中心为10D，出口边界距离边柱2 中心为25D，左右边界距离中柱中心均为10D.D 为边柱直径.

图2 计算域与平面布置图Fig.2 Computational domain and plane configuration

边界条件设定为：入口为速度入口，出口为出流边界，顶部为自由边界，其余边界为对称边界，平台表面采用无滑移固定壁面.

网格划分示意图如图3.由于平台几何结构复杂，因此，将流场分成两种类型的网格进行划分，在平台附近区域采用四面体网格以适应复杂的几何结构，在其他区域采用六面体网格以加快计算速度.为控制网格密度，将流域分成九个区域，在平台附近和尾流区划分密集网格，其他区域则划分较稀疏的网格.根据张新曙等[18]的建议，经过验算，并综合考虑计算耗时成本，确定网格方案为：加密区（平台四周1.4D 的区域）的网格大小为0.12D；尾流区的初始网格大小为0.25D，且随着尾流远离平台，网格逐渐稀疏.

图3 模型网格划分示意图Fig.3 Mesh generation of model

1.4 计算设置

流域采用Ansys Fluent 进行求解.时间离散采用二阶隐式格式，空间离散采用二阶迎风格式，压力与速度耦合方程采用PISO 方法（Pressure-Implicit with Splitting of Operators）.时间步长取0.05 s，每个时间步长内最大迭代次数为20，模拟时长为6 000 ～10 000 s，以保证获得稳定的结果.数据处理时，选取最后约0.6 倍模拟时长内较稳定的计算结果进行处理.

1.5 计算工况

本文选取 3 种典型来流角：θ=0°、30°、60°.根据中国南海的实际情况[19]，选取12 种来流速度：u=0.32 m/s、0.42 m/s、0.53 m/s、0.63 m/s、0.74 m/s、0.85 m/s、0.95 m/s、1.06 m/s、1.16 m/s、1.27 m/s、1.37 m/s、1.48 m/s，来流方向沿x 轴正向.

2 计算结果及分析

2.1 力学特性

流体流经半潜式浮式风机平台时，在水平面内，平台立柱受到两个力的作用，一个是顺流向的阻力Fx，又叫拖曳力，另一个是横流向的升力Fy.定义阻力系数为，升力系数为其中，u 为来流速度，ρ 为海水密度，取为1 025 kg/m3，A 为柱体迎流投影面积.

2.1.1 阻力系数平均值

图4 ～图7 分别为中柱、边柱1、边柱2、边柱3的阻力系数平均值随来流角和流速变化的关系图.从这四幅图中可知，在0°来流下，中柱的阻力系数平均值随流速的变化而小幅波动，其余立柱的阻力系数平均值随流速的改变而基本不变.在30°来流下，边柱2 和边柱3 的阻力系数平均值随流速的改变而基本不变，但中柱和边柱1 的阻力系数平均值随流速的改变而变化较大，离散性较强.这是因为流体流经上游立柱后发展成为湍流，湍流本身具有较强的随机波动性，所以，受上游尾流的影响，下游中柱和边柱1 的阻力系数平均值随流速的变化而波动较大[12].在60°来流下，中柱、边柱2 和边柱3 的阻力系数平均值随流速的改变而基本不变，而边柱1 的阻力系数平均值出现了跳跃波动的现象.这也是因为边柱1 处于下游位置，其阻力系数平均值受上游尾流的影响较大.

图4 中柱的阻力系数平均值Fig.4 Average drag force coefficient of middle column

图5 边柱1 的阻力系数平均值Fig.5 Average drag force coefficient of side column 1

图6 边柱2 的阻力系数平均值Fig.6 Average drag force coefficient of side column 2

图7 边柱3 的阻力系数平均值Fig.7 Average drag force coefficient of side column 3

由图4 可知，中柱的阻力系数平均值随来流角的增大而减小.阻力是漩涡脱落引起的立柱前后表面的压力差和立柱表面的摩擦力共同作用所形成的.随着来流角的改变，立柱和立柱、立柱和底座在垂直于来流方向的相对距离发生了改变.由图2 的平面布置图可知，在0°来流时，中柱后方尾流距离立柱和底座较远，下游柱体的阻挡作用小，上游中柱的尾流漩涡能够充分地发展，所以阻力系数最大.在30°来流时，边柱1 和底座1 靠近中柱尾流，下游柱体的阻挡作用增大，下游中柱后方的漩涡脱落不能充分发展，故阻力系数减小.在60°来流时，边柱1 和底座1位于中柱的正后方，完全阻挡了中柱尾流漩涡的发展，故阻力系数最小.由图5 可知，边柱1 的阻力系数平均值随来流角的增大呈现先减小后增大的趋势.由图6 和图7 可知，边柱2 和边柱3 的阻力系数平均值均随来流角的增大而增大.这些也均是因为随着来流角的改变，柱体之间在垂直于来流方向的相对距离发生了改变，柱体与柱体之间的相互影响所导致的.

2.1.2 升力系数时程曲线及频谱分析

限于篇幅，仅给出了 θ=0°、u=1.06 m/s 时各个立柱的升力系数时程曲线，如图8 所示.从图8 中可以看出，边柱2、边柱3 和中柱的升力系数时程变化表现出较好的谐波特性，而边柱1 的升力系数变化波动较大.边柱2 和边柱3 的升力系数变化相反，中柱的升力系数变化滞后于边柱2.当边柱2 的升力系数达到波峰时，中柱的升力系数稍有增大，且曲线变化趋势发生了改变.以图8 中t1-t2时间段为例，此时间段内中柱的升力系数未按图中粗划虚线趋势变化，而是按实线趋势变化.这可能是因为中柱处于边柱2 的尾流中，当上游边柱2 的升力系数达到波峰时，其后方的漩涡脱落，随后撞击到中柱上，导致中柱的升力系数稍微增大而改变了曲线的变化趋势.

图8 0°来流角下各柱升力系数时程曲线Fig.8 Time history curves of lift force coefficients of four columns under 0°heading

图9 ～图11 为u=1.06 m/s 时，三种来流角下平台各立柱的升力系数频谱图.立柱后方漩涡周期性地脱落是导致升力系数周期性变化的原因.由图9可知，0°来流角下，中柱的频谱除主峰外，还存在多个高频的小峰，小峰的幅度远小于主峰的幅度，说明中柱后方漩涡脱落以主峰频率为主，并伴随多个高频成分.边柱1 的频谱图中存在两个较大的峰值，低频对应的幅度显著大于高频对应的幅度，说明边柱1 后方低频的涡脱引起了较大的升力.边柱2 和边柱3 的频谱图中均只有一个峰值，说明这两立柱后方的的漩涡为单频脱落.此外，4 个立柱的频谱图中均有一峰值对应的频率为f=0.018 5 Hz.综上可知，平台立柱后方漩涡脱落的主频率为f=0.018 5 Hz，部分立柱还伴有高频或低频成分.

由图10 可知，30°来流角下，中柱、边柱 2 和边柱3 的升力系数频谱虽然有多个峰值，但对应的频率极低.受上游尾流的干扰，下游边柱1 的升力系数频谱没有明显的峰值.可见在此来流角下，平台立柱后方漩涡脱落不规则.

图9 0°来流角下各立柱升力系数频谱Fig.9 Lift coefficient spectrum of each column under 0°heading

图10 30°来流角下各立柱升力系数频谱Fig.10 Lift coefficient spectrum of each column under 30°heading

图11 60°来流角下各立柱升力系数频谱Fig.11 Lift coefficient spectrum of each column under 60°heading

由图11 可知，60°来流角下，中柱频谱没有峰值，这是因为流体绕过中柱后，撞击在下游的底座1和边柱1 上，没有形成周期性的涡脱.其余立柱的频谱均有多个峰值，但没有相同频率的成分.由此可见，平台立柱后方的漩涡呈多频脱落模式.其他流速下，各柱的升力系数频谱也具有类似的性质，限于篇幅，不再赘述.

2.2 流场分析

2.2.1 流场相干结构与尾流干扰

湍流场中有组织性的漩涡结构称之为相干结构，其在湍流的发展演化中起着重要的作用[20-21].u=1.06 m/s 时，0°、30°来流角下的涡量等值面图（按 x方向速度 u 着色）分别如图12（a）、（b），涡量等值面采用 Hunt 等[22]建议的 Q 判据（Q-criterion）表示.从图12 中可以看出，流体流经平台后，形成了多种形式和尺度不一的漩涡.图12（a）和（b）中的弧状涡包成月牙形分布在浮箱底座的底面和顶面边缘，反映出流体流经浮箱底座时，在底面和顶面边缘发生流动分离，并在随后的区域形成了马蹄涡.图12（a）中，边柱1、边柱2 和边柱3 后方有成对的、错位排布的管状涡包，说明这些立柱后方形成了交替脱落、错位排布的尾涡，尾涡流经中柱后，逐渐拉长而发展成为流向涡.从图12（b）中可以看出，边柱3 后方形成了斜向下的锥形涡包，说明该处漩涡结构的半径随着水深的增加而逐渐缩小.浮箱底座3 后方呈螺旋状上升的螺旋涡包，反映出其后方形成了发夹涡[21]，该漩涡结构由涡头和两条涡腿组成，涡头受到涡腿旋转向上的抬升力而上扬.靠近浮箱底座的发夹涡的尺度明显小于远离底座的发夹涡的尺度，说明随着湍流的发展，发夹涡逐渐发展壮大.

流体流经平台时，平台的立柱、连杆和浮箱底座的尾流之间会相互干扰.如图12（a）中所示，边柱2后方尺度较大的漩涡结构与下游斜杆2 后方尺度较小的漩涡结构相互融合，最终汇合成为一个整体，相

当i=j 时，表示该点的脉动量自相关，当i ≠j时，表示两点的脉动量互相关.相关系数的大小表征了脉动量在流场不同位置的关联程度，R=±1 表示完全相关，R=0 表示不相关.

图13 为0°来流角下y=0 平面结构的剖面图.为考察上下游立柱尾流间的空间相关性，选择x 方向速度u 作为脉动量，计算了y=0 平面，z=1.17D高度处边柱2（bz2）尾流点与中柱（zz）尾流点的相关系数，计算结果如图14 所示，其中，图14（b）为图14（a）的水平投影图.当于大漩涡“吞并”了小漩涡.如图12（b）中所示，浮箱底座3 后方形成的发夹涡在上扬的过程中，由于涡腿的旋转效应，诱导了边柱3 的尾涡向斜下方运动，最终与浮箱底座3 的尾涡汇合在一起，并导致了次生流向涡的形成.此外，由于该发夹涡涡腿的旋转效应，也诱导了中柱后方的流向涡发生偏转弯折，从而改变了中柱尾流的流动方向.由以上分析可知，尾流干扰的机理正是相干结构之间的相互作用.

2.2.2 空间相关性分析

流场内 i、j 两点的脉动量 Xi、Xj在空间上的相关性可用相关系数R 来表示，其定义[23]表示为：

图13 y=0 平面结构剖面图Fig.13 Stucture section at y=0

由图14 可知，边柱2 尾流与中柱尾流的x 方向速度在一定区间内具有较强的相关性.经计算，x 坐标在 [0.85D，1.1D]区段的边柱2 尾流与x 坐标在[3.1D，3.4D]区段的中柱尾流的x 方向速度的相关系数为0.4 ＜R ＜0.55，表现出较强的正相关性，而x 坐标在[1.7D，2D]区段的边柱2 尾流与x 坐标在[3.1D，3.4D]区段的中柱尾流的x 方向速度的相关系数为-0.55 ＜R ＜-0.4，表现出较强的负相关性.这说明虽然边柱2 后方交替脱落的尾涡流经中柱后，变成了流向涡，但两立柱的尾流之间仍有着较强的空间相关性.

图14 x 方向速度u 的相关系数分布图Fig.14 Correlation coefficient distribution of velocity u in x direction

3 结 论

本文通过 DDES 方法研究了 0°、30°和 60°来流下半潜式浮式风机平台绕流的力学特性，并对平台的尾流进行了深入分析，得出以下结论：

1）由于尾流的干扰，下游立柱的阻力系数平均值随流速变化而波动，离散性大，而上游立柱的阻力系数平均值变化不大.在0°来流时，平台立柱后方的漩涡脱落以某一主频率为主，并伴有其他高频或低频成分；60°来流时，平台立柱后方的漩涡脱落呈多频模式.

2）流场相干结构具有多尺度性，且类型多样，包括发夹涡、流向涡、次生流向涡、马蹄涡等.相干结构之间的相互作用是尾流干扰的内在原因.尾流干扰存在两种模式：处于上下游位置的立柱与连杆之间的尾涡相互交融，导致它们的尾流成为一个整体；浮箱后的发夹涡诱导立柱漩涡改变运动方向，导致立柱尾流汇入浮箱尾流或者运动发生偏转.

3）流体流经成串列排布的立柱时，即使流场的相干结构发生了转变，上下游立柱尾流之间仍有着较强的空间相关性.