海上风机单桩基础的直径与波长之比D/L<0.2，属于小尺度结构物，在计算波浪对单桩基础的作用力可以不考虑波浪反射和绕射的影响，一般使用Morison公式^［13-14］计算波浪载荷，波浪力可以由以下公式表示：

式中：ρ为水的密度（kg/m³）；D为柱体的直径（m）；d为水深（m）；η为波面；u(z)为水质点的速度（m/s）；∂u(z)∂t为水质点的加速度（m/s²）；CM为惯性力系数；CD为拖曳力系数。

线性波理论采用传统的Morison公式通常将波浪力从海底计算到静水面，但对于浅水区非线性波浪，这种计算方式的精度无法满足要求，而考虑波浪拉伸可以将任意时刻波浪力从海底计算到自由波面，实际工程将线性波理论与垂直拉伸^［15］、惠勒拉伸^［16］等经验公式相结合，从而更加真实的模拟波浪。

根据线性随机波理论^［17］，波面η(x,t)可以表示为：

式中：N为用于模拟的总的波列；ωi为每一个波的圆频率（rad/s）；εi为均匀分布的随机初相位；0<εi<2π，Ai为第i个波的振幅。

水平水质点在z处的速度ux为：

波浪拉伸中最简单的一种方式是垂直拉伸。在水线面以下的波的运动学依照公式（3）计算，但是水面线以上水质点的波的运动学均取在水线面点处的对应值，即：

使用惠勒拉伸方法，用以下的公式的等效节点高程z'来代替垂直坐标z，即：

式中：η(x,t)为在点x和时间t时的瞬时表面高程。

国内外学者对于破碎波问题开展了大量研究，目前，求解破碎波波浪力（以下简称破波力）的方法主要采用修正的Morison公式和冲击-绕流法。

美国海岸工程研究中心^［18］根据实验和Morison公式来计算最大破波力的公式如下：

式中：ρ为海水密度；D为桩柱的直径；Hb为破碎波波高。

当海底坡度i≤1/50时，可取Hb=0.78db，db为破碎波形成时的水深；当i>1/50时，根据IEC 61400—3^［19］和GL 2012^［20］建议破碎波波高Hb按照以下公式确定：

式中：T为波浪周期（s）；A和B为计算系数；g为重力加速度（m/s²）。

冲击-绕流法计算破波力是由Wiegel^［21］提出的，该方法将破碎波对圆柱的作用力分为两部分，一部分是由破碎波作用在圆柱上的瞬时冲击力FI，另一部分是由波浪对圆柱的绕流力，绕流力可按Morison公式进行计算，破波力可由以下公式计算：

式中：Fbreaking_wave为破碎波总力（N）；FD和FM为Morison公式拖曳力分量和惯性力分量（N）；FI为破波力的冲击力分量（N）。

波浪在圆柱上发生破碎时的相关参数见图1，C为波速，Hb为破碎波波高，R为圆柱半径，λ为卷曲系数，随圆柱相对于波向的倾斜度而增加，最大值取0.5，ηb为波面的最大高程。

图  1  用冲击-绕流法计算破波力

Figure 1.  The impact-flow method is used to calculate the breaking wave load

ηb的取值可参考《海港工程设计手册》中的公式：

式中：Hb为破碎波波高，可按照公式（7）计算；db为破碎波发生时的水深；L为波长。

Goda^［22］给出了最大冲击力的计算公式：

式中：ρ为海水密度（kg/m³）；ληb为冲击力的作用区域（m）；D为圆柱的直径（m）；Ub为破碎波波峰的运动速度（m/s）。

Ub可以按照以下公式计算：

本文以美国可再生能源实验室（NREL）开发的5 MW^［23］单桩式海上风机模型进行数值分析。该模型的整体参数见表1，单桩基础的外径为6 m，壁厚为0.06 m，水深为12 m，单桩基础与海底为刚性连接，风机在额定风速下为正常运行状态。

