长期循环荷载下海上单桩基础变形特性的数值模拟研究

杨钰荣，孙毅龙，许成顺，翟恩地

（1.北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124；2.温州大学建筑工程学院，浙江温州 325035；3.金风科技股份有限公司，北京 100176）

我国近海岸以大规模的砂质沉积体为主，大直径单桩基础是海上风电系统常用的基础形式之一。海上风电单桩基础在服役期间，长期承受风荷载、浪荷载等引起的循环倾覆力矩，荷载循环作用次数可达到108次[1]。单桩基础长期承受这种低频的循环荷载作用，会导致桩周土体刚度发生变化，造成桩基侧向累积变形增大，严重危及海上风电系统的正常运行。因此，海上风电单桩基础在长期荷载作用下的累积变形是桩基设计中重点考虑的因素之一。

目前，海上风电单桩基础累积变形特性的研究主要采用修正p-y曲线法、室内或现场试验和数值模拟分析方法等。MCCLELLAND B等[2]于1958年基于温克尔地基梁模型提出了计算水平荷载作用下单桩非线性变形的方法，即p-y曲线法。2006年，LESNY K等[3]发现美国石油学会（American Petroleum Institute，API）推荐的p-y曲线法高估了地基初始刚度，针对该问题对p-y曲线法进行了简单的修正。HYUNSUNG L等[4]于2018年发现荷载作用后地基刚度变化影响风机系统频率，在不同的加载频率下，利用模型试验建立拟静力分析的简化动力p-y曲线。2020年，胡安峰等[5]基于刚度衰减模型对p-y曲线进行修正，综合考虑了荷载特性、循环次数、砂土刚度循环弱化等因素。2021年，王卫等[6]针对p-y曲线过于保守等问题，提出一种综合考虑桩径和地层深度影响的黏土p-y修正模型。尽管修正p-y曲线法在一定程度上考虑了土体条件、桩基特点和荷载环境对预测桩基累积变形的影响，但考虑荷载的循环作用次数与实际桩基受荷次数相差较大，对于预测大直径桩基受长期荷载作用的循环累积变形相对保守。

为了进一步研究海上风电单桩的循环累积变形问题，不少学者基于室内外模型试验建立循环累积变形与静态荷载作用下的变形、循环次数的经验函数关系，提出不同类型的水平循环受荷桩的累积变形预测模型。VERDURE L等[7]、BIENEN B等[8]、LI W等[9]基于水平循环加载模型试验研究了长期荷载作用下单桩的动力响应特性，认为桩基的循环累积变形与桩基初始变形呈对数关系。LEBLANC C等[10]、RASMUS等[11]、吴金标等[12]通过开展离心机模型循环加载试验，认为幂函数模型能够较好地预测砂土场地中水平受荷桩的循环累积变形。

由于试验条件、荷载特性等差异，上述学者所建立的经验预测模型不尽相同，对数形式和幂函数形式的模型预测结果存在差异。为此，众多学者基于循环试验，通过建立桩—土相互作用数值模拟，研究长期荷载作用下海上风机单桩基础的累积变形特性。ACHMUS M等[13]于2008年基于循环三轴试验结果提出长期荷载作用下的刚度衰减模型，利用该模型预测N次循环荷载作用下单桩累积侧向变形，并与HETTLER A等[14]、LITTLE R L等[15]分别提出的对数函数型和幂函数型桩基累积变形预测模型进行对比；2015年，DEPINA I等[16]将三维有限元桩土模型与刚度退化模型耦合，模拟了砂土中单桩基础在长期循环荷载作用下的变形、转角和弯矩响应；罗如平等[17]、董爱民等[18]分别于2016年、2017年利用刚度衰减模型，分析长期荷载作用下大直径单桩的循环累积变形特性，并对现有的桩基累积变形预测模型进行修正；2018年，YANG M[19]基于Martin Achmus提出的刚度退化模型运用三维有限差分模型，研究了单桩基础的长期循环特性，定量地分析了循环荷载特性、桩基埋置深度、荷载偏心率对单桩长期累积变形的影响；2020年，LUO R等[20]基于刚度衰减模型建立了修正的桩基累积变形预测方法，开展不同桩径、荷载等对桩基累积变形的影响规律研究。

