海上风电单桩基础-土相互作用特性影响因素分析

2021-02-02孔德森邓美旭

海洋工程 2021年1期

孔德森，刘一，邓美旭，侯迪

(1. 山东科技大学土木工程与建筑学院，山东青岛 266590； 2. 山东科技大学山东省土木工程防灾减灾重点实验室，山东青岛 266590； 3. 中铁建工集团山东有限公司，山东青岛 266590)

海上风电产业在全球的推广范围正在不断加大，对适应全球气候环境的变化和加快低碳能源的转型具有极其重要的意义，目前全球已有50多个国家和地区开始发展海上风电[1]。我国海上风电投资建设速度、安装、工程技术能力均不断获得突破，海上风电逐步走向成熟，风力发电成本将持续降低，是风电产业未来发展的重心[2]。大直径单桩基础凭借其加工制造简便、安装容易、结构受力明确等优势，被广泛应用于近海风力发电工程项目中[3]。大直径单桩基础在服役期间内受到来自风、洋流、波浪等水平荷载的作用，其水平循环受荷特性显著[4]。

国内外学者对桩-土相互作用进行了相关研究，并且取得了一定的成果。Liang等[5]在广义BNWF法中引入逆静力p-y滞回环，分析了循环荷载作用下桩-土相互作用。孙永鑫等[6]通过室内模型试验，研究了水平循环作用下粉土地基中近海风机刚性单桩的桩土相互作用规律。lamo等[7]研究了基于单桩的土-结构相互作用对海上风机动力特性的影响。陈仁朋等[8]进行了饱和粉土地基中单桩和群桩在水平循环荷载作用下的模型试验，得到了单桩和群桩随循环加载的受力变形规律。Arshad等[9]研究了考虑桩土相互作用的桩土应变累积模型，并给出一种具有单桩基础系统的OWT支撑结构的设计实例。隋倜倜等[10]建立了波浪作用下三维单桩-海床动力响应模型，通过定量分析超孔隙水压力和土体初始有效应力的变化，讨论了单桩插入深度对海床液化的影响机制。付鹏等[11]基于流固耦合建立了海上风机桩-土三维数值计算模型，研究了不同波浪荷载耦合作用下近海风电单桩基础与土的相互作用特性。刘修成等[12]以马尔代夫中马友谊大桥为背景，对主墩大直径钢护筒沉桩过程进行了监测，研究了珊瑚礁地质大直径打入桩的承载性能。但是，国内外学者对于水平循环荷载作用下非均质土中海上风电单桩基础-土相互作用的研究比较少。

以海上风电单桩基础为研究对象，采用有限元分析软件ABAQUS建立了非均质土中海上风电单桩基础-土相互作用数值计算模型进行研究，在模型中将波浪、洋流、风荷载等效成双向对称循环荷载，研究了非均质土中海上风电单桩在水平循环荷载作用下不同影响因素对桩身水平位移、剪力和弯矩的影响规律，取得了一些有意义的研究成果。

1 桩-土相互作用数值计算模型的建立

以我国东海某近海海域风电场为例，建立了非均质土中海上风电单桩基础-土相互作用数值计算模型，具体包括桩与土体相关参数的确定、土体本构模型选取、边界条件的设置、网格的划分、荷载的确定和施加、初始地应力的平衡等。

1.1 模型概况

由于荷载和结构的对称性，选择半个物理模型的桩土体系进行建模，数值模型为半圆柱体。为了忽略边界条件对桩土体系的影响，在几何模型上，用大尺寸来模拟半无限空间体，土体直径取20D(D为桩径)，土体高度取2hem(hem为桩的嵌固深度)，该尺度可满足最小边界尺寸和计算精度的要求[13]。

图1 桩-土相互作用数值计算模型

