目前，对于近海40 m水深以内的风电机组，大直径单桩是应用最广泛的基础形式[1]。但在复杂的海洋动力环境下，海床沉积物会出现大规模的迁移，引起基础周围的整体冲刷，此外海上风电单桩基础的扰流作用会导致桩周附近床面剪切应力的显著增大，波浪和潮流可能会在桩基础附近诱导产生严重的局部冲刷问题[2-3]。冲刷的形成会导致以下不利后果：1) 桩基有效埋深减小而荷载力臂增大，极大削弱单桩承载能力；2) 上覆土压力减小促使土体固结状态改变，改变土体的承载性能；3) 风机基础固有频率和桩-土相互作用改变，使风机整体动力响应特性更加复杂[4]。因此，准确评估桩周土体冲刷对单桩基础水平承载能力及动力特性的影响，对海上风电单桩基础承载力设计具有重要意义。
为定量评估桩周土体冲刷对单桩承载力的影响，学者们分别采用模型试验法和有限元法开展研究工作。Qi等[4]采用离心机试验讨论了不同平衡冲刷深度情况下，整体冲刷和局部冲刷对单桩水平承载力的影响，并从工程实用角度出发，将桩基有效埋深的概念引入传统p-y曲线法中，评估冲刷的影响。Chortis等[7]关注到风机刚性短桩基础与柔性桩之间的差异，基于离心机试验，研究了不同冲刷形式(局部窄形冲刷、局部宽形冲刷和整体冲刷)对刚性短桩承载力的影响。Li等[8]进一步基于开发的复合加载离心机试验系统，讨论了冲刷形式、冲刷深度和荷载组合情况对大直径单桩的影响。Malakshah等[9]则基于离心机试验，讨论了冲刷效应对单调和循环加载条件下的水平受荷单桩的影响。上述离心机试验方法虽然能够考虑冲刷对土体固结状态和桩基有效埋深的综合影响，但对试验设备和操作技术要求较高，因此相关研究成果尚不丰富。
数值分析方法作为试验方法的有益补充，能够节省时间和试验成本，同时进一步揭示试验模型内部土体流动状态和桩身变形模式。为此，Mostafa[10]、Chortis等[7]和Li等[11]相继建立冲刷平衡后的桩-土相互作用有限元分析模型，针对砂土海床中的水平受荷桩，研究冲刷深度、冲坑形状、基础埋深、加载高度等因素对单桩水平承载力的影响。此外，为反映冲刷对剩余土体固结参数的影响，Lin等[12]推导了冲刷后剩余土体静止土压力系数K0、有效重度γ′、土体峰值内摩擦角φ′p和地基反力模量Kpy的关系式，通过将考虑冲刷固结效应的土体参数(γ′，φ′p，Kpy)代入Reese等[13]提出的p-y曲线法的极限土反力表达式中，形成考虑冲刷固结效应的修正p-y曲线法。目前，针对桩周冲刷影响下的单桩水平承载力开展的研究大多讨论冲刷对单桩静承载力的影响。海上风电单桩基础在复杂动力环境荷载作用下受到湍流风和不规则波浪等循环荷载作用，对于桩周土体冲刷效应对大直径单桩循环受荷特性及动力响应的影响，目前仍缺乏机理性的认识和快速量化的方法。
为此，本文基于Abaqus有限元软件建立冲刷平衡后的桩-土相互作用有限元分析模型，并采用考虑粒间应变张量的土体亚塑性动本构模型，探讨实际湍流风和不规则波浪作用下，整体冲刷和局部冲刷效应对海上风电大直径单桩基础水平承载力响应的影响，揭示整体冲刷和局部冲刷平衡后水平受荷桩的承载机理，为海上风电单桩基础在冲刷效应影响下的长期在位稳定性评估提供方法参考。
1 有限元分析模型
1.1 桩-土相互作用模型
本文以NREL(National Renewable Energy Laboratory, 2020)[14]发布的IEA-15MW海上风力机大直径单桩基础为研究对象，基于Abaqus软件建立预设冲刷坑的桩-土相互作用有限元分析模型，如图 1所示。根据模型对称性，仅取1/2圆柱形区域进行建模，单桩外径D统一取10 m，基础埋深L=45 m，加载高度选在泥面以上30 m处(对应静水面位置)。桩底至泥面以上5 m，桩体壁厚t1=55.34 mm，此后自下而上每隔5 m壁厚变化一次，分别为t2=53.45 mm，t3=51.51 mm，t4=49.527 mm，t5=47.42 mm，t6=45.52 mm。为增大桩-土相互作用有限元分析模型中桩单元的特征长度，提高模型求解的稳定时间增量和计算效率，有限元模型建立过程中保持钢管桩外径与实际一致，根据桩截面抗弯刚度不变的原则，将钢管桩截面等效为实体桩计算，取钢材弹性模量Ep=2.1×1011 Pa，Poisson比νp=0.3。
预设锥形冲刷坑尺寸由平衡冲刷深度Sd和泥沙休止角φs确定，本研究取φs=32°。为模拟冲刷导致的土体卸载及其引起的土体固结状态的改变，模拟过程分为3步：首先，对完整桩-土相互作用模型指定初始应力场，并施加重力荷载，以达到地应力平衡状态；其次，将桩周预设冲刷坑内的土体单元移除，模拟冲刷卸载过程；最后，对大直径单桩施加风荷载和波浪荷载，获取基础的动力响应。
数值分析中，桩基础和土体均采用C3D8实体单元进行模拟，为避免土体变形过大时网格畸变导致计算过程中断，分别对冲刷坑水平扩展范围、泥面至冲刷坑底部1D及基底上下1D深度范围内的网格进行加密，以局部平衡冲刷深度Sld=1.0D情况下单桩归一化荷载-位移曲线作为验证目标，设置了如表 1所示的3组网格划分策略开展网格收敛性验证。
表 1 桩-土相互作用有限元分析模型网格划分方法 |
| 工况 | 细网格 | 偏置网格 (最大网格) | 粗网格 | 数量/个 |
| 1 | 0.17D | 1.5D | 0.40D | 10 128 |
| 2 | 0.10D | 1.3D | 0.25D | 23 152 |
| 3 | 0.05D | 1.0D | 0.25D | 48 304 |
注：D为单桩外径。 |
计算结果如图 2所示。由图可知，当细化网格尺寸取0.1D，总网格数取23 152个时，进一步加密网格，对计算结果的影响小于1%，由此验证了网格的收敛性。因此最终选取细网格尺寸0.1D，粗网格尺寸0.25D的网格划分方法开展后续计算。
1.2 砂土循环本构模型
基于王磊[19]在饱和的福建标准砂中开展的静力三轴压缩试验和单向循环三轴压缩试验，标定土体的亚塑性本构模型计算参数，验证开发的砂土亚塑性本构模型的有效性，同时为后续桩-土相互作用有限元分析模型提供必要的土体参数输入。福建标准砂为级配不良的均质细砂，土体比重Gs=2.633，中值粒径d50=0.17 mm，最大孔隙比为emax=0.952，最小孔隙比emin=0.607，三轴试验中土体相对密度Dr=70%，土体有效重度γ′=9.5 kN/m3。除了上述基本物理特征参数，还须确定13个亚塑性本构模型计算参数。详细标定过程可参考Hong等[20]研究。由此确定的福建标准砂亚塑性模型计算参数如表 2所示。
表 2 福建标准砂亚塑性本构模型计算参数 |
| 功能分类 | 计算参数 | 数值 |
| 八参数亚塑性模型参数 | φ′ /(°) | 35.7 |
| hs/MPa | 3 060 | |
| n | 0.134 | |
| ec0 | 0.952 | |
| ed0 | 0.607 | |
| ei0 | 1.095 | |
| α | 0.11 | |
| β | 3.5 | |
| 粒间应变张量参数 | R | 1×10-4 |
| mR | 2.8 | |
| mT | 1.4 | |
| βR | 0.1 | |
| χ | 3.5 |
注：φ′为土体临界内摩擦角；hs为颗粒硬度；n为等比压缩参数；ec0为土体临界孔隙比；ei0为土体上限孔隙比，取值近似emax; ed0为土体下限孔隙比，取值近似emin; α为土体峰值强度控制参数；β为土体达到峰值强度时的应变控制参数；R为弹性范围；mR为应变方向改变180°时初始剪切刚度，mT为应变方向改变90°时初始剪切刚度；βR和χ为2个刚度随应变弱化速率参数。 |
1.3 有限元分析模型验证
表 3 单桩水平受荷离心机试验参数设置 |
| 参数 | D/m | L/D | e/m | Ep·Ip/(kN·m2) |
| 数值 | 5.9 | 9.3 | 23.0 | 1.1×107 |
注：D为单桩外径，L/D为单桩长径比，e为加载高度，(Ep·Ip)为单桩抗弯刚度。 |
对于循环加载的条件，开展与Zhu等[21]离心机试验一致的循环加载试验，加载频率f=0.05 Hz，循环加载比Hcyc/Hult=0.1~0.4，由此得到的水平循环受荷单桩y/D-H/γ′D3曲线，如图 3c所示，其中C1-1、C1-2、C1-3和C1-4分别对应Hcyc/Hult=0.1、Hcyc/Hult=0.2、Hcyc/Hult=0.3和Hcyc/Hult=0.4的情况。对比文[21]的离心机试验与有限元分析结果(见图 3b)可知，本研究建立的桩-土相互作用有限元分析模型能够准确预测不同Hcyc/Hult的单桩累积变形和水平反力。因此，后续针对冲刷效应影响下的大直径单桩风机基础动力响应分析将在该桩-土相互作用模型基础上，考虑冲刷地形和土体固结状态的变化，开展分析计算。
2 循环荷载计算
波浪荷载的计算是通过在Abaqus中采用Timoshenko梁单元(即B21类型单元)建立风机塔架和单桩基础的有限元模型，将波浪看作平稳的随机过程，基于线性波浪叠加理论和JONSWAP谱生成有限多个不同周期和不同随机初始相位的波浪时程[26]，并将计算出的不规则波浪通过编辑inp. 文件进行定义，Abaqus/Aqua模块将基于如下所示的Morison公式计算出单桩上作用的波浪荷载：
其中：FH为单桩上作用的波浪荷载；ρ为海水密度；d为水深，本文取30 m；η为波幅；u为水质点速度；t为时间；z为水质点所在的垂直位置坐标，原点(z=0)取静水面，以向上为正，水深为d时，海床位于z=-d; $\dot{x}$ 为单桩运动速度；$\ddot{x}$ 为单桩运动加速度；CD和CM分别为拖曳力系数和惯性力系数，本文根据文[25]取CD=0.65、CM=1.6；Cm为附加质量系数，其大小等于CM-1.0。
为综合考虑风机叶片气动特性和伺服控制策略对风荷载的影响，借助开源的FAST耦合仿真软件计算叶片和塔架上作用的气动荷载。计算过程中，由FAST耦合仿真前处理软件Turbsim根据湍流风谱(如IEC Kaimal谱)的逆Fourier变化得到具有随机特性的湍流风场；风机叶片上的气动荷载主要基于广义动态入流理论和叶素动量理论计算；塔架上的风荷载是基于塔架各点的局部风速和气动特性，计算塔架不同位置的气动阻力，同时考虑风切变、湍流及风轮尾流效应，参考文[27]取塔架阻力系数为1.2；风机伺服控制系统则是通过调用外部动态链接库执行控制策略，实现风机变桨控制、偏航控制和转速调节，并间接影响风机气动荷载计算。
表 4 大直径单桩风机基础动力响应计算工况 |
| 工况 | Vhub/(m·s-1) | Hs/m | Tp/s | Sd/D | L/D |
| NS-1 | 10.66 | 1.54 | 7.65 | 0 | 4.5 |
| NS-2 | 0 | 3.5 | |||
| GS-1 | 1.0 | 4.5 | |||
| LS-1 | 0.5 | 4.5 | |||
| LS-2 | 1.0 | 4.5 | |||
| LS-3 | 1.5 | 4.5 |
注：Vhub为风机轮毂高度平均风速，Hs为不规则波浪有义波高，Tp为不规则波浪谱峰周期，Sd/D为归一化平衡冲刷深度，L/D为单桩长径比。 |
3 冲刷平衡后大直径单桩水平受荷机理
Sd是指特定水流条件下冲刷达到稳定状态时对应的深度。虽然冲刷具有时间尺度，但决定结构破坏程度的往往是最不利的情况，因此本研究以整体平衡冲刷深度和局部平衡冲刷深度为变量，分别计算了风机以额定风速正常运行时，单桩基础埋深L=4.5D无冲刷(NS-1工况)、整体平衡冲刷深度Sgd=1.0D(GS-1工况)、局部平衡冲刷深度Sld=1.0D(LS-2工况)时的y/D和H/γ′D3-y/D，结果如图 7所示。由图 7a可知，相同冲刷深度情况下，整体冲刷平衡后桩基泥面处水平位移响应峰值是局部冲刷情况下的2倍，是无冲刷情况下的2.8倍。水平土反力-位移响应曲线的斜率表示桩-土相互作用的刚度，由图 7b可知，与无冲刷相比，整体冲刷导致的桩-土相互作用水平刚度弱化程度和桩身泥面处变形累积速率，明显高于局部冲刷。
为进一步揭示2种冲刷形式导致桩-土相互作用水平刚度弱化的本质差异，绘制了GS-1和LS-2工况中海上风电单桩基础受风浪联合作用600 s后桩周土体位移矢量图，如图 8所示。
由图 8a可知，对于水平循环受荷的大直径刚性单桩基础，当桩周发生整体冲刷时，桩身变形模式与未冲刷的平坦海床一致，即桩身几乎不产生变形，而是绕某深度处的一点发生刚性旋转，此时桩周土体的破坏模式为浅层区域发生楔形破坏，深层区域发生旋转剪切破坏，并且桩底表现出明显的“反向踢脚”现象。由于循环加载过程中，浅层桩周土体主、被动楔形破坏模式交替变换，因此会形成一个近乎对称的楔形土体破坏区域，并且随着加载的持续进行，楔形破坏区域将不再扩展，仅桩周土体累积变形持续增加。因此，整体冲刷地形中，大直径单桩水平土体抗力主要由上部楔形区和基底旋转区的水平土阻力提供。
由图 8b可知，局部冲刷地形的不同之处在于，每个往复加载过程中，位于主动楔形破坏区域冲刷坑内的土体会发生局部滑移，最终稳定的浅层楔形破坏区域较冲刷坑范围更大，而冲刷坑的深度和横向扩展范围不仅决定着单桩浅埋区的土体约束范围，还影响着剩余土体固结状态。因此，局部冲刷地形中单桩基础水平的土抗力不仅取决于上部楔形区和基底旋转区土阻力的贡献，还须考虑局部冲刷坑的影响。
4 整体冲刷对单桩动力响应的影响
4.1 考虑冲刷固结效应的单桩动力响应分析
对于正常固结海床地基中的单桩风电基础而言，冲刷不仅会减小单桩的埋置深度，冲刷诱发的表层土体运移还会导致剩余土体处于超固结状态。为定量评估整体冲刷导致的剩余土体固结状态改变对单桩水平承载力响应的影响，本研究分别对比了风机以额定风速正常运行时，单桩基础在L=4.5D无冲刷(NS-1工况)、L=3.5D无冲刷(NS-2工况)和L=4.5D且Sgd=1.0D(GS-1工况)时的泥面处水平位移响应和水平土反力-位移响应以及桩-土接触面上土压力，结果如图 9所示。
由图可知，整体冲刷会导致桩身位移响应显著增大，桩-土相互作用水平刚度大幅降低，对于NS-2和GS-1工况，虽然循环加载过程中桩基础埋深一致，但前者为正常固结土体，后者经历卸载过程，为超固结土体。对比图 9a与9b可知，考虑土体固结状态的改变会使单桩泥面处位移响应峰值降低约23%，桩身变形累积速率也明显降低，并且剩余土体由正常固结状态变化为超固结状态，土体将表现出更高的抗剪强度和抗压缩性能，因此相同加载强度下，正常固结土体中的单桩所受的土压力要明显高于超固结土体，如图 9c所示。上述结果也进一步说明，采用冲刷前的土体参数确定的桩-土相互作用刚度评估冲刷后单桩水平位移响应，会使位移峰值高估约23%。
4.2 考虑冲刷固结效应的桩-土相互作用分析方法
其中γw为水的重度。
φ′可表示为[28]
冲刷平衡后砂土的相对密度可由砂土的加卸载曲线推导得到，详细的推导过程可参考文[12]，由此得到的冲刷后砂土相对密度Drsc表达式为
其中：κ为砂土卸载系数，OCR为土体固结系数。
对于海上风机单桩基础整体冲刷，设计过程中可将原始海床面降低一个整体冲刷深度，并根据上述考虑冲刷固结效应的土体参数由API规范[5]确定修正海床面以下的p-y曲线和t-z曲线，然后由Winkle地基梁模型计算桩-土相互作用。
5 局部冲刷对单桩动力响应的影响
其中：Wh为楔形破坏区域高度；Wh0为无冲刷时的归一化楔形破坏区域高度；h为局部冲刷深度影响因子；Wb为楔形破坏区域宽度；Wh0为无冲刷时的归一化楔形破坏区域宽度；b为局部冲刷扩展范围影响因子。
上述结果进一步说明，局部冲刷引起的桩-土相互作用刚度弱化程度主要取决于楔形破坏区域贡献的土阻力的大小。虽然楔形破坏区域会随Sld的增加而扩大，但冲刷深度内土体约束范围及上覆土压力会显著减小，总体表现为水平土阻力的降低，桩-土相对刚度下降。
6 结论
为探究整体冲刷和局部冲刷对海上风电单桩基础水平承载力响应的影响，本文建立了冲刷平衡后的桩-土相互作用有限元模型开展分析计算，主要结论如下：
1) 整体冲刷地形中大直径单桩水平土抗力主要由上部楔形区和基底旋转区的土阻力提供，局部冲刷深度及其横向扩展范围主要通过影响桩周楔形破坏区的桩-土相互作用，进而影响土阻力的大小。
2) 冲刷诱发土体固结状态改变会提升剩余土体抗剪强度和抗压缩性能，进一步考虑冲刷固结效应会使冲刷后桩-土相互作用刚度弱化程度降低，最终表现为单桩泥面处位移响应峰值降低约23%，因此采用冲刷前的土体参数确定桩-土相互作用刚度评估冲刷后单桩水平承载力响应，结果偏大。
本研究采用的有限元方法能够准确量化冲刷效应对循环受荷大直径单桩动力响应的影响，但该方法涉及复杂的有限元二次开发，不利于工程设计普及应用，后续将进一步深入分析，建立能够考虑冲刷效应的桩-土相互作用简化分析模型。
