风、浪、流复合作用下海上风力机基础波浪荷载分布特性研究

2022-05-31柯世堂王同光

船舶力学 2022年5期

王硕，柯世堂，张伟，李晔，王同光

（1.南京航空航天大学a.土木与机场工程系，b.江苏省风力机设计高技术研究重点实验室南京 210016；2.上海交通大学海洋工程国家重点实验室，上海 200240）

0 引言

海上风力机面临着强风、巨浪和急流等极端复杂海洋环境的严峻挑战，多种流场在交互运动中会出现风-浪复合、浪-流复合甚至是风-浪-流复合作用的复杂形态。此类环境导致的波浪荷载效应十分显著，会对风力机基础结构产生波致损伤[1-3]甚至破坏作用。然而，国内缺乏针对风-浪-流复杂环境下海上风力机的设计标准，国际上仅挪威船级社（DNV）给出了较为全面的海上风力机设计手册DNV-OS-J101[4]，现有风-浪-流复合作用下海上风力机设计理论和技术储备相对匮乏。因此，准确描述复杂海洋环境中波浪演化规律及风力机基础结构在其作用下的波浪荷载分布是海上大型风力机设计亟待解决的重要难题之一。

针对大型海上风力机基础结构波浪力研究，主要集中在波浪荷载[5-7]、结构响应[8-10]、疲劳可靠性[11-13]等方面，上述研究主要采用数值模拟方法，相关成果很好地指导了大型海上风力机结构抗波浪设计。然而，已有方法大多对边界条件和材料属性进行了简化，例如传统的风-浪复合造波法，在施加风场时忽略了水质点在入口边界速度的反复性，导致波高出现反常的衰减[7,14-15]，且现有文献鲜有涉及风力机基础结构在风-浪-流复合场作用下波浪荷载沿水深与周向的分布特性研究。

鉴于此，以东海兴化湾6 MW 大型海上风力机基础为研究背景，采用Fluent构建数值波浪水池，首先改进传统风-浪复合造波边界自由液面处风与水质点速度方向的一致性，并在风-浪场中进行对比验证；其次描述风、流与浪复合场中波浪的时空演化规律；在此基础上，研究各工况中基础柱整体水平波浪力和根部弯矩的复合特性，引入并确定极值荷载复合系数；最后系统分析各工况基础柱波浪荷载沿水深与周向分布特性。

1 数值水池模型构建

1.1 控制方程

本文在数值波浪水池建立中，水和空气均采用不可压缩、粘性流体假设。取水平方向为x轴（向右为正），垂直方向为z轴（向上为正），则三维Navier-Stokes连续方程为

为防止波形反射干扰流场，在数值水池尾部添加动量源阻尼以实现消波[16]，添加阻尼项后的动量方程为

式中，u、v、w分别为x、y、z方向速度分量，t为时间，P为流体压力，ρ为流体的密度，β为动力学粘性系数，g为重力加速度，μ（x）为消波系数，x1与x2分别为消波区在x轴上始末位置坐标，α为经验系数（取为8）。

波面在传播过程中时刻发生变化，采用VOF（volume of fluid）方法对自由液面进行追踪处理，即通过引入相体积分数变量aq界定气液两相。当aq=0时，表示网格单元内只含气相；当aq=1时，表示网格单元内只含水相；当0＜aq＜1时，表示该网格单元处于气-液交界面。自由液面处aq应满足

式中：V为流体速度矢量；aq为控制q相流体在网格容积中所占体积分数，当q=1、2 时，aq分别表示气、水相体积分数。

1.2 复合造波方法

数值波浪的模拟方法主要有仿物理造波、源项造波和边界造波，其中仿物理造波依赖动网格技术，带来网格数量急剧上升和计算效率下降；源项造波受到造波源的扰动干扰，波形稳定性相对较差，且对波浪品质影响较高的造波源位置不好控制；而边界造波法适用于多场复合模型同时具有较高的模拟精度和计算效率，在海洋工程数值波浪水池领域的应用更为广泛。与线性波相比，二阶Stokes波存在质点漂移和质量前移，更为接近实际海况，其速度势函数φ（x，z，t）为