ABAQUS中通过桩-土表面定义接触属性以模拟桩与土之间的剪力传递和相对位移，采用主-从接触算法，选择刚度大的桩体为主控面，土体表面为从属面，桩-土法向行为采用硬接触，切向行为采用摩尔-库伦摩擦罚函数形式，界面滑动摩擦系数选取u=tan(0.75φ)(φ为土体内摩擦角)[14]，接触对采用面对面接触与有限滑移。边界条件是约束断面处y方向位移，约束模型侧面x和y方向位移，模型底端为固定约束。桩体和土体都采用8节点6面体线性减缩积分三维实体单元(C3D8R)。为了减小计算误差，同时也为了缩短计算时间，在桩土接触面附近单元网格划分得较细，而在远离接触面的土体，网格划分的相对稀疏。在实际状况中，假定海床泥面处位移为零，土体内部是有应力存在的，因此在施加水平荷载前必须进行初始地应力平衡[15]。采用ODB导入法进行初始地应力的平衡。桩-土相互作用数值计算模型如图1所示。

为了模拟桩-土之间的非线性，根据土层性质采用基于Mohr-Coulomb破坏准侧的理想弹塑性本构模型来模拟，孔位学等[16]针对非关联流动法则下岩土材料的剪胀角选取进行探讨，认为在非关联流动法则条件下采用剪胀角ψ=φ/2所得到的滑移线场与Prandtl理论一致。桩周土及桩端土参数列于表1。桩体采用线弹性模型来模拟，桩的物理力学参数列于表2，其中15 m位于水中，1 m位于水面之上。

表1 桩周土及桩端土参数

表2 桩的物理力学参数

1.2 土体本构模型的选取

Mohr-Coulomb本构模型不仅可以反映土体的抗压强度不同的S-D效应及其对静水压力的敏感性，而且实用简单，土体的黏聚力和内摩擦角比较容易获得，因而得到广泛的应用[17]。

ABAQUS中采用连续光滑的椭圆函数来作为塑性势面，其表达式为：

(1)

式中：ψ为剪胀角；c|0为初始黏聚力，即没有发生塑性变形时的黏聚力；ε为子午面上的偏心率，它用来控制G在子午面上形状与函数渐近线之间的相似度。Rmw则控制了其在π面上的形状，其表达式为：

(2)

式中：Θ为极偏角；Rmc为偏应力系数；e是π面上的偏心率，主要控制了π面上Θ=0～π/3的塑性势面的形状。默认值可根据下式计算：

(3)

由上式计算的e可确保塑性势面在π面拉压角点处与屈服面相切。

1.3 荷载的确定与施加

海上风机单桩基础不可避免的受到波浪、洋流、风等水平循环荷载的作用，这些荷载对桩体产生正向加载-正向卸载-反向加载-反向卸载的不断循环的过程。为了建模方便和获得较强的规律性，将这些荷载等效成双向对称循环荷载[18]的形式来模拟海上风机单桩基础-土相互作用特性，模型中荷载沿x轴方向不断循环加载，在x轴上方时荷载为正，在x轴下方时荷载为负。ABAQUS中采用周期型幅值曲线来定义水平循环荷载[19]，周期型幅值曲线用傅里叶(Fourier)级数表示。

t≥t0时，幅值表达式为：

(4)

a=A0

(5)

式中：N为傅里叶级数向的个数；ω为圆频率，其值为ω=2πf，f为频率；t0为起始时刻；A0为初始幅值；An和Bn为系数。

1.4 接触面的设置

图2 主控面光滑处理示意

如果采用点对面离散方式，为了减少主控面上锯齿状节点穿透从属面而影响从属面上节点滑动的现象，ABAQUS通常对主控面进行光滑处理，如图2所示。光滑处理在有限滑动、点对面的离散分析中尤为重要，否则会引起主控面的法线方向出现不连续的变化，从而出现收敛问题。在ABAQUS中光滑化的程度是通过系数K=m1/l1来控制的，默认值为0.2，其值不能超过0.5。如果采用面对面离散方式，ABAQUS对主控面不会进行自动光滑化处理。由于面对面的离散方式是基于类似的平均意义而建立接触条件的，故在某种程度上认为其是内在的光滑化。接触模拟过程中应该确保从属面位于主控面法线方向所指的一侧，否则计算不能够收敛。在面对面离散中，如果主控面和从属面的法线方向相同，将不会考虑接触[20]。