综上所述，刚度衰减模型是分析单桩基础累积侧向变形较为简便的方法。但目前该模型参数的确定尚存在困难。为此，本文利用竖向—扭转双向耦合循环剪切试验（模拟波浪等循环荷载作用）结果获得刚度衰减模型，并基于Flac3D有限差分平台开展三维桩—土相互作用数值模拟，分析循环荷载幅值、场地密实度和桩径等对长期荷载作用下单桩基础累积变形的影响规律；根据数值计算结果，对常用的幂函数型和对数函数型桩基累积变形预测模型的适用性进行评估，为海上风电桩基设计提供参考。

1 刚度衰减模型

1.1 刚度衰减模型介绍

HUURMAN M[21]基于室内循环三轴试验指出砂土在长期循环荷载下其模量会发生衰减，并将土体的累积塑性应变增量比改写为土体的割线模量增量比，具体表达式如式（1）所示，示意图如图1所示。

图1 砂土刚度衰减模型示意图

式中，EsN为循环荷载作用N次后的割线模量；Es1第一次循环荷载作用后的割线模量；εp，N=1为第一次循环荷载作用后的应变；εp，N为循环荷载作用N次后的应变。

ACHMUS M等[13]基于Huurman提出的土体割线模量衰减公式，建立了刚度衰减模型来表征不同循环次数下的桩基累积变形发展规律，具体公式如下。

式中，N为荷载循环次数；X为定义的循环应力比；b1、b2是根据试验得到的回归参数；C是由参数b1、b2和循环应力比X控制的指数。

由于在实际海上风电桩—土相互作用系统中，土体存在非均匀初始应力，即使未施加循环应力，也会存在初始的非零循环应力比，使式（3）得到一个不符合现实的应变增量，对地基模量计算存在影响。考虑非均匀初始应力的影响，定义特征循环应力比Xc，在数值模型计算中用以替换X。具体表达如下。

式中，X(1)为受荷状态下土体的循环应力比；X(0)为初始应力状态下土体循环应力比；σ1为土体最大主应力；σ3为土体最小主应力；σ1，sf为土体在静态力下破坏时的主应力；σ1，cyc为施加循环力时土体的最大主应力；φ为土体内摩擦角。

1.2 土体刚度衰减特性研究

本文利用竖向—扭转双向耦合剪切仪[22]开展相对密实度Dr=35%、Dr=50%、Dr=70%的饱和砂土（福建标准砂）循环排水剪切试验，得到不同相对密实度砂土在长期循环荷载作用下的模量衰减规律。为了模拟实际海床桩周土的受力状态，在竖向—扭转耦合循环剪切试验中，对试样同时施加轴向荷载及扭矩，并保持二者的相位差为90°；设置固结应力比Kc=2，试样初始固结围压σ1=100 kPa、σ2=σ3=50 kPa。海上风机桩基周围土体主要受风、波浪荷载的影响。常见波浪和风荷载的典型频率为0.1 Hz[23]。砂性土在这种低频的循环荷载作用下易排水。因此，循环试验中的加载频率为0.1 Hz。通过简单的数值计算，选取桩周土的循环应力比X为0.10、0.12、0.15，试验方案见表1。

表1 试验方案

利用式（1）计算每一次循环荷载作用后的割线模量，得到土体的割线模量随循环加载次数的变化规律，如图2所示。由图可知，割线模量在循环荷载作用前100个周期内衰减最为显著，此后逐渐衰减至某一数值，且与荷载幅值没有明显的依赖关系。

图2 不同相对密实度砂土的割线模量衰减规律

