近年来，我国海上风电事业发展迅猛，“十三五”规划要求到2020年，全国海上风电开工建设规模达到1×10⁷ kW^[1]. 海上风电开发逐渐迈向30 m以上深水区，多采用导管架基础. 相比于在浅水环境中被广泛采用的大直径单桩基础，导管架基础具有波流荷载小、变形小、经济性好等诸多优点^[2]. 导管架基础多为群桩结构，在上部传来的倾覆弯矩的作用下，各桩竖向拉拔作用明显^[3]. 在我国东南沿海地区台风频发^[4]，在风、浪、流荷载的耦合作用下，桩-土界面承受大幅循环剪切，桩周土可能产生不可恢复的塑性应变，桩-土相互作用持续弱化，直接威胁海上风机的运行安全.

在桩头竖向力、水平力、弯矩的耦合作用下，桩侧土对桩的侧向作用、竖向作用以及桩端土对桩端的竖向作用可分别用p-y、t-z、Q-z弹簧模型描述^[5]. Wei等^[6]研究p-y弹簧对循环荷载下导管架基础海上风机动力响应的影响，并与无桩-土相互作用的工况进行对比. Abhinav等^[7-8]研究在随机风、浪荷载条件下，导管架基础海上风机在无桩-土相互作用、弹性桩-土相互作用^[5]2种工况下的频率、极限承载力以及动力响应差异. 上述研究均表明，不考虑桩-土相互作用会导致导管架基础的位移响应被低估. 目前对于循环荷载下导管架基础海上风机动力响应的研究较少，尤其是对考虑桩-土相互作用循环弱化效应的动力响应缺乏认识.

本研究针对导管架基础海上风机，开展极端循环荷载下整体结构的动力响应的计算与分析，在数值计算中分别采用美国石油学会（American petroleum institute，API）规范系列模型（p-y、t-z、Q-z）^[5]、可考虑桩-土界面竖向强度和刚度循环弱化效应的弹塑性t-z模型以及桩-土侧向弹塑性p-y模型. 评估各部分弹簧对泥面处抗倾覆弯矩的贡献以及t-z、p-y弹簧的循环弱化效应对导管架基础顶点处转角、位移动力响应的影响.

如图1所示，导管架基础海上风机由转子、轮毂、机舱、塔筒、导管架结构、桩基、附属结构等部分构成. 海上风机的轮毂高度一般在海平面以上约100 m，转子直径超过100 m，塔筒底部直径为5~6 m，为典型高耸动力敏感型结构. 导管架基础在设计中须满足整机频率、杆件极限及疲劳强度、整体变形等要求^[9].

API 规范中推荐的p-y、t-z、Q-z弹性模型^[5]被广泛应用于海洋工程桩基设计. 该方法较全面地考虑了土体非线性、土层变化、土体种类等影响. 其中，砂土中桩的p-y曲线表达式为

(1) $p = A{p_{\rm{u}}}\tanh\; [KXy/\left( {A{p_{\rm{u}}}} \right)].$

式中：A为循环荷载效应系数，在循环荷载下A=0.9；p为泥面以下深度X m处作用于桩上的水平土抗力；K为土抗力初始模量；y为深度X m处桩的侧向位移；p_u为单位桩长的极限水平抗力. 相关参数确定方法请参考API规范^[5].

砂土中桩的t-z曲线表达式为

(2) ${\tau _{\max }} = {K_0}{\sigma _{{\textit z}}}\tan\; \delta .$

式中：τ_max为桩-土界面极限剪切强度，K₀为地基侧压力系数，σ_z为有效上覆土压力，δ为桩-土界面摩擦角. 桩-土界面竖向抗力τ的发挥与桩-土相对位移z有关，API规范^[5]中给出了τ/τ_max与z的关系曲线.

砂土中桩的Q-z曲线表达式为

(3) $q = {\sigma _{{\rm{tip}}}}{N_{\rm{q}}},$

(4) ${Q_{\max }} = q{A_{{\rm{tip}}}}.$

式中：q为桩端单位面积极限承载力， ${\sigma _{{\rm{tip}}}}$为桩端位置有效上覆压力，Q_max为Q-z弹簧极限承载力，A_tip为桩端面积，N_q为无量纲地基承载力系数. 桩端抗力Q的发挥与桩端相对位移z_tip/D（z_tip为桩端竖向位移，D为桩径）相关，API规范^[5]中给出了Q/Q_max与z_tip/D的关系曲线.