波面方程η（x，t）为

波浪水平向流速uwx为

波浪垂直向流速uwz为

式中，H、T、L、d分别为波浪的波高、周期、波长和静水深度，k为波数（k=2π/L），ω为波浪圆频率（ω=2π/T）。

风、流与浪复合流场的入口边界包括自由液面以下水质点速度和上部空气速度，其中水质点速度可由波浪流速和海流流速两部分组成，可根据波流共同作用下流速场理论进行叠加[17]，即

式中，ux为复合场中水质点水平速度，uwx为波浪水平速度，uc为海流水平速度。

自由液面以上的空气流速可采用A类地貌指数率梯度风进行定义，即

式中，η为波面方程，U10为基本风速，z10为风场参考高度（取10 m），α为粗糙度指数（取0.12）。

通过用户自定义函数（user defined functions）组合风速剖面、入射波速、海流流速可构建浪-流、风-浪及风-浪-流复合流场，图1给出了三维风-浪-流复合数值水池模型示意图。水池尺寸为180 m×60 m×20 m（x×y×z），其中水域高度为15 m，空区域高度为5 m，消波区尺寸为60 m×60 m×1.5 m（x×y×z），波浪沿x轴正方向传播。本文选用RNGk-ε湍流模型，基于三维双精度分离式求解器，采用PISO（pressure-implicit with splitting of operators）算法实现压力与动量的解耦。压力方程选用加权体积力格式（body force weighted），自由面重构方法采用Modified HRIC（high resolution interface capturing），湍流动能、湍流耗散率与动量方程均采用二阶迎风格式，控制方程的计算残差设为1×10-6，压力参考值为一个标准大气压，时间步长为0.005 s。边界条件设置为：左边界为速度入口，右边界为静水压力出口，水池底部为无滑移壁面边界，顶部与两侧面为对称边界。

图1 三维风-浪-流复合数值水池模型示意图Fig.1 Schematic diagram of three-dimensional wind-wave-current combined numerical pool model

2 流场模拟及分析

2.1 工程概况及网格划分

以东海兴化湾6 MW大型海上风力机三桩基础为例。基础柱直径为3.3 m，柱距为40 m，三个基础柱采用等边三角形布置，高为18 m，其中静水面以下15 m，水面以上3 m。波形参数选自40海区，波浪周期T=5 s 对应最大频率波高为H=1.5 m，据弥散关系[18]可求得波长L=39.15 m。基本风速U10取为30 m/s，根据波面风海流估算公式uc=0.02U10，取剪切海流流速uc=0.04zm/s（z为高度坐标），风力机基础风浪流工况划分如表1所示。

表1 风力机基础风浪流工况划分表Tab.1 Conditions division of wind turbine foundations under action of wind,wave and current

波浪属于典型的非定常流动问题，模拟精度对网格分辨率的依赖程度极大。为确保数值模拟结果的可信度，本文增加了网格无关性验证。由表2可知，随着网格数的增加，网格质量逐渐提高，波高误差和波高标准差呈现逐渐减小趋势。方案四和方案五的网格质量与模拟结果无明显差异。综合考虑计算精度与效率的平衡，本文选取了网格总数为490万的方案。

表2 不同网格数量下数值波与理论波对比Tab.2 Comparison between numerical wave and theoretical wave under different grid numbers

图2给出了数值水池网格划分示意图，采用混合网格，将整个计算域分为内外两个部分。基础柱核心加密区采用四面体非结构网格，基础柱壁面第一层网格高度取为0.002 m，渐变率为1.1。外围区域较为规整，采用六面体结构化网格，网格分辨率取为Δx=L/120，Δy=H/5，Δz=H/20。为保证波面捕捉和波浪荷载计算精度，对入口边界、基础柱周围和自由液面附近网格进行加密。消波区网格尺寸逐渐稀疏，减少网格数量的同时增强数值耗散，强化消波效果。为了下文研究波浪荷载分布特性，提取不同位面上的波浪荷载，在网格划分过程中将柱表面区域进行分割，如图2（b）所示。