基于该5 MW单桩式风机模型在OpenFASTV 1.0版本中进行建模，OpenFAST采用空气动力-水动力-结构动力-控制系统耦合分析方法对海上风机进行全耦合时域仿真模拟，图2为OpenFAST各模块的耦合分析流程，具体的计算流程如下：

图  2  固定式海上风机耦合流程图^[24]

Figure 2.  Fixed bottom OWT coupling flow diagram

创建单桩式海上风机整体结构模型输入文件，OpenFAST主程序读入该输入文件，并对模型文件的有效性进行检查；空气动力学模块（AeroDyn）和水动力模块（HydroDyn）分别读入风速时程文件和海况文件，并由此确定边界条件和初始条件，基于FAST模型文件分别开展气弹性分析和水动力分析，得到气动力载荷和水动力载荷；依据得到的环境载荷，弹性动力模块（ElastoDyn）和基础结构模块（SubDyn）开展整体结构动力反应分析，同时控制策略模块（SeroDyn）依据风机运行的状态，来判断是否需要启动控制策略，从而得到整体结构动力响应参数，如位移、速度、加速度和内力等；FAST各个模块之间相互传递，互相作用，可实现海上风机的相互耦合作用。

破碎波载荷主要是来源于物理模型实验的实测数据。在物理模型试验中，图3为物理模型试验的布置图，整个风机系统采用1∶80的缩尺比，该模型利用浪高仪实测得到的波面历时曲线，利用弗洛德相似定律将实验得到的波面历时曲线转换成全尺寸模型，从而得到整个系统完整的波面历时曲线。

图  3  物理模型试验布置图

Figure 3.  Physical model test layout

本文使用HydroDyn模块中的Wavemod5选项将实验得到的破碎波数据导入到OpenFAST中计算破碎波，由于破碎波的波面时程是通过实验得到，并且在数值模型中考虑了波浪的拉伸，可以更好地模拟破碎波载荷，所以本文采用HydroDyn模块计算破碎波荷载具有合理性的。

通过将实验实测得到的破碎波波面高程作为输入导入到HydroDyn模块中。HydroDyn可以将波面高程的时间序列通过傅里叶变换，计算并存储每一个频率分量的振幅和相位，并使用这些替代的波浪谱和随机种子内部生成使用的水动力载荷和波的运动学^［24］。HydroDyn模块以过渡点处的浪高的时间历程开始计算，将波高转换成FAST可计算的格式对破碎波进行计算。具体的逻辑示意图如图4所示。

图  4  破碎波在FAST中模拟的逻辑图

Figure 4.  A logical diagram of the breaking wave simulation in FAST

本文选取50年一遇的东海海域实际海况环境，波浪海况选用的是非破碎波和破碎波，水深均为12 m，位于近海岸浅水水域，设定模拟时间1 000 s，设定风机的额定风速为11.4 m/s。模拟的工况如表2所示。

为对比垂直拉伸和惠勒拉伸的区别，在MATLAB里对惠勒拉伸和垂直拉伸进行编程来计算波的运动学，MATLAB程序计算波的运动学和水动力载荷是依照Morison公式进行计算。采用的是非破碎波，时间步长为0.01 s，计算工况采用工况1。

如图5所示，将这两种方法在水线面以上1 m所对应的水质点的速度进行对比，结果表明，惠勒拉伸方法的水质点速度低于垂直拉伸方法得到的相应值，两者的差异主要集中在波峰处，垂直拉伸得到的最大水质点速度3.62 m/s，惠勒拉伸得到的最大水质点速度是3.29 m/s，高于水线面1 m处水质点的速度垂直拉伸比惠勒拉伸更大。每个节点的最大水粒子速度分布的函数如图6所示，当节点为－10 m以上时候，垂直拉伸最大水质点的速度明显大于惠勒拉伸，而且在水线面处的两者的差距最大，高于水线面处的位置，两者的曲线会越来越接近，直到达到节点高程3.6 m处时候，两种方法计算得到的波粒子最大速度是一致的。