根据前文介绍的刚度衰减模型，基于式（4）对不同相对密实度砂土的割线模量进行拟合，得到割线模量对数比与循环加载次数对数之间的关系，如图3所示。由图可知，两者之间呈近似的线性关系，且不同相对密实度砂土的割线模量衰减曲线斜率C不同。此外，由式（5）可知，斜率C由循环荷载影响指数b1、b2和循环应力比X控制。

图3 不同相对密实度砂土的割线模量衰减比

为了将循环试验结果应用于桩基的循环累积变形分析，基于试验结果利用式（6）对参数b1、b2进行反算，得到b1、b2和土体相对密实度Dr之间的关系，如图4所示。由图可知，参数b1、b2与相对密实度Dr呈线性关系式，b1、b2与Dr拟合公式如下（R2>0.99）。

图4 参数b1、b2与砂土相对密实度关系曲线

1.3 刚度衰减模型介绍

应用FLAC3D有限差分平台内置的FISH语言，根据试验结果编制相应的刚度衰减程序，并建立桩—土相互作用整体数值分析模型。将循环荷载作用下土体模量的衰减规律嵌入的数值分析计算模型中，实现土体模量随应力状态及循环次数的更新。在整体数值分析模型中分3个步骤实现桩基累积变形的计算，示意图如图5所示。刚度衰减模型的参数b1、b2根据实际场地的土体密实度，按照式（11）、式（12）计算取值。根据图5中的步骤，更新土体模量后计算的桩基侧向变形即为不同循环次数荷载作用后桩基的循环累积变形。

图5 桩基累积变形计算步骤图

2 数值模型验证

基于上述刚度衰减数值分析模型，模拟ACHMUS M等[24]开展的离心机模型试验及张纪蒙等[25]开展的水平加载模型试验；并与当前两种主要形式的桩基累积变形预测模型计算结果进行对比，具体桩基累积变形预测模型如表2所示。

表2 桩基累积变形预测模型

2.1 数值模型的建立

数值模型中，土体采用摩尔—库伦本构模型，同时将式（9）土体模量与小主应力之间的关系引入到数值模型中，考虑土体弹性模量沿深度的非线性变化，土体模量非线性增长公式如下。

式中，σat为大气压力；σm为土体单元平均应力；κ与λ为土体刚度参数。

桩体采用弹性模型，桩基的建模采用等效刚度法。桩—土之间的相互作用采用接触单元来模拟，接触单元的法向和切向刚度的取值方法按照式（14）计算。

式中，K、G为接触面两侧材料体积模量和剪切模量；ΔZmin为接触面两侧单元最小尺寸。

按上述各式和模型尺寸计算得到kn=ks=1 012 Pa/m[26]，桩—土接触单元参数的粘聚力和内摩擦角取相邻土体0.3~0.5倍[27]。桩基和土体的具体参数见表3、表4。

表3 ACHMUS M等[24]模型试验参数

表4 张纪蒙等[25]模型试验参数

桩—土相互作用数值计算模型水平方向上的土体长度取为12倍桩径，桩端底部的土体模型长度取为10倍桩径；保证在该模型尺寸下桩基性能不受边界的影响；为提高计算效率，在桩基附近划分较密的网格，外侧网格相对疏松；以ACHMUS M等[24]的桩基试验数值模型为例，如图6所示。

图6 ACHMUS M等[24]桩基试验数值计算模型

2.2 桩基累积变形计算对比

图7、图8分别为ACHMUS M等[24]、张纪蒙等[25]模型试验中的桩基累积变形计算结果。本文数值模型计算的桩基累积位移与模型试验实测值基本吻合，较好地反映了单桩基础长期受荷载作用的累积特性，可以用来分析长期循环荷载作用下海上风机单桩基础的循环累积变形。

图7 ACHMUS M等[24]模型试验对比

图8 张纪蒙等[25]模型试验对比