API的p-y模型^[5]未合理考虑循环荷载作用，许多学者针对该问题提出了考虑循环荷载的桩-土p-y模型，主要有p-y曲线循环折减法、土体刚度衰减有限元法、弹塑性增量p-y曲线法. 黄茂松等^[10-12]根据水平循环荷载下单桩离心或室内模型试验，提出考虑循环弱化效应的修正p-y曲线. Achmus等^[13-14]利用土体刚度衰减有限元法对一定次数循环荷载下大直径单桩的水平侧移量进行评估. 苏栋^[15]在边界面模型理论框架下，构建考虑循环荷载效应的增量p-y模型，本研究采用该模型描述循环荷载下的桩-土侧向相互作用，并选取密砂模型（参数见表1）.

对桩-土竖向循环剪切作用，采用边界面模型理论框架^[15-20]下的能够考虑循环弱化效应的桩-土界面t-z模型，模型参数均可通过界面循环剪切试验标定. 如图2所示为模型示意图. 图中，实线为传统的一维边界面模型，虚线表示界面强度的衰减. 在持续的循环剪切作用下，界面附近土体发生剪缩，法向应力不断衰减，界面强度由τ_u0逐渐降低至τ_u，表达式为

(5) ${K_{{\rm{ep}}}} = 1/(1/{K_{\rm{e}}} + 1/{K_{\rm{p}}}),$

(6) ${K_{\rm{p}}} = h{K_{\rm{e}}}[({\tau _{\rm{u}}}/{\tau _{\rm{m}}}) ( \overline \rho /\rho) - 1],$

(7) ${\tau _{\rm{u}}}{\rm{ = }}{\tau _{{\rm{u}}0}}g({{\textit z}^{\rm{p}}}),$

(8) $g({{\textit z}^{\rm{p}}}) \!=\! \max\; \left\{ {\left. {\left[{{\left( {{{\left(1 \!+\! {{\int \frac{\left| {{\rm d}{{\textit z}^{\rm{p}}}} \right|} {{{\textit z}_{50}}}}}\right)}^a}} \right)}^{ - 1}} \!\!-\! b{{\int \frac{\left| {{\rm d}{{\textit z}^{\rm{p}}}} \right|} {{L_{\rm{u}}}}}}\right],\,\frac{\tau _{\rm{r}}}{\tau _{{\rm{u0}}}}} \right\}.} \right.$

式中：K_ep为弹塑性抗力系数；K_e为弹性抗力系数；K_p为塑性抗力系数；h为t-z曲线形状参数；τ_m为强度边界范围内的历史最大界面剪应力； $\bar \rho $=2τ_m为抗力边界的长度；ρ为当前应力点沿应力增量dτ反方向到抗力边界面的距离；τ_u0为截面初始极限强度；g（z^p）为界面强度衰减函数；z^p为塑性位移；a、b为界面强度衰减参数；z₅₀为参考剪切位移，取界面抗力达到极限值一半时对应的界面剪切位移；L_u为单位剪切位移；τ_r为循环界面剪切残余强度.

已有研究^[21-22]表明，等刚度边界可以准确描述桩-土界面法向边界特征. 采用等刚度边界条件下的钢-砂界面剪切试验进行t-z模型参数标定.

1）界面初始极限强度. τ_u0可从界面剪切试验中直接获得，若无相关试验数据，可由式（9）计算得出：

(9) ${\tau _{{\rm{u}}0}} = {K_0}{\sigma _{{\textit z}}}\tan \;\delta .$

2）弹性抗力系数. K_e描述t-z弹簧的弹性行为，表达式为

(10) ${K_{\rm{e}}} = {{{\tau _{{\rm{u}}0}}}/{\left( {2{{\textit z}_{50}}} \right)}}.$

3）曲线形状参数. 曲线形状参数h控制t-z曲线的形状，调整h可使模拟曲线与试验曲线保持基本一致.

4）循环界面残余强度. 在界面循环剪切试验中，界面强度呈现先快速下降后近似线性下降的规律. 在经历足够的循环剪切次数之后，界面强度趋于稳定，可将其视为循环界面剪切残余强度.

5）卸载刚度修正系数. 在界面剪切试验中发现，t-z曲线的初始卸载刚度明显大于初始加载刚度，故须对卸载刚度进行修正. 卸载刚度修正系数B为初始卸载刚度K_un与初始加载刚度K_e之比：

(11) $B = {K_{{\rm{un}}}}/{K_{\rm{e}}}.$

6）界面强度衰减参数. 参数a、b用于描述法向应力衰减，表达式分别为

(12) $a = {\eta _{{a}}}\Delta {\sigma ^{({\rm I})}}/{\sigma _{\rm{n}}}^{(0)},$

(13) $b = {\eta _{{b}}}\Delta {\sigma ^{(\Pi )}}/{\sigma _{\rm{n}}}^{(0)}.$

式中：σ_n⁽⁰⁾为界面初始法向应力，Δσ^(I)、Δσ^(II) 分别为I、II阶段（快速衰减阶段、线性稳定衰减阶段）衰减的法向应力，η_a、η_b为调整系数.

