Research on Pile-Soil Interaction of 10 MW Large Monopile Offshore Wind Turbine

石化能源为主的能源结构是导致当今世界上雾霾等环境问题的主要原因。中国作为世界最大的能源消费国及二氧化碳排放国，亟须探索以“双碳”目标为导向的能源安全和结构优化的新路径^[1]。风力发电作为一种新型发电技术，得到了广泛运用，尤其是近些年随着新技术与新设备的运用，风力发电得到很大发展^[2]。除了在陆地上布置风机外，海上风机由于近海地区风力条件更强、更稳定而发展迅速。根据目前的数据，欧洲海岸已经建造大约1 000多个海上风机。按照目前的发展趋势，到2030年，全球海上风电总容量将达到66.5 GW^[3]。开发和利用海上风能资源，对于治理大气雾霾、调整能源结构具有重要意义，是实现经济低碳、绿色、可持续增长的有效途径。

单桩基础凭借其加工制造简单、安装便捷、结构受力明确等优势，被广泛应用于近海风力发电工程项目中。通过分析大量海洋风力发电工程项目资料，国际上采用单桩基础的海上风电场占总量的53%左右，已建成海上风机基础75%为单桩基础。海上风机单桩基础均属于大直径桩，其直径一般大于5 m，随着海上风电大型化发展，风电单机装电容量的增加，其直径甚至可达到8 m以上。近年来的相关研究表明对于海上风机单桩基础，其大直径效应不可忽略^[4]。

Matlock^[5]最早提出了软黏土的p-y曲线、Reese^[6]、Mcvay^[7]分别提出了硬粘土和砂土的p-y曲线，3种经典p-y曲线已经被列入美国石油协会规范。Brown等^[8]利用三维有限元模型分析了水平荷载作用下桩-土系统的响应规律。Kouda等^[9-11]通过离心机试验研究了水平荷载作用下桩基础的承载特性。Damgaard等^[12]使用半解析频域解决方案来评估桩土系统的动态阻抗函数，以校准集总质量模型集成到气动弹性多体代码中。他们为单桩基础使用了3种不同的模型：表观固定性模型、海床水平的固定支撑模型和耦合参数模型。该研究的作者得出结论，就海床层面的疲劳而言，桩土相互作用现象对于设计很重要。挪威船级社在其海上风机基础设计规范(DNVGL-ST-0126)中指出，传统p-y曲线仅用于细长的桩，其直径不超过1 m，用于大直径单桩的p-y曲线需要进一步验证^[13]。目前，海上风机单桩基础桩土作用主要沿用传统设计方法，如AF法，p-y曲线法。近年来，挪威岩土工程学会（NGI）开发了一种新的单元模型，宏单元法（Macroelement^[14]），较好地解释了海上风机土壤与结构的相互作用。

综上所述，本文针对10 MW大型单桩风机，采用多种不同的桩土作用模型，考虑波浪、风和风机系统的相互作用，揭示10 MW大型海上单桩风机在风浪联合作用下的动力学耦合机理。研究发现宏单元法考虑了非线性的刚度与塑性，使得宏单元法在对比其他传统桩土模型时有很大的优势。

在本文中4种桩土模型将用来参与计算，这4种模型是AF法、p-y曲线法、刚性连接法与宏单元法。AF法假设桩延伸到特定的深度，在海床上产生的侧向位移和旋转与真正的桩埋在土壤中相同^[15-16]。AF法使用的模型用一个虚构的桩取代地基和土壤，桩固定在泥线以下的底部。公式（1）给出了在一个平面内弯曲的轴对称悬臂梁的二维刚度矩阵，忽略了竖向和扭转自由度。刚度矩阵基于欧拉-伯努利梁理论。

式中：

L_AF ——桩的长度（m）；

EI_AF ——桩截面的抗弯刚度（N·m^-1）。

刚度矩阵包括水平刚度K_HH、旋转刚度K_RR以及水平自由度与旋转自由度之间的交叉耦合刚度K_HR=K_RH，其中对角线强度系数(交叉耦合系数)相互匹配，这是单桩式海上风电中的常见做法。但是，必须认识到，如果不匹配非对角线系数，海底的平移和旋转位移就不能精确地再现。除此之外，Løken等^[17]开发了一种新的方法来改进传统的AF模型。所提出的方法可以考虑刚度矩阵的所有项。

