基于三维地震作用下桩-土耦合效应近海超大型风力机动力学响应研究

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

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

海上风能因风速稳定、可开发区域广及对人类影响小等优势逐渐成为最具增长潜力的可再生能源之一[1]。至2019年，全球海上风电新增装机容量27.2 GW，另据国际能源署(International Energy Agency)预测，自2018—2040年，全球海上风电装机容量将增加15倍[2-3]。

为降低度电成本，海上风力机呈大型化发展趋势，其中因技术成熟、安装简单及成本低等优点，多采用固定式支撑结构风力机。其支撑结构可分为单桩式、重力式、三桩式及导管架式。因单桩式支撑结构具有安装简单和成本低等优势，成为近海风力机中应用最广的支撑结构[4]。然而，随单桩式风力机逐渐大型化，其各组成结构尺寸显著增加，多采用降低单个部件重量以减小安装难度，此虽增加了结构柔性但也更易受地震激励影响[5]。由于单桩式风力机多安装于地震频发且海床土质松软的东南沿海地区风场，其地质特征导致桩-土耦合效应明显，桩基极有可能因地震时的土体剧烈震动而遭破坏，导致支撑结构发生局部失稳及属高耸结构的风力机塔架振动加剧。因此，研究准确的桩-土耦合效应模型及地震作用下超大型单桩式风力机动力学响应尤为重要。

单桩式海上风力机桩基直径达9 m，属于超大直径受荷桩，其桩-土耦合效应广受重视。徐海滨等[6]基于莫尔-库伦模型建立土体模型，并模拟了桩径为6 m的大直径桩的桩-土耦合效应。Yang等[7]基于有限元法及现场水平载荷试验进行对比分析。Min等[8]等采用有限元软件建立砂土土体有限元模型，在此基础上研究了桩的受力特点。Byrne等[9]以大长径比的桩为研究对象，提出一种新的桩-土耦合效应模型。Abdullahi等[10]以REL 5 MW单桩式风力机为研究对象，基于分布式弹簧和p-y曲线建立桩-土耦合效应模型，并将其与四弹簧桩-土耦合效应模型对比，发现二者虽有一定误差，但误差较小。

以上针对桩-土耦合效应模型的研究均未考虑地震载荷，由上文可知，地震对近海风力机造成较大威胁。为此，国内外诸多学者开展了地震作用下风力机动力学响应研究。陈法波[11]通过时域法和频域法研究地震作用下3.6 MW风力机动力学响应，发现频域法因无法考虑结构非线性，致其计算结果误差较大。季亮等[12]基于时域法和频域法分析了NREL 5 MW风力机在地震作用下的塔底应力、弯矩及剪力，结果表明：频域法计算结果较时域法偏低。因此，采用时域法分析地震作用下风力机动力学响应十分重要。

地震激励作用下风力机结构响应是准确计算结构承载力和施工成本的重要依据，而桩-土耦合效应是评估地震作用下结构响应的重要因素[13]。刘中胜等[14]基于Wolf方法建立桩-土耦合效应模型，并通过风力机多体动力学仿真软件FAST模拟考虑桩-土耦合效应时风力机在不同地震强度下结构响应，时程结果表明：桩-土耦合效应对风力机在地震作用下的结构响应影响显著，不可忽略。Ma等[15]对比分析地震下是否考虑桩-土耦合效应的3 MW海上单桩式风力机瞬态动力学响应研究发现：桩-土耦合效应明显改变塔架一阶固有频率及塔架振动响应。因此，对于地震作用下风力机动力学响应，必须考虑桩-土耦合效应。

基于p-y曲线法建立的桩-土耦合效应模型因简便、具有一定精度且可较好地描述水平地震作用下桩-土耦合效应，在学术界得到广泛应用[16]。庞文彦[17]以1.5 MW风力机为研究对象，对比分析现有计算桩基水平承载力方法，结果表明：p-y曲线可较好地描述桩基与土体间相互作用。王明超等[18]基于线性p-y曲线确定三种描述桩-土耦合效应的弹簧模型刚度，并研究其对单桩式风力机前两阶固有频率影响，但线性p-y曲线尚无法体现土反力与桩基挠度间非线性关系。Lesny等[19]基于p-y曲线研究5 MW风力机在风、浪载荷作用下的结构动力学响应表明：当桩基受载荷作用产生挠曲时，桩周土反力与桩基挠度呈非线性关系。Reese[20]通过随土深变化的非线性p-y曲线模拟风力机桩基埋于砂土中的桩-土耦合效应，然未考虑地震作用影响。Ju等[21]以非线性p-y曲线建立桩-土耦合效应模型，并以此为基础研究NREL 5 MW 桁架式风力机在风、浪和地震作用下的动力学响应，结果表明：地震作用于风力机时，必须考虑桩-土耦合效应；此外，地震加速度峰值大于0.52g时，较风浪载荷，地震为风力机控制载荷。

