多土层桩-土耦合效应对静止状态下近海超大型风力机地震动力学响应及屈曲的影响分析

2021-09-24李志昊闫阳天岳敏楠薛世成

动力工程学报 2021年9期

李志昊，闫阳天，岳敏楠，李春,2，杨阳，薛世成

(1.上海理工大学能源与动力工程学院，上海 200093； 2.上海市动力工程多相流动与传热重点实验室，上海 200093； 3.Department of Maritime and Mechanical Engineering, Liverpool John Moores University, Liverpool, L3 3AF, UK)

可再生能源在全球节能减排中发挥着至关重要的作用，可再生能源容量系数每提高1%，CO2减排成本平均降低0.7欧元[1]。海上风能因风速高、风向稳定、能量密度高、湍流度低、宜于大型化以及不占用陆地资源等优点，被认为是最具潜力的可再生能源之一[2-3]。根据世界风能论坛数据显示，截至2019年全球累计海上风电装机容量达到27.2 GW，当年新增海上风电装机容量达到5.2 GW[4]。

海上风能虽有诸多优势，但亦面临一些挑战：一方面，为获取更多风能资源、提高风能转化率及降低度电成本，风力机呈大型化趋势发展，塔架和叶片尺寸随之增加，导致其稳定性随柔性增加而大幅削减；另一方面，我国风力机多建于风能资源丰富的东南沿海地区，该地区风速较大且地震活动十分活跃，建于此地的超大型风力机极易受到地震影响，在强破坏力的地震载荷作用(以下简称地震作用)下具有高柔性的塔架剧烈振动，从而导致其发生局部屈曲甚至整机倒塌[5]。因此，对超大型风力机进行地震动力学响应研究及屈曲分析十分重要。

桩-土耦合效应为评估地震作用下风力机动力学响应的重要因素[6]。国内外诸多学者对此展开了研究。Vatanchian等[7]研究了考虑桩-土耦合效应后风力机的地震动力学响应，结果表明桩-土耦合效应会加剧地震作用下风力机动力学响应。Passon[8]研究了不同桩-土耦合效应简化模型，将其总结为3类：底端刚性固定模型、耦合弹簧模型及分布式弹簧模型。Bisoi等[9]利用径向的分布式弹簧模型表示桩-土耦合效应模型，通过p-y(其中p为桩周土水平抗力，y为桩基形变)曲线可求得弹簧刚度。王明超等[10]基于线性化p-y曲线及分布式弹簧模型建立桩-土耦合效应模型，研究其对风力机模态特性的影响，但线性化p-y曲线未能准确反映土反力与桩基挠度间的非线性关系，且未考虑地震载荷。闫阳天等[11]用非线性p-y曲线描述桩-土耦合效应，并在此基础上研究了不同地震作用下风力机动力学响应。

上述研究均以单土层地基模型作为海床来模拟桩-土耦合效应，但其无法反映海床分层特性。目前，国内多采用由黏土与砂土组成的简化模型模拟多土层地基模型。盛振国等[12]将多土层地基模型简化为全黏土、全砂土和半黏土半砂土地基模型，并进行了不同土层土参数对风力机结构动力学响应敏感度研究，结果表明桩土参数对风力机动力学响应的影响显著。刘红军等[13]对比分析了纯软土、上软下硬及上硬下软3种地基模型的桩-土耦合效应模型对桩周土水平抗力的影响，发现不同土质条件下桩周土水平抗力差异较大。

以上针对多土层桩-土耦合效应的研究尚存不足：将多土层地基模型简化为由砂土和黏土构成的两土层地基模型不太妥当，因为沿海地区地基因沉积或人类建设活动等因素影响，往往由多土层构成，利用简化后的地基模型难以获得精确结果。