标定后的t-z模型参数如表2所示. 表中，τ_r/τ_u0为循环后界面残余强度比. 模型模拟结果和桩-土界面剪切试验结果如图3所示，图中，d为剪切位移，τ为剪应力. 可以看出，循环弱化t-z模型可以较好地描述钢-砂界面循环剪切弱化现象.

对比分析各部分弹簧对泥面处抗倾覆弯矩的贡献并评估不同p-y、t-z模型对导管架基础海上风机动力响应的影响. 设计算例如表3所示，其中工况1中的3种弹簧模型均采用API规范^[5]中的弹性模型，工况2中仅有t-z模型为弹塑性模型，工况3中仅有p-y模型为弹塑性模型，工况4中t-z、p-y模型均为具有循环弱化效应的弹塑性模型. 3种API规范^[5]模型的参数取值如表4所示. 循环弱化p-y、t-z模型的参数取标定后参数. 利用有限元软件COMSOL 对导管架基础海上风机进行整体建模，结合MATLAB进行二次开发将上述循环弱化p-y、t-z模型写入有限元模型中，对导管架基础海上风机进行极端循环荷载下的动力时程计算与分析.

在数值计算中选用NREL 5 MW风机，风机详细参数见文献[23]，水深为36.9 m，风机基础形式为四桩导管架基础，转子、轮毂、机舱按照实际质量等效为偏心集中质量块，模型详细参数如表5所示.

考虑海上风机处于台风极端荷载工况下的停机状态. Abhinav等^[8]针对平均风速为50 m/s和脉动风湍流强度为10%的条件，联合TurbSim和FAST软件计算100个湍流风速时程下风机塔筒顶部等效风荷载，在考虑安全系数后，风对塔筒顶部的作用力F_max=1.125 MN，本研究选用该值作为等效风荷载. 作用在塔筒顶部的等效循环集中力表达式为

(14) ${F_{{\rm{wind}}}}{\rm{ = }}{F_{\rm{m}}} + {F_{\rm{A}}}\sin\; (2\pi {t_{\rm{c}}}/T).$

式中：T为波浪周期；t_c为计算时间；F_m为定常风荷载，F_A为脉动风荷载，本研究中取F_m=F_A=0.5F_max=0.562 5 MN.

对于波流荷载，利用海工软件SACS进行计算，采用流函数理论并考虑海流作用的影响. 五十年一遇的1%波高H_1%=10.76 m，对应T=10.10 s，为了保证计算精度，波浪流函数取9阶. 波浪对导管架结构的拖曳力系数C_D=1.2，惯性力系数C_M=2.0^[24]，海流速度为2.15 m/s. 为了方便在有限元模型中施加波流力，采用泥面弯矩等效原则，将作用在导管架结构上的波流力等效为4个集中力，分别施加在导管架结构的4根斜主导管上，加载高度如图1中虚线所示.

沿导管架对角线方向施加风、浪荷载以体现最不利工况，设定波流荷载为不变荷载，风荷载为可变循环荷载. 在数值计算中，为了防止突加荷载对风机动力响应的影响，将不变荷载在5 s之内由0线性增大到设定值，之后再施加循环荷载.

对导管架基础海上风机而言，上部风机承受的风荷载、导管架部分承受的波流荷载可等效为作用在泥面位置基础中心点处的集中水平力H和倾覆弯矩M，如图4所示. 水平抗力由p-y弹簧提供，泥面处抗倾覆弯矩由p-y、t-z、Q-z弹簧共同提供，表达式分别为

(15) ${M_{{{py}}}} = \sum\nolimits_{j = 1}^4 {\sum\nolimits_{i = 1}^N {{P_{{}_{ji}}} {h_{{}_{ji}}}} {H_{{}_{ji}}}} ,$

(16) ${M_{{{t{\textit z}}}}} = \sum\nolimits_{j = 1}^4 {\sum\nolimits_{{{i}} = 1}^N {{\tau _{{}_i}} \pi D {h_{{}_i}}} {L_{{}_j}}},$

(17) ${M_{{{Q{\textit z}}}}} = \sum\nolimits_{j = 1}^4 {{Q_{{}_j}} {L_{{}_j}}} .$

式中：M_py为p-y弹簧力对泥面处的抗倾覆弯矩，M_tz为t-z弹簧力对基础中心点的抗倾覆弯矩，M_Qz为Q-z弹簧力对基础中心点的抗倾覆弯矩，j为桩编号，N为每根桩划分的段数，P_i为第i桩段单位长度的水平力，h_i为第i桩段的长度，H_i为第i桩段中心点到泥面的竖向距离，τ_i为第i桩段桩侧单位面积的竖向力，L_j为第j根桩中心点到基础中心点的水平距离，Q_j为第j根桩的桩端力.