图  5  水线面1 m两种方法的水质点速度

Figure 5.  Wave particle velocities from two methods at 1 m above mean seawater level

图  6  两种不同方法的最大波粒子速度的对比

Figure 6.  Comparison of maximum water particle velocity from two different methods

FAST中计算垂直拉伸的运动学与MATLAB的对比，以非破碎波为例，验证两种程序计算垂直拉伸的波运动学的差异，从而为破碎波的计算提供依据。图7是不同节点处波浪水质点速度在FAST和MATLAB中的对比。从图中可以看出低于水线面1 m处FAST计算得到的水质点速度与MATLAB是一致的，水线面处和高于水线面1 m处水质点的速度，两种程序计算得到的结果也是一致的。从图中的对比情况来看，FAST使用垂直拉伸的波速在MATLAB中得到验证，可进一步验证数值模型的正确性。

图  7  水质点速度在FAST和MATLAB的对比

Figure 7.  Water particle velocities in FAST and MATLAB comparison

图8是水平波浪力的历时曲线，从图中可以看出水平波浪力在FAST和MATLAB中对比良好，曲线几乎一致，与MATLAB计算结果对比，FAST计算的水平波浪力得到验证。

图  8  水平波浪力

Figure 8.  Horizontal hydrodynamic load

本节主要模拟单桩式风机在非破碎波和破碎波作用下风机系统所受到的波浪载荷，验证在不同工况条件下波浪载荷的差异。波浪的方向为0°方向（X轴的正方向），非破碎波和破碎波模拟时长为1 000 s，为消除瞬态效应，去除前200 s的数据。在计算波浪载荷考虑使用垂直拉伸和无拉伸的情况，并将破碎波载荷与计算破碎波的两种方法（修正的Morison公式和冲击-绕流法）进行比较。

图9为不同工况下波面高程的对比，工况1和工况2的区别在于有无使用拉伸方式，而波浪的拉伸并不改变波面，所以工况1和工况2的曲线是重合的，工况3和工况4也是如此。由图8可知：非破碎波最大波高为7.2 m，破碎波最大波高为10.49 m，波浪破碎会引起波浪爬升。非破碎波的曲线十分平滑，破碎波波面峰值是非破碎波的两倍左右，并且波峰处很陡。由于单桩的存在会改变附近流场，单桩附近的波面高程会由于破浪的破碎而引起波浪爬高。它在单桩附近能够产生将近两倍的波高。

图  9  波面历时曲线

Figure 9.  Surface elevation record

图10为单桩式风机在不同工况下所受到的水平波浪力的历时曲线，截取的时段为400~420 s，并且将工况1-工况4数值模拟的结果与两种计算最大破波力的公式进行对比。表3为单桩式风机在不同工况所受到的最大水平波浪力的统计值。其中两种计算公式是通过在MATLAB中进行编程计算得到的结果。

图  10  不同工况的水平波浪力的对比

Figure 10.  Comparison of horizontal hydrodynamic load in different LCs

由图10和表3可知：

单桩式风机受到非破碎波（工况1、工况2）作用下的水平波浪力较小，在没有使用拉伸（工况1）和使用拉伸（工况2）时候，两者最大水平波浪力相近，工况2的最大水平波浪力约为1.560 MN，而工况1的最大水平波浪力约为1.420 MN，使用波浪拉伸会对非破碎波水平波浪力的计算提高了约10%。单桩式风机受到破碎波（工况3、工况4）的水平波浪力相比于非破碎波比较大，在没有使用拉伸，即工况3计算得到的最大水平波浪力为3.250 MN，使用拉伸时，即工况4计算得到的最大水平波浪力为8.170 MN，大约是工况3的2.5倍左右，即使用波浪拉伸会对破碎波水平波浪力的计算提高约150%。采用公式（6），即美国海岸防护手册建议修正的Morison公式可知破碎波波高为9.36 m，计算得到最大破碎波力为7.928 MN。采用公式（8）和公式（9），即冲击-绕流法计算得到最大破波力为9.408 MN。