2 数值计算结果分析

在模型中定义参考点，将参考点与水面处桩身横截面建立分布耦合约束，在参考点上施加水平循环荷载。为了节约计算时间，且不考虑循环次数的影响，故对第20次循环时不同因素对桩身水平位移、剪力和弯矩的影响规律进行了分析。

2.1 水平极限承载力的确定

图3 桩顶荷载位移曲线

由于海上风电单桩基础桩径较大且采用钢管桩，桩身强度非常大[21]，且桩的水平荷载-位移曲线为缓变型，该曲线没有明显的拐点，可以认为是渐进式破坏，所以海上风电单桩基础的水平极限承载力主要由桩体的水平变形控制。采用位移控制法，对桩顶施加0.3 m的水平位移，在有限元软件ABAQUS后处理中提取水平支反力和桩身水平位移的相关数据，然后利用绘图软件Origin绘制水平支反力与水平位移之间的关系曲线，得到桩顶荷载位移曲线，如图3所示。

根据允许变形法[22]，将桩顶水平位移达到0.02D时对应的水平荷载，确定为海上风电单桩基础的水平极限承载力。根据荷载位移曲线，当位移为0.1 m时对应荷载1.31 MN。为了反映单桩基础所受的荷载水平，单桩上施加的水平循环荷载幅值的大小为其水平极限荷载Fu的一定比值，荷载频率根据工程资料取值。

2.2 循环荷载比的影响

取循环荷载比为0.2、0.4、0.6、0.8、1.0，荷载频率为0.1 Hz，对不同循环荷载比下桩身水平位移、剪力和弯矩沿埋深的变化规律进行了研究。

2.2.1 桩身水平位移分析

由于第20次循环过程中193 s和198 s时桩身水平位移值最大，且193 s时和198 s时桩身水平位移沿埋深曲线关于y轴大致呈对称分布，故只对193 s时不同循环荷载比下桩身水平位移沿埋深的变化规律进行了分析，其沿埋深分布曲线如图4所示。

由图4可知，由于桩侧土抗力的抵抗作用，不同循环荷载比下桩身水平位移在193 s时沿埋深均出现零点，桩身正向位移最大值均出现在泥面处，负向最大位移均出现在桩底端，桩身的变形主要集中于埋深30 m以内，说明此范围内桩身挠曲明显，稳定性较差。桩身位移绝对值随荷载幅值的增大而增加，且增加的幅度有所增大，这是因为不同循环荷载比作用下桩周土体塑性变形逐渐累积。随着荷载的增加，桩身位移零点埋深分别为33.31 m、33.59 m、33.83 m、33.86 m、35.23 m，桩身位移零点沿埋深逐渐下移，究其原因，当水平荷载较小时，主要是浅层土体提供土抗力，随着水平荷载的逐渐增加，浅层土体开始塑性屈服，深层土体的土抗力逐渐发挥，导致位移零点逐渐下移[23]。由零点位置可知，循环荷载比为0.2～0.8时，零点下移不明显，循环荷载比为1.0时，位移零点下移显著，说明水平荷载接近正常使用状况下极限荷载时，浅层土体与之前相比产生较大的塑性区。

图4 193 s时不同循环荷载比下桩身水平位移变化曲线

图5 193 s时不同循环荷载比下桩身剪力变化曲线

2.2.2 桩身剪力分析

由于第20次循环过程中193 s时桩身正向剪力值最大，故对193 s时不同循环荷载比下桩身剪力变化规律进行分析，其沿埋深分布曲线如图5所示。

图6 193 s时不同循环荷载比下桩身弯矩变化曲线