因此，笔者实测地震位移数据为地震载荷，采用我国东海某风电场实测地质数据构建多土层地基模型，同时用纯砂土及纯黏土单土层地基模型分别模拟桩-土耦合效应，并以丹麦科技大学(DTU)研发的近海10 MW超大型单桩式风力机为研究对象，对比分析基于纯砂土、纯黏土和多土层3种桩-土耦合效应下动力学响应，并对其进行屈曲分析，为风力机结构抗震设计提供参考。

1 地震载荷

为准确还原地震作用下的地表运动，选用太平洋地震工程研究中心数据库(PEER)提供的1999年发生于台湾里氏7.62级地震地表位移数据，地震位移时域曲线见图1，其中X为前后方向，Y为侧向，Z为垂向。

图1 地震位移时域曲线Fig.1 Seismic displacement time domain curve

2 模型建立

2.1 DTU 10 MW超大型单桩式风力机

研究对象为近海DTU 10 MW超大型单桩式风力机，如图2所示，主要参数见表1[14]。采用壳单元建模并将轮毂、叶片及机舱等结构简化为置于塔顶的偏心质量点。模型中忽略油漆、焊接及法兰等因素，为避免因模型简化导致的质量损失，计算模型有效密度取8 500 kg/m3。网格尺寸为0.5 m，网格数为23 913，节点数为24 137。

图2 DTU 10 MW超大型单桩式风力机模型Fig.2 DTU 10 MW monopile wind turbine model

表1 DTU 10 MW超大型单桩式风力机主要参数Tab.1 Main parameters of DTU 10 MW monopile wind turbine

以瑞利阻尼反映结构阻尼，瑞利阻尼系数α、β按照下式确定：

α=2ωiωjζ/(ωi+ωj)

(1)

β=2ζ/(ωi+ωj)

(2)

式中：ζ为阻尼比，取1%；ωi、ωj为自振频率，分别取0.2 Hz和0.21 Hz；α为0.002;β为0.05。

2.2 桩-土耦合效应模型

土壤与桩基间的刚度和弹性模量存在明显差异，二者变形能力差别较大，当地震波通过土层作用于桩基时，桩基与土体间发生不同程度变形，所产生的相互作用力称为桩-土耦合效应[15]。通过美国石油协会(API)规范中p-y曲线及弹簧建立桩-土耦合效应模型，该方法因应用简便、可综合反映土体非线性及分层特性等优势得到广泛应用。埋土桩基长30 m，每隔5 m布置垂直于桩身的弹簧以模拟桩周土水平抗力，因其与桩基形变呈非线性关系，故选取非线性弹簧单元(combin 39)模拟弹簧约束，其刚度通过定义力(桩周土水平抗力)与变形(桩基形变)间关系确定。

2.2.1 土壤参数

海床由多土层构成，不同土层土壤参数均不相同。因此，选用位于我国东海海域的江苏某海上风电场通过试桩试验得到的地质资料[16]，土层分布及参数如图3所示，其中γ为土体有效容重，Cu为原状土不排水抗剪强度，φ为内摩擦角。

图3 海床土壤分层示意图Fig.3 Diagram of seabed soil stratification

为对比分析多土层桩-土耦合效应对风力机地震动力学响应及屈曲的影响，选取图3中②层及③层土壤参数分别建立纯黏土及纯砂土单土层地基模型。

2.2.2p-y曲线

API提供了一系列不同土质条件下的p-y曲线，黏土p-y曲线表达式如下：

HR=6D/(γD/Cu+J)

(3)

pu=3Cu+γH+J(HCu/D),0≤H≤HR

(4)

p=pu,y/yc≥8

(5)

p/pu=0.5(y/yc)1/3,0

(6)

p=Aputanh[kHy/(Apu)],0

(7)

yc=2.5εcD

(8)