在浅水海域和波浪易破碎区域单纯使用Morison公式计算破波力会较低，而使用波浪拉伸方式可以更精确地计算破碎波的波速，可以提高计算破波力的精度，从而可以更准确的描述破碎波。单桩式风机受到破碎波作用时，波浪由于破碎，会释放大量的能量，自由表面的水质点的速度较之非破碎波有很大的差异，波的运动也会表现出很强的非线性，使用拉伸计算方法可以将水平波浪力积分到瞬时表面，所以破碎波的计算使用拉伸会更加接近实际情况下单桩基础受到的破波力。通过非破碎波与破碎波的对比可看出，破碎波产生的波浪力远远大于非破碎波产生的水平波浪力，可能会对结构造成很大的损坏，故单桩式海上风机的设计需要考虑破碎波的影响。

本节针对研究在不同波浪和风载荷组合情况下单桩式风机塔筒顶部结构的动力响应。其中风机结构的整体坐标系如图11所示。

图  11  整体坐标系

Figure 11.  Global coordinate system

图12给出了单桩式风机塔筒顶部X方向的剪力和弯矩的历时曲线。模拟了四种不同的工况（工况1、工况3、工况5、工况6），对比了单独浪和风浪联合作用下单桩式风机的耦合响应。由图12可知：

图  12  不同工况下塔筒顶部X方向剪力和弯矩的时域图

Figure 12.  Time-series of tower top shear force and bending moment in the X direction under various load conditions

图  12  不同工况下塔筒顶部X方向剪力和弯矩的时域图

Figure 12.  Time-series of tower top shear force and bending moment in the X direction under various load conditions

1）单桩式风机塔筒顶部X方向剪力和弯矩在破碎波和定常风（工况6）作用下的时域运动响应最大。剪力的最大值为776 kN，最小值为733 kN，均值745.1 kN；弯矩的最大值为2.570 MNm，最小值为1.320 MNm，均值1.838 1 MNm。非破碎波和定常风（工况5）作用下塔筒顶部X方向剪力和弯矩的运动响应比工况6小，剪力和弯矩的最大值分别为755 kN和2.470 MNm，均比工况6的值要小。工况1和工况3分别是单独非破碎波和单独破碎波的工况，工况3下剪力和弯矩的动力响应大于工况1。破碎波对于塔筒顶部X方向的剪力和弯矩的运动响应影响要大于非破碎波。

2）定常风作用下单桩塔筒顶部X方向剪力和弯矩的动力响应远远大于单桩未受到风载荷的作用。由于塔筒顶部更接近风机的轮毂中心，此时风载荷远远大于波浪载荷在此处的贡献，所以，有无风载荷会对单桩塔筒顶部的剪力有较大影响。同样，由工况3和工况6的对比可得知，有风载荷和无风载荷的对比也可验证这一结论。

图13为单桩式风机在四种不同工况作用下塔筒顶部X方向剪力和弯矩的功率谱密度曲线。由图13-（a）可知：

图  13  不同工况的塔筒顶部X方向剪力和弯矩的功率谱密度曲线

Figure 13.  Power spectral density of tower top shear force and bending moment in the X direction under various load conditions

图  13  不同工况的塔筒顶部X方向剪力和弯矩的功率谱密度曲线

Figure 13.  Power spectral density of tower top shear force and bending moment in the X direction under various load conditions