由图5可知，193 s时不同循环荷载比下桩身剪力沿埋深均出现反弯点，这是由于桩侧土体的抵抗作用引起的。反弯点以上桩身剪力值均为正且随埋深增加逐渐减小，泥面处桩身剪力均达到最大值，反弯点以下剪力值均为负，且剪力绝对值随埋深的增加先增大后减小，均在埋深31～32 m范围内达到负方向最大值，桩底侧剪力值比较小。随着循环荷载比的增加，剪力零点以上范围内剪力减小幅度变大，究其原因，桩身位移随循环荷载比的增加而变大，引起桩周土体压缩变范围变大，能够提供较大的土抗力，故剪力变化幅度较大。随着循环荷载比的增加，桩身剪力反弯点逐渐下移，这主要是由于随着荷载的增加，浅层土体发生屈服，桩侧土抗力的发挥沿埋深有所下移。

2.2.3 桩身弯矩分析

由于第20次循环过程中193 s时桩身正向弯矩值最大，故对193 s时不同循环荷载比下桩身弯矩变化规律进行分析，其沿埋深分布曲线如图6所示。

由图6可知，193 s时随着循环荷载比的增加，桩身弯矩均没有出现反弯点，弯矩均沿埋深先增加后减小，桩身最大弯矩主要位于距泥面7～9 m范围内，即弯矩最大值发生在浅层土体。桩身弯矩值随着循环荷载比的增加逐渐增大，且增大的程度大致相同，这是因为循环荷载比增加，桩身承担的土抗力变大，故桩身承担的弯矩增加。

2.3 荷载频率的影响

取循环荷载比为0.6Fu，荷载频率为0.05 Hz、0.10 Hz、0.20 Hz、0.25 Hz、0.40 Hz，研究了不同荷载频率下桩身水平位移、剪力和弯矩沿埋深的变化规律。

图7 第3 s时不同荷载频率下桩身水平位移变化曲线

2.3.1 桩身水平位移分析

对第3 s时不同荷载频率下桩身水平位移变化规律进行了分析，其沿埋深分布曲线如图7所示。

由图7可知，不同荷载频率下第3 s时桩身位移均出现零点，零点以上桩身位移均沿埋深逐渐减小，桩身水平位移随荷载频率增加而不断增大。零点以下桩身位移绝对值均沿埋深逐渐增加，不同荷载频率时桩身位移曲线几乎重合。随着荷载频率的增加，桩身位移零点沿埋深逐渐下移，究其原因，荷载频率增大，浅层土体循环弱化明显，强度降低，深层土体土抗力开始发挥，导致桩身位移零点下移。

2.3.2 桩身剪力分析

对第3 s时不同荷载频率下桩身剪力变化规律进行了分析，其沿埋深分布曲线如图8所示。由图8可知，不同荷载频率下第3 s时桩身剪力均出现反弯点，主要位于埋深8～9 m范围内。反弯点以上不同荷载频率时桩身剪力曲线沿埋深几乎重合，反弯点以下桩身剪力曲线在埋深约28 m处出现分界点，分界点以上桩身剪力绝对值随着荷载频率的增加而减小，且变化不大，分界点以下桩身剪力绝对值随着荷载频率的增加而变大，桩身剪力变化相对其它位置较大，主要集中在埋深32～44 m范围内。桩身剪力负方向最大值的绝对值随荷载频率的增加有所增大。

图8 第3 s时不同荷载频率下桩身剪力变化曲线

图9 第3 s时不同荷载频率下桩身弯矩变化曲线

2.3.3 桩身弯矩分析

对第3 s时不同荷载频率下桩身弯矩变化规律进行分析，其沿埋深分布曲线如图9所示。由图9可知，不同荷载频率下第3 s时桩身弯矩均沿埋深先增加后减小，均没有出现反弯点，随着荷载频率的增加逐渐增大，且变化不明显。桩身弯矩最大值主要分布在埋深8～9 m范围内，桩底端弯矩均比较小，接近于0。

2.4 加载方式的影响

海上风电单桩基础所受的荷载条件比较复杂，荷载形式并不唯一，因此开展了加载方式对海上风电单桩基础-土相互作用特性的影响规律研究。不同荷载施加方式如图10所示。

图10 不同荷载施加方式示意

取荷载幅值为0.6Fu，荷载频率为0.1 Hz，对不同加载方式下桩身水平位移、剪力和弯矩沿埋深的变化规律进行了研究。

2.4.1 桩身水平位移分析