图2 数值水池网格划分示意图Fig.2 Schematic diagram of wind turbine foundation grid division

2.2 模型验证

为验证二阶Stokes 波模拟结果的可靠性，图3 给出了的数值波与理论波对比。由图可知，水体由静态到动态平滑过渡，波形稳定后数值波与理论波形锲合度较高，验证了造波程序的正确性。

图3 数值波与理论波对比Fig.3 Comparison of numerical wave and theoretical wave

浪-流复合造波已有广泛的理论支撑[17,19]，为验证风-浪复合造波的可靠性，图4给出了传统风-浪复合场波形变化。由图可得，风场的加入使波速增大，相比纯波浪场更早达到稳定状态，然而波峰波谷水位均有降低，波浪振幅略有减小，这有悖于风浪发育成长的运动结果[20]。

图4 传统风-浪复合场波形变化示意图Fig.4 Schematic diagram of waveform change of traditional wind-wave combined field

为探究传统风-浪复合造波法导致波浪振幅衰减的原因，图5 给出了传统风-浪复合场入口边界附近速度矢量云图。由图可知：波浪水质点的速度方向存在周期性的正负交替，倘若在造波边界波面以上区域加入定向风，当风与波浪质点速度同向时，高速风场对波浪输入能量，风浪成长并出现一定倾角；当风与水质点速度逆向时，正向风被逆向水体阻挡，两者相遇结果表现为近水面风向发生折射，流场发生紊乱，局部水体出现风致凹陷，此时风与波浪之间存在能量抵消并导致波浪振幅衰减。

图5 传统风-浪复合场速度矢量云图Fig.5 Velocity vector contour of traditional wind-wave combined field

2.3 风-浪复合造波改进

本文对传统风-浪复合技术进行改进。定义波面以上波峰以下区域的用户自定义函数，将波浪水质点速度化作单位矢量加入风剖面公式之中，从而保证气液交界面处风与水质点速度方向的一致性，避免风与水质点速度相反时造成波能损失，更加符合真实风-浪复合时的能量传递规律，同时可以获得更好的收敛性。流场入口边界水平方向速度计算公式为

图6 给出了两种风-浪复合模拟结果对比示意图。对比分析可知，改进前后风-浪复合场波幅差异明显。具体表现为波峰抬升、波谷下沉，波浪起伏增大约15%，体现出风对波浪的能量传递。改进后的风-浪复合场在长时间的演化计算中能够保持质量和动量守恒，未出现波形突变。

图6 两种风-浪复合模拟结果对比Fig.6 Comparison of the two wind-wave combined simulation results

图7给出了风-浪复合场压力分布对比。对比分析可知，传统风-浪复合场中无明显的压差区域，而改进后风-浪复合场波峰两侧压力不对称，迎风面压力大于背风面，导致压差环境形成，风浪得以成长，这符合风浪发育成型的“遮拦机制”[21]，也验证了本文风-浪复合改进方法的有效性。

图7 风-浪复合场压力分布对比Fig.7 Comparison of pressure distribution in wind-wave combined field

3 波浪传播时空演化规律

3.1 时程演化

图8 给出了复合场中波形变化时程。由图可知：海流的加入使波浪场的波高降低约50%，随流场的演化在30 s 时出现波群现象，波群周期约为60 s；在波群前半周期水位上升而后半周期水位下降；风与浪-流场的复合效应相比波浪场更加显著，波群出现后会改变其上方的风场结构，从而致使大量能量由风向海浪转移[22]，导致浪-流场的波峰抬升、波谷下沉，波高增大接近一倍；剪切海流和风场的加入对波浪周期的影响很小。

图8 复合场中波形变化时程Fig.8 Time course of waveform changes in combined fields

3.2 空间演化