显然，地震对结构的空间影响是三维的，地震不仅引起水平方向土体运动，也将引发垂向土体运动。但，上述研究未考虑地震作用影响或仅分析水平而忽略垂向地震作用下桩-土耦合效应，无法完整反映三维地震作用下桩-土耦合效应。Kjørlaug等[22]对比分析是否考虑垂向地震作用下桩-土耦合效应的陆上风力机动力学响应，发现考虑垂向地震激励下桩-土耦合效应使塔架加速度响应放大2倍～5倍，但未对海上风力机垂向地震作用下桩-土耦合效应展开研究。此外，随风力机日趋大型化，而目前针对超大型10 MW单桩式风力机研究较少，且上述文献均仅以地震加速度数据为地震作用于风力机，其未能真实反映地表运动，而能准确反映地表运动的地震位移数据作为地震作用于风力机却鲜有研究。

为此，选取实测地震位移数据为地震作用于风力机，基于p-y、t-z及q-z曲线建立桩-土耦合效应模型，并以超大型DTU 10 MW单桩式近海风力机为研究对象，对比分析不同桩-土耦合效应下动力学响应，为风力机结构抗震设计提供参考。

1 地震作用

为准确描述地震发生时地表运动，选取Pacific Earthquake Engineering Research(PEER)[23]提供的1999年发生于台湾里氏7.62级地震位移数据，地震位移峰值为0.146 m，地震位移时域曲线如图1所示，图中X为流向，Y为横向，Z为垂向。

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

2 模型建立

2.1 超大型单桩式风力机

选取近海超大型DTU 10 MW单桩式风力机为研究对象，如图2，主要参数见表1[24]。分析模型中，采用壳单元建模并将叶片、轮毂及机舱等按实际质量等效为塔顶偏心质量点，该模型采用结构化网格，网格尺寸为0.5 m，网格数为25 074，节点数为25 423，边界条件为流向、横向及垂向地震发生时地表位移数据，此外，鉴于油漆、焊接及法兰等因素，计算模型有效密度取8 500 kg/m3。以瑞利阻尼反映结构阻尼，瑞利阻尼系数按下式确定

(1)

(2)

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

图2 风力机及桩-土耦合效应示意图Fig.2 Schematic diagram of wind turbine and pile-soilcoupling effect

表1 DTU 10 MW风力机主要参数Tab.1 Main parameters of DTU 10 MW wind turbine

2.2 桩-土耦合效应模型

地震发生时，水平和垂直方向地震波通过土层传递并作用于结构，结构由于振动所产生水平及垂直方向惯性力反作用于土层，称为三维地震作用下桩-土耦合效应[26]。本文采用美国石油协会(American Petroleum Institute，API)推荐的p-y、t-z及q-z曲线弹簧模型建立三维地震作用下桩-土耦合效应模型，分别模拟水平及垂直方向地震作用下的桩-土耦合效应，即桩侧土水平抗力与桩挠度间关系、桩周土摩擦力及桩周土摩阻力与桩沉降距离间关系，该模型因综合考虑土体非线性、不同土体种类及各土层特性等因素广泛应用于近海工程[27]。桩基长30 m，以1 m及5 m为间隔布置沿桩身方向和垂直于桩的弹簧模拟桩侧土水平抗力和桩周土摩擦力，桩端则通过垂向弹簧模拟桩周土摩阻力，选取非线性弹簧单元(Combin 39)模拟弹簧约束，见图2。该弹簧单元在一、二及三维应用中具有拉伸和扭转能力且可模拟力与变形间非线性关系轴向单元，其需定义力-变形(F-D)曲线以确定弹簧受力，如图3所示。

图3 非线性弹簧单元力-变形曲线Fig.3 Nonlinear spring element force-deformation curve

如图3所示，非线性弹簧单元将根据输入的离散点以拟合成力-变形曲线，从而实现弹簧非线性变形。

为观察桩周土摩阻力及三维地震作用下桩-土耦合效应对风力机的影响，以单桩式风力机为基础建立三种桩-土耦合效应模型，详情如表2所示。

表2 不同桩-土耦合效应模型

2.2.1p-y曲线

按API推荐，砂性土p-y曲线表达式为