1）单桩式风机在单独非破碎波（工况1）作用下运动响应较小，功率谱密度曲线呈现“单峰”状，最大的频率响应在0.28 Hz左右，而塔筒一阶固有频率的值为0.28 Hz，单独非破碎波可以激发塔筒的一阶固有频率。单独破碎波（工况3）的功率谱密度曲线呈现“三峰”状，频率为0.22 Hz、0.28 Hz和0.33 Hz，分别对应着波浪的2倍波频，塔筒一阶固有频率和波浪的3倍波频，可看出3倍波频的峰值最大，说明破碎波波浪主要的能量集中区域在3倍波频峰值附近。

2）非破碎波和定常风（工况5）作用下塔筒顶部

X方向的功率谱密度曲线只有一个峰值为0.6 Hz，对应NREL 5 MW单桩式风机的3P频率，从而体现出风载荷的对塔筒的频率响应要大于波浪载荷。破碎波和定常风（工况6）作用下结构的频率响应的功率谱密度曲线，该工况下的最大频率效应为波浪三倍频率0.33 Hz，其次是波浪的2倍波频0.22 Hz和风机的3P频率0.6 Hz。单独浪能激发塔筒的1阶固有频率，而风浪联合作用下，3P效应更加明显。同时破碎波更能激发波浪的2倍波频和3倍波频，由于波浪的破碎可能会改变波浪的波要素，所以更容易激发波浪的2倍波频和3倍波频。

图13-（b）为塔筒顶部X方向弯矩的功率谱密度曲线，由图可知：只有工况5和工况6有峰值，且在0.6 Hz处，也就是风机的3P频率，图中只存在弯矩的3P频率响应，说明波浪对弯矩的频率响应的贡献很小，弯矩的频率运动响应效应主要依赖于风载荷，在有风的情况下，能量都集中在3P频率附近。

泥面线处剪力和弯矩的极限值对于固定式风机，尤其是单桩式风机有重要影响，是评估风机系统稳定性的重要指标。泥面线的结构响应的极限值在设计范围内可保证风机正常工作。因此，针对单桩式风机在不同海况作用下的运动响应问题具有十分重要的研究意义。本节针对研究在不同波浪和风载荷组合情况下单桩式风机泥面处的动力响应。图14是单桩式风机泥面X方向的剪力和倾覆力矩的历时曲线。由于风和波浪都是沿着X方向作用在基础上，所以X方向的内力会承受最大的载荷。由图14可知：

图  14  不同工况的泥面X方向剪力和弯矩的时域曲线

Figure 14.  Time-series of shear force and bending moment in the X direction at the mudline under various load conditions

图  14  不同工况的泥面X方向剪力和弯矩的时域曲线

Figure 14.  Time-series of shear force and bending moment in the X direction at the mudline under various load conditions

1）在破碎波和定常风（工况6）作用下，泥面X方向的剪力的运动响应最大，最大响应在X的负方向上，达到-4.400 MN，耦合效应最明显，在单独破碎波（工况3）作用下，剪力的最大运动响应也在X的负向，最大剪力达到-3.670 MN，两者相差了730 kN，由于风载荷引起的差异，刚好也是额定风速下风机所受到的推力。此时，波浪载荷在泥面剪力处的贡献大于风载荷的，波浪载荷成为主要的控制载荷，所以，泥面剪力的运动响应波浪载荷明显大于风载荷。通过对比工况1和工况3可知，在受到单独非破碎波影响下，单桩式风机泥面处最大值为1.440 MN，最小值为-1.440 MN，幅值为2.880 MN，而受到单独破碎波作用下，泥面处最大值为2.870 MN，最小值为-3.670 MN，幅值达到6.540 MN，破碎波作用在泥面处的剪力的运动响应的幅值是破碎波的两倍多，说明破碎波对单桩基础泥面处剪力的影响要明显大于非破碎波。