图9给出了各工况流线空间分布示意图。由图可知：工况1流线由多组“U+V”型区域间隔排序组成，其中波谷下部为顺时针旋转U型区，波峰下部为逆时针旋转V型区。工况3中，波面以上空气回流被正向风场破坏，风掠过波峰之后在波谷区域积聚旋转，风能向波浪能转移导致U 型区增大、V 型区减小；海流加入后，波浪场的流线结构被破坏，水质点保持正向运动导致U 型区消失，工况2中上部流线以波峰为中心逆时针旋转，下部流线在水平方向上平直；风作用于浪-流场后，波峰下部V型区与流线曲率同时增大，说明下部水体受风场的影响运动变得更加剧烈。

图9 各工况流场线空间分布Fig.9 Spatial distribution of flow field line in each working condition

4 波浪荷载模拟结果分析

4.1 荷载验证

当基础柱直径/波长＜0.2 时，基础柱对波浪场的影响微弱，可不考虑绕射力，根据《港口与航道水文规范》[23]用Morison方程计算基础柱的波浪荷载。单位长度上基础柱水平波浪力计算公式为

式中，CD为拖曳力系数（取1.2），CM为惯性力系数（取2.0），D为基础柱直径，A为柱体截面面积，a为水质点轨道运动的水平加速度。

规范规定当桩距大于4 倍桩径时，不考虑群桩效应。以距离入口最近的基础柱为例，模拟值与规范值对比如图10 所示。由图可知，模拟值与规范值周期相位变化基本吻合，水平波浪力正负峰值略小于规范值且误差不超10%，再次验证本文数值模拟的有效性。

图10 波浪场中基础柱水平波浪力对比Fig.10 Horizontal wave force comparison of foundation column in wave field

4.2 基础整柱荷载分析

图11给出了各工况基础柱波浪荷载时程曲线。由图可知：水平波浪力与根部弯矩正负峰值同步出现，水平总波浪力作用点（M/F）在海床以上2/3 水深附近；海流对波浪荷载的贡献较小，浪-流场的波浪荷载总体上大于波浪场和浪、流叠加值，伴随波群发展个别波周期波浪荷载出现陡降；风-浪复合结果表现为正峰值增大约7%～11%，负峰值增大约4%；风-浪-流复合作用下，波浪荷载正峰值相对于风-浪和浪-流场又进一步增大，且最大增幅超过25%；三场复合荷载比浪、流叠加值增大约10%～35%，鉴于此，建议在结构设计过程中综合考虑不同的复合工况，对荷载的安全系数予以放大。

图11 各工况基础柱波浪荷载时程曲线Fig.11 Time history of wave load of foundation column under various working conditions

为探究流场复合对荷载的放大效应，对比分析复合场作用下基础柱受力与各单场作用线性叠加的关系，定义每个周期的极值荷载复合系数γ为

式中，P复合为复合场每周期水平波浪力极值，∑P单场为对应周期单流场水平波浪力叠加值。

图12给出了各工况基础柱极值荷载复合系数时程变化曲线，其中横坐标为基础柱经历波浪冲击的波周期个数。由图可得，各复合场的极值荷载复合系数γ普遍大于1，风、流与浪之间的复合作用总体对荷载呈现放大效应，其中风-浪-流三场复合的放大效应最大达135%。这主要因为：

图12 极值荷载复合系数时程曲线Fig.12 Time history of the combined coefficient of the extreme load

（1）工况2 中，水质点正方向的运动速度沿水深方向要比波浪场单独作用时大，导致拖曳力正峰值增大。U 型区消失后，流场丧失加速度缓冲过程，波谷相位需要更大的逆流向加速度来克服顺流向海流的影响，导致惯性力负峰值增大。

（2）风加入浪-流场后，三者在复合运动过程中，相对高速的风场对浪-流场不断输入能量，波浪表面的水质点速度随着风生流显著增大，同时风把波浪的波面高程抬升，增大波浪冲击时的受力面积，最终对波浪的作用力产生增益效果。

（3）从波能角度来看，一个波的总能量Ew=ρgH2L/8，其中H为波高，L为波长。风使波-浪场、浪-流场的波高和波长增大，波浪能变大，因此会出现极值荷载复合系数γ＞1的情形。