式中：HR为极限桩周土水平抗力的转折点深度；H为土体深度(以下简称土深)；D为桩径；J为无因次常数；pu为土深H时桩周土水平抗力；yc为50%极限桩周土水平抗力的桩基形变；εc为不排水压缩试验时50%最大应力对应的应变；A为经验参数，循环载荷时取0.9，定载荷时取(3-0.8H/D)且A≥0.9；k为初始土体切变模量。

砂土p-y曲线表达式为：

pus=(C1H+C2D)γH

(9)

pud=C3DγH

(10)

pu=min{pus,pud}

(11)

式中：pus为浅层土壤极限桩周土水平抗力；pud为深层土壤极限桩周土水平抗力；C1、C2、C3为砂砾摩擦角函数。

纯黏土p-y曲线计算结果如图4所示，多土层p-y曲线计算结果如图5所示，纯砂土p-y曲线计算结果如图6所示。

图4 纯黏土不同土深处p-y曲线Fig.4 p-y curve at different depths of pure clay soil

图5 多土层不同土深处p-y曲线Fig.5 p-y curve at different depths of multiple soil layers

图6 纯砂土不同土深处p-y曲线Fig.6 p-y curve at different depths of pure sand soil

3 有限元方法

3.1 模态分析

模态分析通过求解结构自由振动微分方程以获得结构动态特性，结构动态特性包括模态振型及固有频率等。模态分析为最基本的动力学分析，亦为其他动力学分析的基础，如瞬态动力学分析、响应谱分析及谐响应分析等。

假设该结构为线弹性结构，且不考虑阻尼时，结构自由振动微分方程为：

(12)

线弹性结构发生简谐振动时：

u=ficos(wit)

(13)

联立式(12)与式(13)可得：

(14)

式(14)有非零解的充要条件为系数行列式为零，即

(15)

式中:wi和fi分别为第i阶固有频率和模态振型;t为时间。

3.2 瞬态动力学分析

瞬态动力学分析可计算随时间变化的载荷作用下结构动力学响应，其本质为求解结构运动方程。结构运动方程为：

(16)

式中：C为阻尼矩阵；F(t)为载荷向量。

瞬态动力学分析方法主要可分为振型叠加及完全积分法，笔者选取可考虑结构非线性的完全积分法，有限元软件提供的积分方法分为HHT(Hilber Hughs and Taylor)及New-mark积分法，笔者选用求解精度更高的New-mark积分法。

3.3 屈曲分析

结构在外载荷作用下未达到材料强度极限而失效称为屈曲。通过屈曲分析可确定外载荷作用下结构的屈曲振型及临界屈曲载荷。

稳定平衡状态下，结构平衡方程为：

(KL+λKσ)Δu=0

(17)

式中：KL为线性刚度矩阵；Kσ为初应力刚度矩阵；Δu为位移特征矢量；λ为特征值，即屈曲因子。

4 结果及分析

4.1 模态

采用精度高、收敛速度快及适用于大型对称特征值问题的Block Lanczos方法提取风力机塔架的固有频率和模态振型。塔架纵向和侧向前两阶模态固有频率见表2，模态振型见图7。

表2 塔架固有频率

图7 塔架模态振型Fig.7 Tower modal shape

如图7所示，塔架一阶前后及侧向模态振型为前后弯曲振动且相对位移峰值位于塔顶，二阶前后及侧向模态振型为侧向弯曲振动且相对位移峰值位于塔架中部。此外，因未考虑门洞及爬梯等设备，塔架可近似视为对称结构，并由表2和图7可知，一、二阶模态前后及左侧向固有频率和模态振型相近，因此可进一步证明计算结果准确。

4.2 时域

为研究多土层及不同土质桩-土耦合效应对地震作用下风力机塔架动力学响应的影响，输出如图8～图10所示的塔顶位移时域曲线。

图10 塔顶侧向位移时域曲线Fig.10 Time domain curve of lateral displacement of tower top