p=Aputanh[kHy/(Apu)]

(3)

pus=(C1H+C2D)γH

(4)

pud=C3DγH

(5)

pu=min{pus,pud}

(6)

式中:pus为浅层土壤极限土抗力；pud为深层土壤极限土抗力；C1、C2、C3为砂砾摩擦角函数；D为桩径；γ为土体有效容重；H为土体深度；A为经验参数，循环载荷时取0.9，定载荷时取(3～0.8H/D)≥0.9；pu为深度为H时土体极限抗力；y为桩基挠度；k为初始土体切变模量。土壤参数如表3所示。

表3 土壤参数

p-y曲线计算结果如图4所示。

2.2.2t-z曲线

砂性土t-z曲线可由式(7)、(8)和表4得到，其表达式为

fs=KP0tanδ

(7)

tmax=fs·AS

(8)

式中:K为无因次土层侧压力系数；P0为计算点处有效复土压力，P0=γH，γ为土体有效容重；δ为桩土间摩擦角；AS为桩身面积；fs为单位面积极限桩周土摩擦力；tmax为极限桩周土摩擦力。

图4 不同土深处p-y曲线Fig.4 p-y curves at different depths

表4 砂性土t-z曲线参数Tab.4 Parameter of t-z curve of sandy soil

表4中:z为桩沉降距离；D为桩直径；t为桩周土摩擦力。t-z曲线计算结果如图5所示。

图5 不同土深处t-z曲线Fig.5 t-z curves at different depths

2.2.3q-z曲线

将式(9)、(10)及表4联立可得q-z曲线，计算式为：

qp=NqP0

(9)

qmax=qp·Ap

(10)

式中：Nq为无量纲承载力系数；qp为桩端单位面积极限阻力；Ap为桩端横截面积；qmax为极限桩周土摩阻力。

表5 砂性土q-z曲线参数

计算q-z曲线结果由图6所示。

图6 q-z曲线Fig.6 q-z curves

2.2.4 三维实体土体模型

为衡量三种桩-土耦合效应模型的准确性，采用基于粘弹性理论的等效线性模型作为土体的本构模型并建立三维实体土体模型以模拟桩-土耦合效应，土体模型如图7所示。

3 验 证

原有整机在进行叶片及轮毂的简化后，其将造成一定误差，为此对不考虑桩-土耦合效应的整机模型与简化模型分别进行模态计算，以此判断简化模型与整机结构动力学特性差异，前四阶固有频率如表6所示。

图7 三维实体土体模型Fig.7 3-D solid soil model

表6 简化模型与整机结构前四阶固有频率Tab.6 Simplified model and the first four natural frequencies of the machine structure

由表6可知，简化模型与整机结构固有频率误差范围较小，均在5%以内。因此，可判断简化模型与整机结构动力学特性在误差允许范围内及二者一致性较好。

4 结果及分析

4.1 模 态

为研究垂向桩-土耦合效应对塔架动力学特性的影响，采用Block Lanczos方法分别对模型1、2和3模态分析，将得到的模态振型和固有频率对比分析。模型1、2和3前四阶模态振型如图8所示，前六阶模态固有频率如表7所示。

(a) 模型1

表7 模型1、2和3前六阶模态固有频率

由图8可知，模型1和2前四阶模态振型相似，模型3前两阶模态振型与模型1及2类似，三阶及四阶模态振型发生较大改变。模型1和2塔架前两阶模态振型发生前后弯曲振动，三阶及四阶模态振型为侧向弯曲振动，而模型3塔架前四阶模态振型均为前后弯曲振动，且模型1和2三阶及四阶模态振型相对位移峰值位于塔架中部，模型3三阶及四阶模态振型相对位移峰值位于塔顶。此外，愈加复杂的桩-土耦合效应下，塔架模态振型振动相对位移峰值愈大。因此，桩周土摩阻力未对塔架前四阶模态振型产生较大影响，仅增加塔架前四阶模态振型相对位移，而垂向桩-土耦合效应改变塔架三阶及四阶模态振型及相对位移峰值位置，同时亦加剧前四阶模态相对位移，且其增幅更大。