存在海流的复合流场受波群影响，水位降低使基础柱的受力面积减小导致基础整体波浪荷载减小。同时海流的存在会改变复合流场的速度与加速度时空分布特性，基础柱所受拖曳力、惯性力分量的空间分布特性随之改变，当复合场中拖曳力、惯性力分量沿水深方向的积分小于单流场叠加值时，会出现γ＜1的情形。

4.3 水深分布

为探究水平波浪力的水深分布，将基础柱从柱底到柱顶划分为9个分段，分别为Cy1-Cy9，每段2 m，前7段在液面以下，Cy8处于气液交界面处。图13给出了波浪场中水平波浪力沿水深分布，其中图13（b）横坐标a代表水平波浪力沿水深方向的比重。由图可知：水平波浪力沿水深方向呈指数型变化规律，上部水体占比较大，Cy8 未被水完全浸没，水平波浪力占比小于Cy7，Cy9 分段处于空气中，水平波浪力几乎为0；Cy1～Cy7分段水平波浪力正负峰值同步出现，Cy8分段荷载的正负峰值出现时刻相比其他分段存在相位差，表现为正峰值推迟出现，负峰值提前到来，由负到正的爬升曲线存在拐点；同时可以发现，水下各分段正负峰值并不相等，Cy1～Cy7分段水平波浪力负峰值略大于正峰值，而Cy8分段正峰值大于负峰值。

图13 波浪场中水平波浪力沿水深分布Fig.13 Horizontal wave force distribution along the water depth in the wave field

为探究波浪荷载沿水深方向的复合效应，采用式（13）进行无因次化分析，定义波浪力系数λ为

式中，F为各分段水平波浪力，Fmax、Fmin分别为波浪场中对应分段的水平波浪力正负最值。

图14给出了基础柱各分段波浪力系数时程曲线。由图可知：工况2与工况4中Cy1～Cy7相比工况3复合效应更显著，水平波浪力系数随水深增加而增大；风对波浪场荷载影响较小，而对浪-流场影响较大，主要表现为正负峰值均有增大且正峰值增益更为明显。位于自由液面的Cy8 复合效应最为强烈，风-浪复合导致荷载正峰值增大约20%；工况2中的波浪荷载极不稳定，波群周期前后波浪荷载相差接近一倍；风场加入浪-流场后，荷载正峰值增大而负峰值减小。

图14 基础柱各分段波浪力系数时程曲线Fig.14 Time history of wave force coefficient of each section of foundation column

4.4 周向分布

为探究自由液面附近波浪荷载的分布特性，选取Cy8分段基础柱，观察各工况不同波浪相位水平波浪力的周向分布。因Cy8分段并未完全浸入水中，波群现象会造成水位的上下波动，波浪荷载会有较大悬差，本文综合考虑了波群全周期，以0°为入流角，图15 与图16 分别给出了波浪周期相位和两个典型波浪周期的水平波浪力周向分布。由图可知：

图15 波浪周期相位图Fig.15 Periodic phase diagram of waves

图16 两个典型波周期水平波浪力周向分布Fig.16 Circumferential distribution of horizontal wave forces in two typical wave periods

波浪力在各工况中沿0°波向角呈现对称分布规律，正负极值分别出现在0°和160°，并在90°接近于0；基础柱迎浪面的波浪力始终保持为正值，最小值（接近于0）出现在T1 相位，最大值出现T5 相位，背浪面始终保持为负值，最小值（接近于0）出现在波谷（T2），最大值出现在波峰（T6）。

Cy8 在波谷（T2）周向波浪力接近于0，处于“卸载”状态；波峰（T6）水平波浪力接近于0，但周向分布的波浪力幅值较大，说明基础柱在波峰相位迎/背浪面两侧出现“加速度对流”，此时基础柱承受两侧波浪的惯性力挤压作用显著。同时可以发现，在波群前半周期，工况4波浪荷载的周向分布处于包络工况，在后半个周期，工况3 波浪荷载的周向分布包络其他流场，此类现象在T4～T7 相位表现更为明显。