Winckler^[18]在研究海上油气相关结构时提出了一种方法，即用弹性梁取代桩，用一组独立弹簧取代土。这些弹簧分布在单桩嵌入的底部。经过了这段时期的发展，它现在被称为p-y曲线法。实际使用的p-y曲线为API^[19]推荐的曲线，API推荐的砂土中桩的p-y曲线如下:

式中：

p_u ——泥面线以下z深度处的极限水平阻力（N·m⁻¹）；

$k$ ——最初的地基反力系数随深度的增加速率（N·m⁻³）；

A ——静荷载或循环外荷载条件下的经验系数；

Y ——z深度处的横向位移（m）。

本文参考NGI^[20]的相关规定，采用一系列p-y曲线来表征下部单桩在横向荷载下的相对位移。沿着单桩向下方向，每隔0.75 m有61条p-y曲线。每条p-y曲线由22对点来定义。重要的是，每单位长度桩的横向土阻力是与单位长度桩的横向位移相对应的。桩的横向土阻力p与桩横向位移y之间存在对应关系，桩横向位移y随桩长呈非线性变化。

Bazeos^[21]和Lavassas等人^[22]采用刚性连接法，在他们的研究中，上部单桩风机模型是由塔筒的梁单元和上部风机的集中质量建模的。在这个模型中，完全忽略基础和周围的土壤，桩与海底刚性连接。

宏单元法出自REDWIN项目^[14]中，它是在多面塑性框架内制定的具有随动硬化的弹塑性模型^[23]。它将桩和周围土壤简化为一组线性和非线性荷载-位移关系，在海床面的6个自由度中，将桩和结构的其余部分分开。它可以表示在实验测试和现场观察到的非线性滞回荷载-位移响应，包括水平荷载和弯矩之间的耦合响应。该模型已被证明与土壤和单桩的有限元分析结果、大型桩试验结果以及安装在北海的海上风机的全尺寸现场测量结果一致^[24]。

为了更好地研究土壤与结构的相互作用，本文对上部风机进行了简化。仅对支撑结构进行建模，并使用集中质量和惯性矩表示转子机舱组件（Rotor-Nacelle Assembly）。为了更好地建立模型和相关结果的输出，坐标系的建立如图1所示。转子机舱组件采用集中质量点进行建模，质量点的质量和转动惯量如表1所示，其采用的是国际能源署（IEA）开发设计的10 MW风机。

Figure 1.  Selection of coordinate system

本文中使用的塔筒采用的是海上DTU 10 MW风机中的塔筒。图2显示了塔筒和单桩的示意图和高度信息。对于塔筒和单桩，两者的材料特性相同，其中杨氏模量（E）为210 GPa，剪切模量（G）为80.8 GPa，密度（ρ）为8 500 (kg·m⁻³)。塔顶直径为5.5 m，壁厚为30 mm，塔基直径为8.3 m，壁厚70 mm。塔顶部到塔基的距离为105.63 m。单桩从塔基底部到海床的总长度为40 m，壁厚为110 mm。结构部分水深为30 m。对于AF模型，建模期间需要在海床底部延伸一段。延伸部分的外径与单桩的外径相同。其刚度矩阵基于欧拉-伯努利梁理论（见公式1）。表2总结了AF方法的参数和海床等效刚度。

Figure 2.  Schematic diagram of tower and monopile

本文使用了SACS与OpenFAST两种软件，其中AF法，p-y曲线法刚性连接法（Rigid）通过在SACS软件中进行建模，宏单元法（REDWIN）在OpenFAST中进行建模。SACS是一个集成的有限元结构分析套件，专门为海上结构(包括石油平台和海上风力)的设计、制造、安装、操作和维护提供服务。在这项工作中，利用SACS中的Extract Mode Shapes模块，可以获得每个桩-土模型的特征频率和模态。利用Module Time History Wave, Wind and Force模块，得到了各种桩-土模型在风和波条件下的响应。此模块采用了模态叠加方法，通过在软件中对结构施加不同荷载，可以得到结构在关键节点位置的响应。OpenFAST也是由各个模块组成，比如有Aerodyn、Subdyn和InflowWind、Beamdyn等等模块。通过对各个模块进行使用，即可得出相应的结果。图3显示了4种桩土模型的示意图。为了重点验证桩土模型与上部单桩风机的相互作用行为，本文对上部风力机进行了简化，仅对支撑结构进行建模，并使用集中质量点和转动惯量来表示转子机舱组件(RNA)。在模型中，各节点信息及单桩位置如图4、表3所示。

Figure 3.  Schematic diagram of pile-soil model

Figure 4.  Schematic diagram of node information