由表7可知，较之模型1，模型2和3塔架前四阶固有频率均有不同程度变化，其中模型3变化幅度更大，而三种模型的五阶及六阶固有频率差异较小。故前四阶固有频率受桩周土摩阻力和垂向桩-土耦合效应影响显著，受后者影响更甚，而高阶模态固有频率受二者影响较小。由于未考虑门洞和爬梯等设备，塔架为几乎完全对称结构，因此相邻两阶固有频率较为接近，进一步说明计算结果的准确性。此外，工程实际中，常以塔架一阶固有频率是否避开风轮旋转频率及叶片通过频率(1P：0.230～0.288 Hz及3P：0.779～0.865 Hz)判断塔架是否将与风轮及叶产生共振，通常认为1P～3P为安全区，为更直观地观察塔架一阶固有频率是否在安全区内，将其与安全区范围绘于图9中比较[28]。

由图可知，模型1、2和3塔架均位于安全区内，但较之模型1和2，模型3塔架固有频率安全裕度更大。因此，考虑垂向桩-土耦合效应后，塔架固有频率远离共振频率，更不易与上部结构发生共振。

4.2 动力学响应

4.2.1 塔 架

风力机大型化导致塔架柔性大幅增加，使其在极具破坏力的地震作用下易发生过大振动造成风力机损坏，为监测地震过程中不同桩-土耦合效应下的塔架振动响应，采用求解精度高的New-mark积分法对其进行动力学响应分析[29]，其中塔顶加速度为反映塔架振动的重要特征参数。为判断及分析地震作用下三种桩-土耦合效应模型与三维实体土体模型塔架的动力学响应情况，输出如图10所示的塔架加速度时域曲线。

图9 模型与塔架固有频率安全区对比示意图Fig.9 Comparison diagram of the model and the naturalfrequency safe zone of the tower

图10 塔顶加速度时域曲线Fig.10 Time domain curve of tower top acceleration

由图10可知，模型3与三维实体土体模型的塔顶加速度响应最为接近，因此模型3与实际情况最为接近。此外，不同桩-土耦合效应模型塔顶加速度时域曲线差别较大。模型1、2和3塔顶加速度分别于29.7 s、16.175 s及15.5 s达到峰值，后因结构阻尼和地震能量衰减，塔顶加速度响应逐渐减小，峰值分别为2.82 m/s2、4.21 m/s2及6.31 m/s2，标准差分别为0.54 m/s2、0.55 m/s2及0.82 m/s2。较模型1，模型2及3位移响应峰值分别增大49%和124%，标准差分别增加2%和52%。因此，较之仅考虑水平桩-土耦合效应，桩周土摩阻力和三维地震作用下桩-土耦合效应将导致塔顶加速度响应更早达到峰值，且峰值和波动均增加，其中三维地震作用下桩-土耦合效应增幅更甚。综上所述，若仅采用p-y曲线及p-y、q-z曲线模拟桩-土耦合效应，将导致对地震作用下塔架动力学响应评估不足。

图11、12及13为三种桩-土耦合效应塔顶流、横和垂向加速度响应，由图可知桩周土摩阻力及三维桩-土耦合效应不可忽略。较之横和垂向，流向塔顶加速度响应峰值更大，并流和横向塔顶加速度响应均以0为中线上下振荡。为详细分析流、横及垂向加速度波动情况，输出模型1、2及3的流、横和垂向加速度变化范围如表8所示。

图11 流向塔顶加速度时域曲线Fig.11 Time domain curve of acceleration to the top of flow tower

图12 横向塔顶加速度时域曲线Fig.12 Time domain curve of horizontal top acceleration

图13 垂向塔顶加速度时域曲线Fig.13 Vertical tower top acceleration time domain curve

较之模型1，模型2和3三个方向塔顶加速度响应波动均有所增大，其中横向塔顶加速响应变化范围增幅最大，分别增大209%和317%，且三维地震作用下桩-土耦合效应相较于垂向地震作用中仅考虑桩周土摩阻力，流、横和垂向塔顶加速度响应变化范围均有增大，其中流向塔顶加速度响应变化范围增幅最大，增加169%。

4.2.2 支撑结构

单桩式风力机塔架、叶片、轮毂和机舱等设备完全由支撑结构支撑，而支撑结构在地震作用下易发生应力集中现象，若应力过大将导致其局部屈服甚至结构破坏，使风力机存在倒塌风险，因此支撑结构处的等效应力分析十分重要，输出如图14所示的不同桩-土耦合效应时最大等效应力响应云图。

表8 模型1、2和3流、横和垂向加速度变化范围

图14 模型1、2和3支撑结构最大等效应力响应云图Fig.14 Models 1, 2 and 3 support structures respond tomaximum equivalent stress clouds

由图14可知，桩周土摩阻力和三维地震作用下桩-土耦合效应未改变最大等效应力位置，其为过渡段和入水段连接处。对比不同模型最大等效应力可得，较模型1，模型2和3最大等效应力分别减少283%和566%，且模型2及3应力集中区出现明显缩减，模型3应力集中区最小。

