大直径单桩海上风机横向自振频率实用计算方法

黄永健，余云燕，付艳艳

（兰州交通大学土木工程学院，甘肃兰州 730070）

风能作为一种优质、清洁的可再生能源，目前正被欧盟国家大范围推广使用。截至2020年，欧洲已在12个国家累计建造了116 个海上风电场，共计5 402 台海上风力发电机，总装机容量达25.01 GW［1］。与欧洲国家相比，中国海上风电技术起步较晚，在海上风电功率与装机容量上仍存在较大差距，目前正处于高速发展阶段。

在海上风机的设计与运营期间，风电系统会受到涡轮机和叶轮扫掠过程中的谐振作用，因而需避免风电系统一阶横向自振频率与涡轮机的转动频率（1P 频率）和叶轮扫掠频率（2 P/3 P 频率）发生重叠而引起共振效应。而海上风电系统设计一般采取“软-刚”模式［2］，要求风电系统的一阶横向自振频率介于1 P 频率与2 P/3 P 频率之间。为安全起见，作为行业标准的DNV 规范［3］建议，在风机系统设计时，应在1 P 频带与2 P/3 P 频带之间留有10%的安全余度，使得频率设计容许区间受限，风机系统频率计算精度要求进一步提高。

对大直径海上风机而言，塔筒变截面特性的考虑方式以及桩、土相互作用模型的选择是风机系统频率准确预估的关键。对于塔筒变截面的特性，常见的研究方法是通过等效截面的方式［4-7］加以考虑，或者辅以有限元软件对塔筒进行仿真模拟。采用等效截面的方式在一定程度上反映塔筒的变截面特性，但并不全面，相比之下，有限元建模的方式模拟效果较好，但在设计前期由于风机结构尺寸要素的频繁改动，常需重复建模，工作量较大。同时，为准确计算单桩海上风机系统横向自振频率，研究人员采取了不同的计算模型对桩、土相互作用加以考虑，包括等效桩长法［8］，独立弹簧法［5，7，9］、耦合弹簧法［4，10］、Winkler 地基模型［6］以及p-y曲线法［11-12］等。用弹簧替代复杂的桩、土相互作用能较大程度简化计算，但弹簧刚度的选择依赖于对土体刚度的准确预估，而对于成层土，地基刚度的取值尚未形成系统的认识。而p-y曲线法是规范［3］推荐的分析方法，适用于各种土体的桩、土相互作用模拟，常配合有限元复杂建模使用。然而，在目前风机系统前期设计过程中，由于桩身尺寸及土体参数的不确定性，常需多次建模进行调试，过程较为繁琐。因此，海上风机的前期设计中，风机系统横向自振频率的计算仍需一种既能全面考虑塔筒变截面特性，又能考虑桩、土相互作用的快速计算方法。

文中基于Timoshenko 梁理论，以Winkler 地基作为桩、土相互作用模型，采用多段均匀梁组合的方式建立振动方程，旨在为大直径单桩海上风力发电机的前期设计提供一种思路简单、计算精确、利于提升工作效率的频率预估方法。该方法能精准模拟风机塔筒的变截面特性，同时兼顾考虑土体分层的特点，适用于求解类似结构的各阶横向自振频率。

1.1 计算模型的建立与简化

如图1（a）所示，典型的大直径单桩海上风机支撑结构主要由叶轮—机舱组合件（Rotor Nacelle Assembly，即RNA）、塔筒、连接段与大直径单桩4部分组成。风机顶部为风机叶片、轮毂和涡轮机组成的叶轮-机舱组合件。RNA 与风机平台之间为变截面塔筒，塔筒截面直径在平台处最大，随高度的增加逐渐减小。泥面以下为大直径单桩，塔筒与单桩之间采用连接段和刚性法兰连接。

以桩底为原点，沿高度增加方向为x轴，考虑剪切变形和转动惯量的影响建立Timoshenko梁横向振动方程进行求解，风机整体计算模型如图1（b）所示。风机整体计算模型由塔筒、连接段以及全埋入单桩3部分组成，其基本假设如下：

（1）风机顶部采取集中质量mRNA的方式对叶轮-机舱组合件进行简化；

（2）将塔筒分为NT段，每一段选用所截取圆台塔筒的中位线所在截面作为该段等效截面，考虑塔筒截面直径沿高度方向线性变化，高度为LT，壁厚恒定为tT，底部和顶部外径分别用Db和Dt表示；

（3）假定连接段为等截面段，截面直径与桩身直径相同，连接段高LS，外径为DS，壁厚为tS；

（4）风机泥面以下考虑桩、土相互作用的影响，以Winkler 地基模型进行模拟，桩身按土体分层情况分成NP段，单桩埋深LP，外径为DP，壁厚为tP。