由图8可知，地震作用时不同桩-土耦合效应下风力机塔顶位移时域曲线随机性较强，且波动形式有略微差别，位移峰值差异较大。地震作用时，多土层、纯砂土及纯黏土桩-土耦合效应下，塔顶位移峰值分别为0.347 m、0.891 m和0.007 27 m，标准差分别为0.387 m、0.118 m和0.002 64 m。相较纯砂土和纯黏土桩-土耦合效应，多土层桩-土耦合效应下塔顶位移峰值分别减小了61.1%和增大了46.73倍，标准差分别减小了69.5%和增大了43.70倍。这是由于黏土内摩擦角小，其土壤颗粒间内摩擦力远小于砂土，土壤刚度也小于砂土，当地震载荷作用于风力机时，纯黏土桩-土耦合效应远小于纯砂土桩-土耦合效应，其因桩基形变而产生的水平反力远小于纯砂土，地震作用时由砂土及黏土组成的多土层桩-土耦合效应小于纯砂土而大于纯黏土。因此，地震作用时多土层桩-土耦合效应下风力机塔顶位移峰值及标准差分别小于纯砂土桩-土耦合效应并大于纯黏土桩-土耦合效应下的风力机塔顶位移峰值及标准差，纯砂土桩-土耦合效应塔顶位移峰值最大，波动最剧烈，多土层其次，纯黏土最小。但与纯黏土桩-土耦合效应相比，地震作用时多土层桩-土耦合效应下塔顶位移更接近纯砂土桩-土耦合效应下塔顶位移。

图8 塔顶位移时域曲线Fig.8 Time domain curve of tower top displacement

由图9和图10可知，地震作用时多土层桩-土耦合效应下塔顶前后和侧向位移均发生较大变化。多土层、纯砂土及纯黏土桩-土耦合效应下塔顶前后位移峰值分别为0.253 m、0.809 m和0.005 74 m，相较纯砂土和纯黏土，多土层桩-土耦合效应下塔顶位移峰值分别减小了68.7%和增大了43.08倍；三者的侧向塔顶位移峰值分别为0.303 m、0.881 m和0.007 07 m，相较纯砂土和纯黏土，多土层桩-土耦合效应下塔顶位移峰值分别减小了65.6%倍和增大了41.86倍。因此，相较纯砂土或纯黏土桩-土耦合效应，多土层桩-土耦合效应下塔顶前后及侧向位移均有不同程度的减小或增大，其中相较侧向，多土层与纯黏土和纯砂土桩-土耦合效应下的塔顶前后位移相差更大。

图9 塔顶前后位移时域曲线Fig.9 Time domain curves of forward and backwarddisplacement of tower top

地震载荷将造成塔架局部较大弯曲，可能会导致局部剪应力激增，若局部剪应力过大，塔架将遭到破坏。图11给出了塔架最大剪应力时域变化曲线。由图11可知，地震作用时不同桩-土耦合效应下塔架最大剪应力峰值和波动幅度存在很大差异。多土层、纯砂土及纯黏土桩-土耦合效应下，塔架最大剪应力峰值分别为23.9 MPa、75.9 MPa和0.556 MPa，标准差分别为6.74 MPa、19.6 MPa和0.143 MPa。相较纯砂土和纯黏土，多土层桩-土耦合效应下塔架最大剪应力峰值减小了68.5%和增大了41.99倍，标准差分别减小了65.6%和增大了46.13倍。显然，纯砂土桩-土耦合效应下塔架最大剪应力响应最剧烈，多土层次之，纯黏土最小。因此，风力机抗地震结构设计时，若以纯砂土或纯黏土单土层计算桩-土耦合效应，将导致对塔架最大剪应力响应过度估计或估计不足。

图11 塔架最大剪应力时域变化曲线Fig.11 Time domain variation curve of maximumshear stress of tower