此外，两种预测模型计算的桩基累积变形与试验实测值的趋势基本一致；幂函数模型的预测结果与实测值更加接近；而对数型模型的预测结果在循环次数较大时，与试验实测值有一定偏差。

3 桩基累积变形特性分析

3.1 数值模型建立

基于前文建立的数值模型，考虑中密、密实、紧密砂土3类场地中，钢管桩桩长L=30 m，桩径D=4 m、5 m、6 m、7 m、8 m的单桩基础，桩基受水平侧向荷载与泥面高1 m，数值模型如图9，具体的桩、土参数如表5、表6所示。

图9 单桩基础三维数值模型

表5 数值模型土参数[28]

表6 数值模型桩参数

3.2 桩基水平极限承载力

为了确定数值模型中对桩基施加的荷载大小，先计算静力荷载下桩基水平极限承载力，通过引入无量纲荷载幅值系数φb以确定不同工况下的水平荷载大小。荷载幅值系数φb的表达式如式（15）所示。

式中，Hmax为施加的水平静力侧向荷载；Hu为计算的桩基水平极限承载力。

桩基水平极限承载力不是桩基破坏时的承载力，而是被定义为发生某一变形值时的水平荷载。AHMED S S等[29]、BARARI A等[30]提出单桩泥面水平位移为0.1D时所受的水平力为极限承载力。本文采用该极限承载标准，同时计算了不同条件下单桩基础的水平极限承载力。数值计算结果如图10所示。

图10 不同桩径桩基水平极限承载力

当场地土体的密实程度相同时，桩基水平极限承载力随桩径增加而增大，且二者呈现近线性增长的关系；这是由于随着桩径的增大，土体的破坏模式逐渐向半刚性桩，刚性桩的破坏模式发展，即土体受力区范围增大，因而桩基的水平承载力增加。此外，当桩径相同时，桩基的水平极限承载力随着场地密实度的增加而增大；例如，当D=5 m时，密实砂土、紧密砂土场地中桩基水平极限承载力分别为中等密实砂土场地中的1.25倍、1.52倍；这是因为土体的强度、刚度与密实度紧密相关，砂土越密实其承载能力越强，因而桩基极限承载力越大。

3.3 荷载幅值的影响

以中等密实场地中D=4 m的工况为例，计算桩基不同荷载大小下桩基的累积变形。设置4组不同的荷载幅值，φb=0.10、0.20、0.25、0.30，最大循环次数N=106，该工况下的桩基水平极限承载Hu=26.8 MN。

以荷载幅值系数φb=0.20、0.25、0.30的3种工况为例，绘制桩身侧向变形随循环次数变化曲线，如图11所示。由图可知，当桩基受静态侧向力作用时（N=1），荷载幅值系数越大，桩基的初始侧向变形越大；随着循环次数的增加，桩基侧向变形不断增长，且荷载幅值系数越大桩基位移累积越快；此外，在不同大小的荷载长期作用下，桩基转动中心始终位于泥面下24.8 m处，约为0.7倍的桩基埋深。

图11 不同荷载幅值下桩身位移发展

为了对比不同条件下桩基累积变形随循环荷载作用次数增加的发展规律，评价桩基累积变形预测模型的可靠性，对桩基泥面位移进行无量纲处理。根据数值模型计算结果，得到桩基泥面位移累积率yN/y1与循环荷载次数N的关系。最后，对比分析幂函数型、对数函数型桩基变形预测模型和数值模型计算的桩基累积位移发展规律。不同循环荷载幅值下桩基泥面位移累积率的发展规律如图12所示。

图12 不同荷载幅值下桩基泥面位移发展