（5）风机支撑结构总共分为(NP+NT+1)段，以r表示风机支撑结构中的第r段，r=1,2,···,(NP+NT+1)，第r段梁长为lr，风机分段如图1（c）所示。

图1 风机系统计算模型Fig.1 Analysis models of OWTs

1.2 桩身部分求解

文献［13］指出，土体阻尼对风机系统一阶横向自振频率的影响较小，可忽略不计。泥面以下的大直径单桩，假设桩径沿深度方向不变，采用Winkler地基模拟桩、土相互作用，并采用分段的方式考虑分层土体的不同力学特性，建立桩身部分的横向振动方程［14］为：

式中，i表示所在的土层序号，i=1,2,3···,NP，vpi(x,t)和θpi(x,t)分别为泥面以下单桩的横向位移和截面转角；
Ap和Ip分别为泥面以下单桩的截面面积和截面惯性矩，下标p表示单桩的相关变量，弹性模量E、剪切模量G与密度ρ为材料常数，取决于材料属性，κ为截面剪切系数，可通过κ=2(1+υ)/(4+3υ)取值［15］，其中v为泊松比，钢材取0.3。ki为桩周土弹簧刚度，与土体的弹性模量ESi有关，可根据Markris和Gazetas［16］的方法ki=1.2ESi取值。

令vpi(x,t)=ϕpi(x)eiωt，θpi(x,t)=φpi(x)eiωt，变量分离后可得特征方程：

式中：ω为风机系统横向自振频率；
Fpi=s1p(s1ps2ps3p-s2pηQi-1)+ηMi；
。

式（2）的通解按照参数的不同存在以下4 种情况：（1）Hpi＞＞Fpi＞0；
（2）Fpi＜0；
（3）Hpi＜0，Hpi2＞Fpi＞0；
（4）Fpi＞Hpi2。算例试算结果表明，大直径海上风机系统的一阶横向自振频率属于第4 种情况，可解得相应的振型函数为：

1.3 塔筒及连接段部分求解

风机运行过程中，塔筒和叶轮—机舱组合件会受到气动阻尼的影响，连接段也会受到水动力阻尼的影响。但文献［17］指出，气动阻尼和水动力阻尼对风机系统一阶横向自振频率的影响并不显著，由此建立塔筒及连接段部分等效截面的无阻尼横向振动方程：

式中：j表示第j段，j=(NP+1),(NP+2),(NP+3),···,(NP+NT+1)，vj(x,t)、θj(x,t)、Aj和Ij分别为第j段梁的横向位移、截面转角、等效截面面积和等效截面惯性矩。

令vj(x,t)=ϕj(x)eiωt，θj(x,t)=φj(x)eiωt，根据文献［18］可得梁横向位移和截面转角分别为：

式中，C1j,C2j,C3j,C4j为待定系数，其中，

1.4 横向频率特征方程

如图1（c）所示，由第r段梁与第r+1段梁在连接点xr处的位移、转角、弯矩和剪力的连续关系［19］可得：

将式（3）和式（6）分别代入式（8），整理可得：

式中：Cr,C(r+1)分别为第r段梁与第r+1段梁振型函数的待定系数；
Zr为待定系数的系数矩阵。

考虑式（9）中r分别为1,2,···,(NP+NT+1)，作进一步推导可得和C1关系为：

式中，Z=。

如图1（b）所示，风机桩底边界条件按铰接考虑，即位移和弯矩为零：

塔筒顶部自由，集中质量的惯性力与塔筒内力平衡，有：

将式（3）代入式（11）中，整理可得：

将式（6）代入式（12），整理可得塔筒顶部边界条件：

式中，

联立式（10）和式（14），可得塔筒顶部边界条件表示如下：

进一步联立（13）和式（16），求解可得：

式中，R=[BTZ]T。

为使式（17）有非零解，其系数矩阵行列式必为零，即：

式（18）为包含系统横向自振频率ω的超越方程，基于matlab软件，通过二分法求解即可得到系统的横向自振频率。

如图1（b）所示，采用Winkler 模型考虑桩、土相互作用，对风机系统的一阶横向自振频率进行求解。选取文献［10］中的4 个海上风机作为算例，各风机具体规格参数如表1 所示。各风机地层情况及土层预估模量见表2。土层信息可通过表中数据来源作进一步了解。根据各风机地层信息，配合文献［20-21］中的经验数据对地层土体弹性模量进行预估，数值根据砂土的内摩擦角ϕ或黏土的不排水抗剪强度ccu线性插值确定。采用Arany的三弹簧法［10］所得结果及风机现场实测频率［10］对本方法的有效性进行评估，结果如表3所示，表中相对误差均由计算结果与现场实测频率比较得到。为观察塔筒变截面特性对风机系统一阶横向自振频率的影响，文中通过将塔筒分为1段，4段和16段的情况对变截面影响进行量化分析。