开展了第20次循环正向荷载和负向荷载作用结束时不同加载方式下桩身水平位移沿埋深变化规律研究，其分布曲线如图11所示。

图11 不同加载方式下桩身水平位移变化曲线

由图11(a)可知，单向循环荷载作用下桩身水平位移最大，双向对称循环荷载作用下桩身水平位移最小，究其原因，桩身水平位移在循环荷载作用下发生累积，而单向循环荷载作用时荷载方向始终为正，桩身位移累积程度最大，双向对称循环荷载作用时，正向荷载和负向荷载作用时长相等，桩身正向累积位移和负向累积位移存在相互抵消，故桩身水平位移最小。不同加载方式时桩身位移均出现零点，零点以上桩身位移均为正，沿埋深逐渐减小，零点以下桩身位移为负，绝对值均沿埋深逐渐增加。随着桩身位移的增加，桩身水平位移零点位置沿埋深下移，究其原因，桩身位移越大，上部土体变形越大，开始发生屈服，深层土体土抗力开始发挥，导致桩身水平位移零点下移。由图11(b)可知，单向循环荷载作用下桩身水平位移显著大于另外两种循环荷载作用下的桩身水平位移，这主要是由于，双向循环荷载作用时桩身正向累积位移和负向累积位移存在相互抵消。

2.4.2 桩身剪力分析

开展了第20次循环正向荷载和负向荷载作用结束时不同加载方式下桩身剪力沿埋深变化规律研究，其分布曲线如图12所示。

图12 不同加载方式下桩身水平位移变化曲线

由图12(a)可知，正向荷载作用结束时，双向不对称循环荷载作用下桩身剪力最大，双向对称循环荷载作用下桩身剪力最小。不同加载方式下桩身剪力均出现反弯点，反弯点以上桩身剪力值均为正，沿埋深逐渐减小，反弯点以下桩身剪力值均为负，剪力绝对值沿埋深先增大后减小。不同加载方式时，桩身剪力均沿负方向出现最大值，主要集中在埋深32～34 m范围内。桩身剪力反弯点位置随着桩身剪力的增大而沿埋深有所下降。

由图12(b)可知，负向荷载作用结束时，单向循环荷载作用下桩身剪力绝对值最大，双向不对称循环荷载作用下桩身剪力绝对值最小。负向荷载作用结束时不同加载方式下桩身剪力变化规律与正向荷载作用结束时不一致，这是由于循环过程中桩身剪力发生累积引起的。

2.4.3 桩身弯矩分析

开展了第20次循环正向荷载和负向荷载作用结束时不同加载方式下桩身弯矩沿埋深变化规律研究，其分布曲线如图13所示。

由图13(a)可知，正向荷载作用结束时，不同加载方式下桩身弯矩均没有出现反弯点，均沿埋深先增大后减小，双向不对称循环荷载作用下桩身弯矩最大，双向对称循环荷载作用下桩身弯矩最小，桩身最大弯矩位于埋深8～9 m范围内，桩底端弯矩比较小，且相差不大。

由图13(b)可知，负向荷载作用结束时，双向对称循环荷载下桩身弯矩出现反弯点，其余荷载作用下桩身弯矩没有出现反弯点，单向循环荷载作用下桩身弯矩绝对值最大，双向不对称循环荷载作用下桩身弯矩绝对值最小。负向荷载作用结束时桩身弯矩变化规律与正向荷载作用结束时不一致，这是由于循环过程中桩身弯矩发生累积引起的。

图13 不同加载方式下桩身水平位移变化曲线

图14 193 s时不同桩壁厚下桩身水平位移变化曲线

2.5 桩壁厚的影响

取荷载幅值为0.6Fu，荷载频率为0.1 Hz，对不同桩壁厚时桩身水平位移、剪力和弯矩沿埋深的变化规律进行了研究。

2.5.1 桩身水平位移分析