由图12可知，桩基累积变形随荷载作用次数先增加随后逐渐趋于稳定，在循环加载前期（N<105），幂函数模型的预测结果偏小，与数值计算结果的误差在18.88%以内，对数函数模型的预测结果与数值计算结果吻合较好，其最大误差在13.61%以内；随着循环加载次数不断增加，当N>105时，幂函数模型的计算结果与数值计算值最大误差仅为1.42%，其预测结果更加准确；而对数函数模型预测值的最大误差为4.68%；随着荷载幅值系数增大，幂函数模型预测结果的偏差逐渐减小，而对数函数模型的预测结果偏差逐渐增大。分析表明，从桩基累积变形的长期发展来看，在分析不同荷载大小作用下的桩基累积变形时，幂函数模型的预测结果更加准确。

根据式（16）、（17）反算得到不同荷载大小下的模型参数α与tb，分析幂函数和对数函数预测模型中参数α与tb受荷载幅值影响的规律。

式中，y1为桩基的初始侧向位移；yN为循环荷载作用N次后桩基的侧向位移。

图13绘制了α和tb与φb的关系曲线。由图可知，参数α和tb均与荷载幅值系数φb呈正比。此外，对数函数模型参数拟合曲线的斜率大于幂函数模型参数拟合曲线的斜率，这说明对数函数预测模型的准确性受荷载幅值的影响更大，与图12所得结论一致，预测不同荷载幅值下桩基的循环累积变形规律时，幂函数预测模型的计算结果较为准确。

图13 模型参数与荷载幅值系数的关系

3.4 场地密实度的影响

为研究场地土体密实度对桩基循环累积变形的影响，考虑一桩径D=5 m的单桩设置于中密、密实和紧密砂土3种场地中，在桩顶施加循环的水平侧向荷载，分析桩基长期荷载作用下桩基累积变形规律。桩基的循环累积变形计算属于疲劳承载状态，基于DNV（Det Norske Veritas）规范，桩基的循环荷载幅值取为水平极限承载力的25%，即φb=0.25[31-32]。水平荷载幅值大小如表7所示。

表7 不同砂土中桩基承受的荷载幅值

图14为不同密实度砂土中桩身变形发展。由图可知，3种场地中的桩基侧向变形均随荷载作用次数增加而增大；场地土体越密实，桩基受静态荷载（N=1）作用后的初始位移越小，受循环荷载作用后的位移累积越小。此外，3种密实度砂土中桩基转动中心均位于泥面以下24.8 m位置处，约为0.7倍的桩基埋深，这说明桩基转动中心受场地密实度的影响较小。

图14 不同密实场地中桩身位移发展

为分析不同密实场地中桩基累积变形预测模型的适用性，根据式（16）和式（17）反算得到不同密实度砂土中的模型参数α与tb，并对比中密、密实、紧密砂土中，两种预测模型计算的桩基泥面位移累积率，如图15所示。由图可知，3种密实度砂土中桩基的泥面位移累积率在循环加载初期（N<105）增加显著随后逐渐趋于稳定，紧密砂土中桩基的位移累积率最小。

图15 不同密实度下桩身泥面位移发展

两种模型的预测结果显示，当N<105时，幂函数模型的预测结果与数值计算结果的最大误差约为32.42%，对数函数预测结果与数值计算结果的最大误差约为27.47%。当N>105时，对数函数预测模型结果与数值计算结果的最大误差为13.86%；幂函数模型的预测结果与数值结果的最大误差为4.23%，其预测结果更为准确。由上述分析可知，在不同密实程度的砂土场地中，使用幂函数模型预测桩基累积变形的长期发展更具优势。

为了分析土体密实程度对幂函数模型参数α和对数函数模型参数tb的影响规律，图16绘制了α和tb与土体密实度Dr的关系曲线。由图可知，幂函数模型和对数函数模型中的参数α和tb均随场地土体密实度的增加而降低。这是因为场地土体越密实，桩基位移累积率越小。此外，α与Dr关系曲线的斜率接近为0，说明密函数预测模型受土体密实度的影响较小。

图16 模型参数与土体相对密实度的关系

3.5 桩径的影响

考虑桩径D=4 m、5 m、7 m、6 m、8 m的单桩设置于中密砂土场地中，计算不同桩径的桩基在长期荷载作用下的累积变形，不同桩径桩基所施加的荷载大小如表8所示。