2）破碎波和定常风（工况6）作用下，单桩式风机泥面X方向的倾覆力矩达到-104 MN·m，明显大于其他三种工况，非破碎波和定常风（工况5）作用下泥面最大倾覆力矩达到-84.900 MN·m，相差了将近10.100 MN·m，所以破碎波对于泥面处的倾覆力矩有着重要的影响，在结构设计中应不容忽视。在单独非破碎波作用下（工况1），泥面处倾覆力矩的最大值为12.000 MN·m，最小值为-9.190 MN·m，均值为1.476 1 MN·m，而在工况5作用下，泥面处倾覆力矩的最大值为-65.600 MN·m，最小值为-84.900 MN·m，均值-75.976 6 MN·m，所以风载荷对于泥面处倾覆力矩的影响很重要，是设计中不可忽视的载荷之一。

图15为单桩式风机在四种不同工况作用下泥面处X方向剪力和倾覆力矩的功率谱密度曲线。由图15可知：

图  15  不同工况的泥面X方向剪力和弯矩的时域图和功率谱密度图

Figure 15.  Time-series and power spectral density of shear force and bending moment in the X direction at the mudline under various load conditions

图  15  不同工况的泥面X方向剪力和弯矩的时域图和功率谱密度图

Figure 15.  Time-series and power spectral density of shear force and bending moment in the X direction at the mudline under various load conditions

单独非破碎波（工况1）作用下剪力和倾覆力矩的功率谱密度曲线只有单峰，为波频处，即0.33 Hz处有峰值，同样的其他三种工况也能观察到波频处有峰值。但是在0.22 Hz处，即波浪的2倍波频观察到破碎波（工况3、工况5）在此处有峰值，非破碎波并不能激发波浪的2倍波频，同样的3倍波频处，即0.33 Hz也只能观察到破碎波的峰值。所以破碎波作用在单桩基础上，剪力和弯矩的频率响呈现出高倍波频能够被破碎波激发，而非破碎波不容易激发波浪的倍频。

本文采用FAST和MATLAB软件，主要针对单桩式海上风机在我国东海浅水海域的海况条件下，在非破碎波和破碎波工况下（有风和无风），同时考虑波浪的拉伸的方法应用到破波力的计算，对单桩式海上风机的动力响应进行了对比分析，可以得出以下结论：

1）研究结果表明：不同的模拟水质点运动学的拉伸方法会导致不同的响应值，基于惠勒和垂直拉伸方法的极端响应的概率分布可能存在显著差异，导致设计中应采用何种方法存在不确定性。因此，需要进行进一步的研究来解决这个问题。

2）使用垂直拉伸方式来模拟水质点的运动学通过FAST和MATLAB进行比较，比较结果表明在不同节点处两者显现出一致性，说明FAST的数值模型可以得到了验证。使用波浪拉伸方式对波浪载荷的计算有很大的影响，对于非破碎波，使用波浪拉伸可以使得波浪载荷计算提高约10%左右，而对于破碎波，使用波浪拉伸可以使波浪载荷计算提高近150%，可以大大提高计算破碎波波浪力的精度。但是考虑波浪拉伸的影响只分析了波浪力的差异，在耦合分析计算中未考虑波浪的拉伸影响，在以后的工作中会考虑波浪拉伸对结构的动力响应。

3）通过在塔筒顶部的内力反应可以知道，风浪联合作用下风机结构动力响应更加显著，荷载作用在塔筒上产生的内力，风载荷处于主导地位，相比于风浪共同作用，单独波浪更容易激发塔筒的一阶固有频率。从泥面处的内力反应可以得知，破碎波和定常风作用下泥面处剪力和倾覆力矩的动力响应最大，破碎波作用下，结构泥面处内力的频率响应的高倍波频更容易被激发，而非破碎波不容易激发波浪的倍频。破碎波对海上风机的设计有着重要的影响，在结构设计中应不容忽视。

本文的研究成果对我国浅海海域尤其是波浪破碎区域的单桩式海上风电的设计和开发提供了一定的参考。此外，破碎波对于单桩式海上风机有着重要的影响，在浅水海域海上风机的设计和研究应该考虑破碎波载荷。