此外，由结构钢的屈服极限为250 MPa，强度极限为460 MPa可知，模型1超过材料强度极限460 MPa，结构发生破坏，模型2和3均未超过材料强度极限，其中模型2等效应力大于屈服极限250 MPa，发生局部屈服。因此，桩周土摩阻力降低了支撑结构最大等效应力，三维地震作用下桩-土耦合效应降幅更大。

4.2.3 桩 基

桩基为桩-土耦合效应直接作用对象，因其深埋土体中对其修复异常困难且花费巨大。但对于桩基分析，尤其是考虑垂向地震作用桩-土耦合效应下的桩基分析常被忽略，为分析垂向地震作用桩-土耦合效应对桩基的影响，对模型1、2及3的桩基进行动力学分析。

由图15可知，模型2和3较模型1桩基位移响应明显增大，模型3位移响应最剧烈。此外，模型1和2地震开始作用后，桩基位移响应迅速增大，达到峰值后，桩基位移响应迅速减小，而模型3桩基位移响应激增至峰值后，桩基位移开始衰减但衰减速度远小于模型1和模型2。由此可见，桩周土摩阻力和三维地震作用下桩-土耦合效应均使桩基位移响应显著增加，三维地震作用下桩-土耦合效应下桩基位移响应最剧烈，同时，考虑三维地震作用下桩-土耦合效应地震后期桩基位移响应衰减速度大幅下降。

图15 桩基位移时域曲线Fig.15 Time domain curve of pile foundation displacement

由图16可知，较之模型1，模型2和3桩基各高度处位移响应均有不同程度增大，模型3增幅更大。同时，模型1、2及3桩基位移响应峰值均位于桩基顶部，且模型1和3桩基位移响应随桩基高度增大而减小。综上所述，桩周土摩阻力及三维地震作用下桩-土耦合效应未改变桩基位移峰值位置，均位于桩顶，但其增大了桩基各高度处位移响应，三维地震作用下桩-土耦合效应增幅最大。同时，桩基位移响应随桩基高度增加而减少。

图16 不同高度处桩基位移响应Fig.16 Pile foundation displacement response at different heights

综上所述，三维地震作用下桩-土耦合效应与仅考虑水平桩-土耦合效应风力机动力学响应有很大不同，若仅以p-y曲线及p-y、q-z曲线模拟桩-土耦合效应，将导致对地震作用下风力机动力学响应的评估严重不足。

5 结 论

基于p-y、t-z及q-z曲线描述桩-土耦合效应，以超大型DTU 10 MW单桩式风力机为研究对象，分析不同桩-土耦合效应下风力机动力学响应，结论如下：

(1) 采用集中质量点模拟叶片等上部结构的简化模型与风力机整机结构动力学特性差异较小；桩周土摩阻力和垂向桩-土耦合效应使前四模态相对位移增加，垂向桩-土耦合效应增幅更大，且桩周土摩阻力未明显改变塔架前四阶模态振型及位移峰值位置，考虑垂向桩-土耦合效应后，三阶及四阶模态振型由侧向弯曲振动变为前后弯曲振动且位移峰值位置由塔架中部移至顶部，此外，塔架固有频率安全裕度更大。

(2) 桩周土摩阻力和三维地震作用下桩-土耦合效应增大了塔顶加速度响应波动及峰值，峰值分别增加49%和124%，标准差分别增大2%和52%，并使得流、横及垂向加速度响应变化范围激增，其中横向增幅最大，分别增加209%和318%，且较桩周土摩阻力，三维地震作用下桩-土耦合效应使流向塔顶加速度响应增幅最大，增大169%；若仅以p-y曲线及p-y、q-z曲线模拟桩-土耦合效应，将导致对地震作用下风力机动力学响应的评估严重不足。

(3) 桩周土摩阻力和三维地震作用下桩-土耦合效应未改变支撑结构最大等效应力位置，位于过渡段和入水段连接处，且其降低了最大等效应力，峰值减小了283%及566%，并减小了应力集中区。

(4) 桩周土摩阻力和三维地震作用下桩-土耦合效应未改变桩基位移峰值位置，位于桩顶；此外，较之仅考虑水平载荷桩-土耦合效应，桩周土摩阻力和三维地震作用下桩-土耦合效应下桩基及其各高度处位移响应均有不同程度增大，其中三维地震作用下桩-土耦合效应下桩基位移增幅更大，且使地震后期桩基位移响应衰减速度大幅下降。