3种桩-土耦合效应下塔架最大剪应力达到峰值时的响应云图见图12。由图12可知，3种桩-土耦合效应下塔架最大剪应力峰值均位于支撑结构处但集聚区域明显不同。相较多土层和纯黏土桩-土耦合效应，纯砂土桩-土耦合效应下支撑结构处最大剪应力集聚区域面积最大。因此，相较多土层和纯黏土桩-土耦合效应，纯砂土桩-土耦合效应下支撑结构处最大剪应力集聚最严重。

图12 塔架最大剪应力峰值响应云图Fig.12 The peak shear stress of the tower responds

综上所述，地震作用时，若以纯砂土单土层模型模拟桩-土耦合效应将导致对塔架地震动力学响应预估不准确。

4.3 频域

因塔架动力学时域响应为非线性载荷(地震载荷)作用结果，需通过频域分析进一步研究不同桩-土耦合效应下塔架振动特性。采用快速傅里叶变换(FFT)将塔顶位移时域转换为频域数据，如图13和图14所示。

图13 塔架前后振动频域特性Fig.13 Frequency domain characteristics of before andafter vibration of tower

图14 塔架侧向振动频域特性Fig.14 Frequency domain characteristics of lateralvibration of tower

由图13和图14可知，对于3种桩-土耦合效应，塔架前后及侧向一阶模态的固有频率处均存在明显的位移峰值，且纯砂土桩-土耦合效应下塔顶位移峰值最大，多土层次之，纯黏土最小。因此，地震作用时，纯砂土桩-土耦合效应下塔顶振动最强，多土层略弱，纯黏土最弱。此外，3种桩-土耦合效应下塔架一阶模态均被地震载荷诱发。

4.4 屈曲

目前，大型风力机塔架多采用钢制圆锥筒壳结构，该结构易发生局部屈曲，因此有必要对塔架进行屈曲分析。对塔架进行屈曲分析，发现塔架末端为塔架发生局部屈曲的区域，若该区域网格过大易导致有限元计算结果不准确，需细化该区域的有限元网格以获得较精确的计算结果，但网格细化将耗费巨大的计算资源和时间。因此，采用子模型方法将局部网格细化，网格尺寸为0.1 m，网格数为38 267，节点数为38 509。将塔架末端最后15 m作为子模型，对其进行屈曲分析。

表3给出了地震作用时3种桩-土耦合效应下塔顶位移达到峰值时子模型前六阶屈曲因子。

表3 地震作用时3种桩-土耦合效应下子模型前六阶屈曲因子

由表3可知，地震作用时3种桩-土耦合效应下塔架子模型屈曲因子均随阶次升高而逐渐增大，此外，纯砂土桩-土耦合效应下塔架子模型屈曲因子最小，多土层略大，纯黏土最大。因此，桩-土耦合效应将影响地震作用下塔架屈曲因子，不同桩-土耦合效应下塔架屈曲因子不同。考虑纯砂土桩-土耦合效应的地震载荷最接近塔架临界屈曲载荷，多土层次之，采用纯黏土模拟桩-土耦合效应的地震载荷与临界屈曲载荷相差最大。

根据风力机发电机组安全要求[17]，当稳定性分析为特征值屈曲分析时，应附加安全系数1.35，而对地震作用时3种桩-土耦合效应下塔架屈曲因子附加1.35后的安全系数仍大于1，即塔架初始载荷未达到临界屈曲载荷，可判断地震作用时3种桩-土耦合效应下塔架未发生失稳，满足屈曲稳定性设计要求。

为分析塔架局部屈曲区域，输出塔架子模型前三阶屈曲振型，如图15所示。由图15可知，随着阶次升高，塔架局部屈曲区域面积有一定程度的增大，且集中在塔架下端，此处为塔架与支撑结构连接处，二者厚度存在一定差异，属于风力机最薄弱区域，受地震作用时易发生局部屈曲失稳。因此，在风力机抗地震结构设计时，应重点关注该处屈曲区域，通过采用一定优化设计，增强塔架整体结构的抗屈曲能力。