表8 不同桩径桩基的荷载大小

图17以桩径D=4 m、6 m、8 m的桩基为例，绘制了长期荷载作用下桩身侧向变形的发展。由图可知，桩径越大，桩基的侧向变形越大；当桩基受静力荷载作用时，3种桩径桩基的泥面初始位移分别为5.3 cm、7.7 cm、9.1 cm；当N=106时，桩径D=4 m、6 m、8 m的桩基侧向变形分别增加了83.5%、203.3%、334%。此外，3种桩径桩基的转动中心分别位于泥面下23 m、24.5 m、25 m，桩基转动中心随桩径的增大而呈现下移趋势，但不随循环次数的增加而变化。

图17 不同桩径单桩基础的桩身侧向位移发展

根据式（16）和式（17）反算得到桩径不同时的模型参数α与tb，根据两种预测模型计算不同桩径桩基的位移累积率，将其计算结果与数模拟结果进行对比，如图18所示。由图可知，随荷载作用次数的增加，不同桩径桩基的泥面位移累积率均先增加后逐渐趋于稳定；且桩基的泥面位移累积率随桩径的增加而增大，当N=106时，较小直径（D=4 m、5 m）桩基的泥面位移累积率逐渐稳定于某一值，而直径较大（D=6 m、7 m、8 m）的桩基泥面位移累积率仍有进一步增加的趋势。对比桩基累积预测模型的计算结果发现，在桩径不同的条件下，幂函数型预测模型的计算结果与数值模型的结果在整个加载周期内的最大误差在2.92%以内，与数值计算值基本一致；在桩径较大的情况下（D=7 m、D=8 m），当N<105时，对数函数模型的预测结果与数值计算结果最大误差为4.01%，存在明显的偏差，在长期荷载作用后（N>105），对数函数模型的预测结果与数值计算结果的最大误差为7.73%。分析表明，在计算较大直径桩基的循环累积变形时，幂函数模型仍然具有较大优势。

图18 不同桩径单桩基础泥面累积变形

为进一步分析桩径对预测模型中参数α和tb的影响规律，图19分别绘制了两种形式预测模型中的参数α和tb与桩径D的关系。由图可知，幂函数模型和对数函数模型中的参数α和tb均与荷载幅值系数呈正比，但两条关系曲线的斜率不同，对数函数预测模型受桩基直径的影响偏大。对于桩径较大的桩基，利用对数函数模型预测桩基累积变形可能偏大。

图19 模型参数与桩径的关系

4 结 论

本文基于FLAC3D有限差分平台，根据循环试验结果编制刚度衰减程序，建立海上风机三维桩—土相互作用数值分析模型，研究单向水平循环荷载作用下荷载幅值、砂土相对密度和桩径大小对单桩基础累积变形特性的影响；根据整体数值分析模型的计算结果评价了幂函数和对数函数形式的桩基变形预测模型的适用性。相关结论如下。

（1）在长期水平循环荷载作用下，桩基的累积变形随循环荷载幅值的增加而增大；随场地砂土相对密实度的增加而降低；当无量纲荷载幅值系数相同时，桩基累积变形随桩径的增加而增大；总体上，桩基的累积变形随循环次数的增加先增加后逐渐趋于稳定。

（2）当循环荷载幅值、场地密实度增大时，桩基转动中心的埋深位置基本不变；当荷载无量纲系数相同时，桩基转动中心位置随桩径的增大呈现下移趋势。

（3）幂函数和对数函数模型的模型参数随荷载幅值和桩径的增大而增加，随场地土体密实度的增加而减小。

（4）在循环荷载作用前期（N<105），对数函数模型预测桩基变形较为准确，幂函数模型的计算结果偏小；随着循环荷载加载次数的增加，幂函数预测模型的预测结果更为准确；从桩基累积变形的长期发展来看，幂函数型预测模型更具有优势。