本文重点验证桩土模型与上部单桩风机的相互作用行为，表4列出了本文中模拟的工况，其中包括静态分析、模态分析、单独风、单独浪和风浪联合条件下的模拟。需要注意的是，涉及风的荷载情况(即单独风和风浪联合)的响应仅用于验证目的，且在使用的模型中缺乏气动阻尼，故不能被认为代表实际运行条件下的风机。在第3节中，将对表4所示的研究内容进行分析，来对比各个桩土模型的特点。

LC1是静态分析，其各种桩土模型是仅以重力加速度作为唯一的外部荷载。通过这种方式来保证结构模型的正确性，并且对特定位置处的偏转和位移进行比较，也可以得到各种模型之间的关系。图5显示了在仅限重力条件下沿x轴（前后方向）的塔顶水平位移，这种位移是由于单桩式风机偏心结构引起的，这受到结构整体刚度的影响。从图5中可以看到刚性连接法在重力条件下的水平位移是最小的，这是由于刚性连接法在海床位置处是完全固定的，这对于整个塔筒-单桩模型在运动时有非常大的限制。其次AF法的水平位移排在第二低，这是由于AF法属于线性方法，这使得其塔顶水平位移比非线性的p-y曲线法得到的水平位移略低。在图中我们看到宏单元法（REDWIN）得出的结果是4个模型中最大的，这是因为此模型考虑了非线性的刚度与塑性，这使得宏单元法得出的塔顶水平位移最大。

Figure 5.  Horizontal displacement of tower top

模态分析是通过对结构系统的特征频率，阵型等来对桩土模型进行一个对比。表5显示了4种模型得到的在前后方向（Fore-Aft）与两侧方向（Side-Side）的基频。从表中可以得到，刚性连接法的特征频率与其他3种模型相比具有一定的差距，这与海床位置处的限制有很大关系。同时可以看到AF法，p-y曲线法和宏单元法得出的一阶频率差距不大，有着很好的一致性，宏单元法得到的一阶频率最高。而对于二阶频率来说，其得到的结果不如一阶频率一样有着较好的相似性，这与各种模型在海床位置处的限制方式有很大关系。

图6显示了在前后方向和两侧方向上4种桩土模型的阵型对比，其位移进行了归一化处理，有利于更加直接的来进行模型之间的对比。从图中可以看到，一阶阵型中最大的位移在塔顶位置处，4种桩土模型表现出了一致性。对于二阶阵型，最大位移则是出现在约整个结构高度的三分之二高度处，在这种情况下，刚性连接法得出的结果在下部产生的偏转位移较小，而在塔顶位置处AF法与p-y曲线法有着最大的偏转位移。

Figure 6.  Comparison of four model formations

在分析了4种桩土模型的静态荷载情况和特征频率后，本节将对单独风情况下的桩土模型结果进行分析。如第1节所述，上部风机被简化成了一个集中质量点，故在本节和接下来的风浪联合工况下将作用在塔顶处的风推力作为外力对整个模型进行一个加载。图7塔顶在水平方向（沿x方向）的加速度时域图，从图中可以看到，4种桩土模型具有一定的一致性，同时刚性连接法与宏单元法在幅值上比其他两种模型稍大。图8是塔顶在水平方向的加速度PSD图，其主要的激励频率和特征频率在图中有所显示（红色垂线中间部分为一阶FA频率范围，蓝色垂线中间部分为二阶FA频率范围，虚线为3倍转子频率3P）。从图中可以看到，除了结构自身存在的特征频率之外，3P成为了对塔顶加速度影响的主要激励因素。图9对比了不同桩土模型塔顶水平加速度功率谱密度在特定频域区间内的总和，一个是0.2 Hz到0.3 Hz这个频率区间，另一个是1.0 Hz到1.5 Hz这个频率区间，考虑到在PSD图中4个桩土模型在两个范围内具有不小的响应，故可以认为是一阶FA频率和二阶FA频率的激励效果。从图中可以看到，刚性连接法虽然在一阶FA频率附近的总和与其他三种模型较为接近，但是在二阶FA频率附近的总和却较为低。宏单元法比AF法与p-y曲线法在一阶FA频率附近的总和大，这是因为宏单元法考虑了非线性的刚度与塑性。从图中可以看到一阶FA频率中的幅值要比二阶FA频率中的要大的多，这在整个运动中占主要作用。

Figure 7.  Time domain diagram of horizontal acceleration of the tower top

Figure 8.  PSD diagram of horizontal acceleration of the tower top