针对工况1，分析各部分弹簧力（p-y、t-z、Q-z）对导管架基础泥面处抗倾覆弯矩贡献的时程响应. 如图5所示为在静荷载施加结束5 s之后的计算结果. 图中，W为各部分弹簧对泥面处抗倾覆弯矩贡献占比. 可以看出，t-z、p-y弹簧的贡献占比波动较大，且t-z弹簧的贡献明显大于p-y弹簧，两者基本关于0.5贡献占比参考线上下对称. 在风机受循环荷载的过程中，p-y弹簧力在泥面位置形成的弯矩可能为负，说明此时浅部土体的水平抗力形成的负弯矩大而深部土体的水平抗力形成的正弯矩小，p-y弹簧对风机的抗倾覆发挥负面作用. Q-z弹簧贡献占比小于0.01，可以忽略不计. 因此，在3种弹簧模型中，t-z弹簧对导管架基础抗倾覆性能发挥着最重要的作用，应在工程设计中予以重点关注.

海上风机基础顶点处的转角是风机能否正常运行的重要评价指标. 根据挪威船级社海上风机设计规范^[9]，25 a服役期内基础顶点处的累积转角须小于0.25°. 本研究重点关注不同桩-土相互作用对导管架基础顶点处的转角、位移动力响应的影响.

如图6所示为导管架基础顶点处的转角响应时程图. 图中，R为基础顶点处转角. 在循环荷载作用下，1）工况1转角响应最小且最快达到稳定，工况2转角响应最大且稳定较慢，说明弹塑性t-z弹簧的循环弱化对基础顶点处转角响应的影响显著，使转角增大109.4%；2）工况3转角响应与工况1基本一致，说明弹塑性p-y弹簧的循环弱化对基础顶点处转角响应的影响不明显，使转角增大1.2%；3）工况4转角响应介于工况1、2之间，t-z、p-y弹簧的循环弱化耦合作用使该点处的转角增大66.1%，说明相比于仅考虑t-z弹簧的循环弱化，同时考虑t-z、p-y弹簧的循环弱化可使基础顶点处的转角响应有所减小. 这是由于桩基的水平抗力由p-y弹簧提供，工况4中弹塑性p-y弹簧随荷载循环不断弱化，导管架结构平动量明显大于工况2. 工况2中导管架结构转动更加明显，工况4中导管架结构平动更加明显，受力模式的不同导致工况2转角响应大于工况4.

如图7所示为导管架基础顶点处水平位移响应时程图. 图中，u为基础顶点处水平位移. 可以看出：1）工况1中水平位移响应幅值最小且最快达到稳定；2）工况2中弹塑性t-z弹簧的循环弱化导致该点处水平位移增大135.4%；3）工况3中弹塑性p-y弹簧循环弱化导致该点处水平位移增大 60.7%；4）工况4中弹塑性p-y、t-z弹簧的循环弱化耦合作用使得该点处水平位移增加更加明显，为247.8%. 如图8所示为导管架基础顶点处竖向位移响应时程图. 图中，w为基础顶点处竖向位移. 在循环荷载作用下，1）工况1竖向位移响应最小且最快达到稳定；2）工况2中该点竖向位移响应有较明显的增加，t-z弹簧的循环弱化导致该点竖向位移增大8.9%；3）工况3中竖向位移响应与工况1基本一致，说明p-y弹簧的循环弱化对基础顶点处竖向位移响应的影响不明显；4）工况4 中竖向位移响应最大且无稳定趋势，说明t-z、p-y 弹簧的循环弱化耦合作用会大大增加基础顶点处的竖向位移响应，使得该点处竖向位移增大158.1%.

（1）导管架基础海上风机在承受极端循环荷载的过程中，t-z弹簧对泥面处抗倾覆弯矩的贡献最大，p-y弹簧对泥面处抗倾覆弯矩的贡献可能为负，其主要作用为抵抗水平力，Q-z弹簧对泥面处抗倾覆弯矩的贡献较小，可忽略不计.

（2）相比于API规范桩-土相互作用弹性模型，考虑桩-土竖向相互作用的循环弱化效应，即考虑t-z弹簧的循环弱化效应会大大增加基础顶点处的转角响应.

（3）若同时考虑桩-土竖向、侧向作用的循环弱化，即t-z、p-y弹簧循环弱化的耦合作用，导管架基础的平动更为明显，基础顶点处的位移响应明显增大，但转角响应相比于仅考虑t-z弹簧循环弱化的情况有所减小. 因此，在导管架基础海上风机设计中，须充分考虑桩-土作用的循环弱化效应.