表1 风机规格参数表Table 1 Specifications of four offshore wind turbines

表2 地层属性及土体参数表Table 2 Ground profile and soil parameters

表3 整体计算模型一阶横向自振频率与实测频率对照表Table 3 Comparison on fundamental transverse frequencies of four turbines among results of proposed model and measured frequencies Hz

计算结果显示，在考虑桩、土相互作用的情况下，忽略塔筒变截面特性的影响，即当NT=1 时，计算结果与现场实测频率相比普遍偏小，误差范围在6.88 %～13.86 %而当NT=4 时，计算结果与实测频率相近，最大误差为4.07 %，计算精度基本满足要求。随着分段数目的增大，当塔筒按16 分段考虑时，本方法计算结果与Arany 的方法所得结果差异较小，与现场实测频率相对误差均在3.50%以内，对于风机系统前期设计而言是可以接受的。且从算例的计算结果可以看出，相比NT=4 时所得结果，NT=16 时所得结果在计算精度上更加稳定，在设计时适宜采用塔筒16 分段的方式进行计算。综上分析，在考虑了塔筒变截面特性及桩、土相互作用后，文中结果与现场实测频率基本一致，对风机系统一阶横向自振频率的预测是准确有效的。

风机系统频率控制要求严苛，而风机本身的规格参数、土体模量以及运营期间的海床平面变化等因素都会对其造成一定的影响。且在风机运营过程中，一阶横向自振频率起主要控制作用，其余各阶对风机的振动也存在一定的贡献，为进一步量化各因素对风机系统各阶横向自振频率的影响，现以算例中的Burbo Bank风电场风机为例，就各因素对风机前三阶横向自振频率的影响进行探讨。

3.1 塔筒、连接段高度对风机系统横向自振频率的影响

考虑泥面以上部位塔筒和连接段高度对风机系统横向自振频率的影响，保持风机其余尺寸要素不变，以塔筒高66 m与连接段高23 m为基准展开变量分析，计算结果显示风机系统横向自振频率随塔筒和连接段高度变化较为均匀。如图2 所示，塔筒高度分别变化-10 m，-5 m，5 m，10 m，即塔筒高度分别为56 m，61 m，71 m，76 m 时，风机系统一阶横向自振频率较原一阶自振频率分别变化了20.0 %，9.3 %，-8.2 %，-15.3 %；
二阶横向自振频率分别变化17.0%，8.0%，-7.2%，-13.8%；
三阶横向自振频率分别变化19.9%，9.2%，-8.0%，-15.0%。而当连接段高度分别变化-10 m，-5 m，5 m，10 m，即连接段高度为13 m，18 m，28 m，33 m 时，风机系统一阶横向自振频率分别变化10.8%，5.3%，-5.1%，-10.0%；
二阶横向自振频率分别变化23.6%，11.0%，-9.5%，-17.7%；
三阶横向自振频率分别变化17.9%，8.1%，-6.9%，-13.0%。分析显示，在塔筒和连接段既有高度变化量相同情况下，风机系统一阶横向自振频率对塔筒高度变化敏感性较强，而对连接段高度变化敏感性稍弱，风机系统二阶和三阶横向自振频率对连接段高度变化敏感性较一阶有所增强，在风机系统设计时应加以综合考虑。

图2 塔筒与连接段高度对风机系统前三阶横向自振频率的影响Fig.2 The first three order frequency affected by height of tower and transition piece

3.2 土体模量对风机系统横向自振频率的影响

分析土体模量变化对风机系统横向自振频率的影响，土体按均质土考虑，结果如图3所示。可以发现，当土体为黏土、粉质黏土、粉土、粉砂、细砂等土质时，地基模量在2 MPa到40 MPa之间［20］，此时随着地基土模量的增加，风机系统一阶横向自振频率从0.199 Hz迅速增加到0.309 Hz。而当土体模量大于40 MPa，即土体为密砂、岩性材料等土质［17］时，风机系统一阶横向自振频率逐渐趋于稳定，随地基土模量的增加仅有轻微的上涨。类似的趋势在土体模量变化对风机系统二阶和三阶横向自振频率的影响曲线中也有体现。且对比图3可以发现，土体模量变化对风机系统第三阶横向自振频率影响最大，对第二阶影响中等，对第一阶影响最小。总体而言，土体模量对前三阶风机系统横向自振频率影响极大，系统设计时必须考虑桩、土相互作用对风机系统横向自振频率偏移的影响。为减小风机运营过程中系统一阶横向自振频率的偏移，应在风电场选址时尽可能选择高模量土层地区进行风电场建设，或对软弱地基采取适当的加固措施以提高地基的土体模量。