Figure 9.  The sum of the acceleration power spectral density of the tower top in a specific frequency domain

图10显示了海床倾覆力矩（沿FA方向）时域图，图11显示了海床倾覆力矩PSD图，其得到的结果与塔顶加速度的情况类似。其主要的区别是在0 Hz时PSD中显示的值不为0，这是4种桩土模型在受到风荷载时弯矩的平均值，这主要是作用在结构上的风推力导致的。图12同样对特征频率范围内的功率谱密度进行统计，得到的结果仍然是宏单元法的值最大，这与前面的结论相似。

Figure 10.  Time domain diagram of sea-bed overturning moment

Figure 11.  PSD diagram of sea-bed overturning moment

Figure 12.  The sum of the power spectral density of sea-bed overturning moment in a specific frequency domain

LC4不规则波浪的模拟，其有义波高为1.25 m，波浪周期为5.5 s。在模拟过程中，拖曳力系数${C_D}$,${C_M}$分别取为1和2。由于波浪的峰值频率（$1/{T_{\rm{p}}}$）落在一阶FA频率范围内，故在此波浪条件下，我们只考虑0到0.5 Hz之间内的响应。图13、图14显示了塔顶在水平方向（沿x方向）和海床倾覆力矩（沿FA方向）的加速度PSD图。从图中可以看到，刚性连接法在一阶FA频率范围内中有着最小的响应，这可能是与其他3种桩土模型相比，其自身的频率较高。图15同样对比了不同桩土模型塔顶水平加速度，海床倾覆弯矩功率谱密度在特定频域区间内的总和（0.2~0.3 Hz）。从图中可以看到刚性连接法的总和值是最小的，远远小于AF法与p-y曲线法，同时宏单元法得到的结果是最大，这个可能是到宏单元法在海床位置荷载限制方式有关。

Figure 13.  PSD diagram of horizontal acceleration of the tower top

Figure 14.  PSD diagram of sea-bed overturning moment

Figure 15.  The sum of the power spectral density in a specific frequency domain

风浪联合模拟结合了LC3工况中的单独风条件和LC4工况中的单独浪条件。图16显示了塔顶加速度沿x轴的PSD图，图17是相应的PSD之和，它可以与前面单独风条件下的情况进行比较。对比可以发现，两者之间的主要差别是，风浪联合作用下一阶FA频率处的响应会增大，这是由于波浪叠加作用产生的影响。同时宏单元法二阶FA频率会有一些轻微程度的偏移，这可能是由于波浪作用条件下水面变化引起的。图18显示了海床倾覆力矩PSD图，图19是相关的PSD总和。和前面单独风，单独浪得到的结果类似，宏单元法在一阶FA频率和二阶FA频率附近的响应最大。

Figure 16.  PSD diagram of horizontal acceleration of the tower top

Figure 17.  Sum of the acceleration power spectral density of the tower top in a specific frequency domain

Figure 18.  PSD diagram of sea-bed overturning moment

Figure 19.  The sum of the power spectral density of sea-bed overturning moment in a specific frequency domain

在本研究中，4种不同的桩土模型用于10 MW大型单桩风机的比较研究。本文采用SACS与OpenFAST软件进行建模，并对上部风机进行了简化。通过在单独风、单独浪和风浪联合条件下输出典型位置的信息，得出以下结论：

1）刚性连接法在重力条件下的水平位移是最小的，这是由于刚性连接法在海床位置处是完全固定的。AF法的水平位移排在第二低，这是由于AF法属于线性方法，这使得其塔顶水平位移比非线性的p-y曲线法得到的水平位移略低。宏单元法（REDWIN）得出的结果是4个模型中最大的，这是因为此模型考虑了非线性的刚度与塑性，这使得宏单元法得出的塔顶水平位移最大。

2）刚性连接法的特征频率与其他3种模型相比具有一定的差距，这与海床位置处的限制有很大关系。同时可以看到AF法，p-y曲线法和宏单元法得出的一阶频率差距不大，有着很好的一致性。而对于二阶频率来说，其得到的结果不如一阶频率一样有着较好的相似性，这与各种模型在海床位置处的限制方式有很大关系。

3）宏单元法比AF法与p-y曲线法在特征频率附近的功率谱密度总和大，这是因为宏单元法考虑了非线性的刚度与塑性，这表明了宏单元法在对比其他传统桩土模型时有很大的优势。同时一阶FA频率中的幅值要比二阶FA频率中的要大的多，这在整个运动中占主要作用。