193 s时不同桩壁厚时桩身水平位移沿埋深分布曲线如图14所示。由图14可知，不同桩壁厚下桩身位移均出现零点，零点之上桩身位移曲线出现分界点，大约位于埋深12 m处，说明不同桩壁厚下桩身发生不同程度的屈曲。分界点以上桩身位移基本随桩壁厚增加而减小，而分界点以下正好相反，为了减小桩身位移，此范围可以考虑设计“上厚下薄”的钢管桩。零点之下桩身位移曲线也出现交点，大约位于距桩端4 m处。随着壁厚的增加，桩身位移零点随桩壁厚的增加逐渐下移，究其原因，桩身刚度随壁厚增加有所增大，导致桩身承担的土抗力值增加，进而引起桩身位移零点下移。从图14中还可以看出，随着桩壁厚的增加，桩身水平位移变化程度有所降低，说明桩壁厚较小时对桩身沿埋深水平位移影响较大。

2.5.2 桩身剪力分析

对193 s时不同桩壁厚时桩身剪力变化规律进行分析，其沿埋深分布曲线如图15所示。

由图15可知，不同桩壁厚下桩身剪力沿埋深均出现反弯点，反弯点以上桩身剪力曲线几乎重合，说明此范围内桩壁厚对桩身剪力影响不大。反弯点以下桩身剪力曲线在埋深约30 m处出现交点，交点以上桩身剪力绝对值随桩壁厚的增加而有所减小，而交点以下桩身剪力绝对值随桩壁厚的增加而增大。究其原因，桩身剪力大小与桩侧土抗力有关，交点以上土抗力随桩壁厚增加有所减小，而交点以下土抗力随桩壁厚增加有所增大。随着桩壁厚的增加，桩身剪力反弯点沿埋深逐渐下移。桩身剪力正方向最大值均位于泥面处，负方向最大值随着桩壁厚的增加绝对值逐渐增大，主要位于距泥面29～32 m范围内。

图15 193 s时不同桩壁厚下桩身剪力变化曲线

图16 193 s时不同桩壁厚下桩身弯矩变化曲线

2.5.3 桩身弯矩分析

对193 s时不同桩壁厚时桩身弯矩变化规律进行分析，其沿埋深分布曲线如图16所示。由图16可知，193 s时不同桩壁厚下桩身弯矩均没有出现反弯点，沿埋深先增加后减小，基本随桩壁厚的增加而有所增大，这是因为桩壁厚越大，桩身刚度越大，承担的土抗力增加，桩身弯矩从而变大。仅在埋深8 m和下部砂土层中变化不明显，此范围内桩壁厚对桩身弯矩影响不大。

3 数值计算结果对比分析

由于模型试验中桩的直径相对较小，无法直接与数值模型中的大直径桩进行对比验证，故与已有的有限元分析结果进行对比分析。陈新奎[24]采用ABAQUS对大直径桩的水平受荷特性进行了模拟，并将数值模拟结果与风电场现场试桩结果进行了对比分析，验证了ABAQUS数值模拟的可靠性。取文献[24]中4 MN (此荷载与本文荷载幅值最接近)水平力作用下桩身水平位移、剪力、弯矩与193 s时桩身水平位移、剪力、弯矩(此时桩身水平位移、剪力、弯矩最大)进行对比分析，如图17、图18和图19所示。

图17 桩身水平位移对比分析

图18 桩身剪力对比分析

从图17～图19可以看出，泥面以下范围内桩身水平位移、剪力、弯矩沿埋深变化规律类似。由图17可知，桩身水平位移曲线沿埋深均出现零点，均呈现出非线性的变化，桩均绕桩身轴线上某一点转动，表现出刚性桩的性质。由图18可知，泥面以下范围内桩身剪力均出现反弯点，图18(a)中桩身剪力变化幅度较大，究其原因，桩身水平位移较大，引起桩周土体压缩范围较大，能够提供较大的土抗力，导致剪力变化程度较大[25]。由图19可知，泥面以下范围内桩身弯矩沿埋深均没有出现反弯点，均沿埋深先变大后变小，最大值均发生在浅层土体，桩底端弯矩均接近于零。

综上可知，运用ABAQUS建立的数值模型是可靠的，且数值模拟结果是准确的。