图3 土体模量变化对风机系统前三阶横向自振频率的影响Fig.3 The first three order frequency affected by soil modulus

3.3 海床平面高度变化对风机系统横向自振频率的影响

海洋中由于海流的冲刷、搬运、堆积作用以及海底滑坡、地震等自然灾害的影响，海床平面高度时有变化。而海床高度变化对桩埋入深度和连接段高度均有直接影响，对海床高度变化进行定量分析，结果如图4。其中，当海床平面高度变化为正值时，意味着海床平面上升，桩埋入深度增加，而连接段高度减小。反之则为海床平面下降，桩埋深减小，连接段高度增加。从图4 可以看出，随着海床平面从降低6 m 到升高6 m，风机系统前三阶自振频率均出现类似线性的变化。其中，一阶横向自振频率从0.268 Hz增大至0.320 Hz，涨幅达18.4%，二阶横向自振频率涨幅达38.6%，三阶横向自振频率涨幅达21.2%。分析表明，风机系统前三阶横向自振频率受海床平面高度变化影响较大，在海洋滑坡等自然灾害频发地区，在设计时必须充分考虑海床平面高度变化对风机系统横向自振频率的影响。

图4 海床平面高度变化对风机系统前三阶横向自振频率的影响Fig.4 The first three order frequency affected by seabed height

3.4 大直径单桩的桩长与桩径对风机系统横向自振频率的影响

在风机系统设计过程中，桩的直径与埋入深度是影响风机系统横向自振频率的关键要素之一。为分析桩长与桩径的影响，现将土体考虑为单一土层，考虑风电场土体多为密砂层，故假定地基土模量为40 MPa，其结果如图5 所示，图中给出桩径分别为2.5 m、3.5 m、4.5 m、5.5 m、6.5 m 和7.5 m 时风机前三阶横向自振频率随桩长的变化曲线。变化曲线由上升段和平直线段两部分组成。在桩长较短时，系统横向自振频率随桩长的增加呈现迅速增加的趋势，而后随着桩长的增加，系统横向自振频率受桩长变化的影响逐渐减弱，而当桩长超过某一临界点后，桩长的继续增加对系统横向自振频率几乎不产生影响，曲线趋于平直线。同样的频率变化趋势在桩径上也有所体现。当桩径较小时，桩径变化对系统横向自振频率的影响较大，但随着桩径的继续增加，其影响程度在逐渐降低。风机系统前两阶横向自振频率在桩长变化下的影响曲线较为类似，而对于第三阶横向自振频率，当桩长较短时，拥有更大直径单桩的风机系统频率更小，但不同直径系统的频率差距不大，且桩径越大，其曲线临界点对应的桩长越长，在桩长足够长的情况下，对应的横向自振频率也更高。曲线表明，大直径长桩具有较好的抗频率偏移能力，在设计中应在保证经济合理性的情况下尽量予以采用。

图5 桩长与桩径的变化对风机系统前三阶横向自振频率的影响Fig.5 The first three order frequency affected by diameter and length of pile

基于Timoshenko 梁理论，考虑桩、土相互作用和风机的塔筒变截面特性，建立了大直径单桩海上风机横向自振频率的实用求解方法，并通过4 个实际风机的现场实测频率对文中简化方法的有效性进行验证。该方法在完成首次推导及matlab 程序录入后，仅需在程序中改动相应的风机尺寸以及土层参数，即可获得相应的风机系统横向自振频率，因而有利于大幅减少设计前期有限元重复建模所需要的时间，从而达到提升工作效率的目的。

同时，文中方法不仅适用于自身结构部位的频率敏感性分析，还可用于计算分析风机运营期间产生的频率偏移现象，经过算例分析，得到如下结论和建议：

（1）塔筒的变截面特性对风机系统一阶横向自振频率的准确预估意义较大。考虑桩、土相互作用及塔筒变截面特性后，计算结果与实测数据基本吻合，具有良好的一致性。

（2）土体模量较小时，风机系统横向自振频率受土体模量影响很大，自振频率随土体模量的增加呈现较大的增长；
土体模量大于40 MPa后，土体模量的变化对风机系统横向自振频率影响较小，几乎可忽略不计。

（3）在相同变化量的情况下，塔筒的高度变化对风机系统横向自振频率有较大的影响，连接段高度对风机系统横向自振频率影响较弱。海床平面高度变化通过影响桩的埋入深度与连接段的高度，进而导致风机系统横向自振频率产生较大的偏移。在桩径与桩长的选取上，大直径柔性长桩具有较好的抗频率偏移能力。

（4）文中方法适用于土体模量恒定的超固结土、岩层及土层厚度不大的砂土，但对于土层厚度较大的砂土层，土体模量随深度线性增加的现象不可忽视，应对土层模量适度加权或作分层处理。