Large Diameter Monopile Pile-Soil Interaction in Sandy Seabed

近些年来，中国的海上风电产业快速发展。不占用陆上土地资源，风资源稳定，离能源消纳中心距离近等诸多优势使得海上风电成为了新能源开发的重要领域。相较于陆上风电，海洋环境复杂，风浪荷载作用大，运维难度高，安全可靠的海上风电基础是海上风电机组安全运行的关键保障^[1-2]。在所有的海上风电基础形式中，单桩基础应用最为广泛。在我国已建成的海上风电项目中，超过60%的项目全部或者部分选用了单桩基础^[3-4]。单桩基础加工制造容易，施工技术成熟，成本相对较低，在海上风电行业向深水区挺进的进程中大直径单桩基础仍将是工程人员的首选^[5-8]。

尽管大直径单桩基础在实际工程中得到了广泛应用，目前砂土中大直径单桩桩-土相互作用机理的研究尚不够透彻。行业设计人员主要使用API规范推荐的p-y曲线法（下简称API方法）对砂土中水平受荷桩进行分析，该方法是基于小直径柔性长桩的试验结果建立的^[9-10]，研究发现该方法在海上风电的大直径刚性短桩设计中高估了大直径单桩桩-土相互作用初始刚度^[11-15]。此外，API方法假设砂土中大直径单桩在加载过程是完全排水的^[9]。由于试验桩尺寸和加载速率的限制，大多数的试验研究如PISA项目的现场试验^{[6, 8]}和离心机试验^[13-14]也仅在完全排水条件下进行。随着研究的深入，越来越多的学者认为砂土中的水平受荷大直径单桩在单个荷载循环内的加载过程是不排水的^[16-18]。对砂土中大直径单桩的水平受荷过程，不同排水条件对桩-土相互作用产生何种影响，现有的研究尚未达成共识。

本研究选用有限元方法对大直径单桩水平受荷过程进行分析。有限元方法能模拟土的非线性特性，桩土接触界面和各类边界条件，较真实地对实际工程中可能产生的情况进行预测^[19-20]。然而，有限元法仅提供了桩-土相互作用分析的通用工具，选用合理的本构模型是确保有限元数值模拟成功的关键。

SANISAND（Simple ANIsotropic SAND）模型是Dafalias及其合作者开发的一系列基于边界面理论的弹塑性砂土模型的统称^[21-25]。SANISAND本构模型融合极限状态土力学的思想，其显著优点是可以使用同一组模型参数来模拟处于不同应力状态和相对密实度的砂土材料的力学性质，具有较好的工程实用潜能。

基于上述讨论，文章将采用SANISAND高级岩土本构，开展有限元分析，着重探索砂土中大直径单桩桩-土相互作用的两个关键问题：（1）砂土中大直径单桩桩土相互作用初始刚度以及与API方法对比；（2）排水条件对大直径单桩桩-土相互作用的影响。为此，论文首先介绍研究所选用的SANISAND本构模型的理论框架，利用室内土工试验数据标定代表性密实砂土的本构模型参数，将该模型应用于大直径单桩基础的有限元模拟，对上述两个问题开展深入研究。

采用Dafalias和Manzari于2004年发表的SANISAND模型^[22]，该版本采用了幂函数形式的极限状态线和指数形式的相变面和边界面方程。

本节将首先介绍p-q平面上的3个状态面：极限状态面、边界面和相变面，然后介绍三维空间中的弹塑性应力-应变关系，在本节最后列出了SANISAND模型的本构参数及具体的物理意义。

因SANISAND模型是有效应力模型，如不做特殊说明，文章中所有应力均指有效应力。

SANISAND模型基于极限状态土力学的框架，引入了边界面（Bounding Surface）和相变面（Dilatancy Surface）的概念，如图1所示。在持续的剪切作用下土最终达到极限状态，在极限状态下土的体积和剪应力不再随剪应变变化。在p-q平面，土体的极限状态面投影成一条直线，相应的相变面和边界面也是直线，三条线的斜率分别记为M，M_d，M_b。

Figure 1.  Yield, dilatancy, critical state and boundary surfaces on p-q plane^[25]

本模型选用Li和Wang^[26]提出的e-p空间幂函数形式的极限状态线，如式（1）所示：

式中：

p ——当前平均主应力（kPa）；

P_atm ——大气压强，取100 kPa；

e₀ ——p=0时的极限状态孔隙比；

λ_c、ξ ——极限状态面参数。

本模型中状态变量ψ=e－e_c为砂土孔隙比和对应于当前平均主应力的极限状态孔隙比的差值，表示当前状态与极限状态之间的距离，用以表征土的密实程度。当状态变量ψ=0，土体处于极限状态；若ψ>0，土体比极限状态更松散，剪切时将发生剪缩；若ψ<0，土体比极限状态更密实，剪切时将发生剪胀。

以排水三轴压缩试验为例，土调动应力比$\eta = {q \mathord{\left/ {\vphantom {q p}} \right. } p}$在剪切初期不断增加，直到达到峰值即边界面应力比M_b。边界面应力比是极限状态应力比M和状态变量ψ的函数，由式（2）表达：

式中：

n_b——边界面参数。

土的应力状态到达边界面之前（$\eta < {M_{\mathrm{b}}}$），土不断硬化至峰值抗剪强度；超过边界面之后（$\eta > {M_{\mathrm{b}}}$），土开始软化至残余抗剪强度。在这个过程中，土的应力状态最终趋于极限状态，最终M_b和η达到极限状态应力比M。

土的塑性硬化模量由调动应力比η和边界面应力比M_b之间的差值决定，如式（3）所示：

式中：

h——塑性硬化系数。

塑性硬化模量H的符号代表土的塑性响应状态，$H > 0$代表硬化阶段，$H < 0$代表软化阶段，$H = 0$表示土体达到了极限状态。

相变面的定义方式与边界面相似，相变面应力比M_d表示为极限状态应力比M和状态变量ψ的函数，如式（4）所示：

式中：

n_d——相变面参数。

剪胀变量表示剪胀发生的程度，模型中剪胀变量D的定义方式与塑性硬化模量H相似，与调动应力比η和相变面应力比M_d之间的差值有关，如式（5）所示：

式中：

A_d——剪胀状态参数，包含了组构效应的影响。

土的应力状态达到相变面之前（$\eta < {M_{\mathrm{d}}}$），剪胀变量D为正值，土处于剪缩状态；应力状态达到相变面之后（$\eta > {M_{\mathrm{d}}}$），剪胀变量D为负值，土进入剪胀状态。

如图2所示，在三维应力空间中，模型的极限状态面、边界面和相变面由在π平面上的映射点表示：这些应力点是由从静水压力轴出发的，与${\boldsymbol{n}}$方向相同的映射向量${{\boldsymbol{\alpha}}_{{i}}}$唯一定义的。由此边界面和相变面由映射向量来表示，如式（6）所示：

Figure 2.  Relationship between mapping points and mapping vectors^[27]

式中：

b和d ——分别代表边界面（Bounding Surface）和相变面（Dilatancy Surface）；

${\theta _{n\;}}$ ——向量${\boldsymbol{n}}$方向对应的应力洛德角，代表应力路径各向异性的影响，由对应的应力不变量J₂和J₃计算，如式（7）所示：

三维空间中当前应力状态与边界面和相变面之间的距离分别由向量${{{\boldsymbol{\beta}}}_{\text{b}}} = {{{\boldsymbol{\alpha}}}_{\text{b}}} - {{\boldsymbol{\alpha}}}$和向量${{{\boldsymbol{\beta}}}_{\text{d}}} = {{{\boldsymbol{\alpha}}}_{\text{d}}} - {{\boldsymbol{\alpha}}}$表示，该距离向量被用于计算塑性体应变和剪胀程度。

三维空间中SANISAND模型的弹性应力-应变关系由平均主应力$p = {{\left( {{\sigma _{11}} + {\sigma _{22}} + {\sigma _{33}}} \right)} \mathord{\left/ {\vphantom {{\left( {{\sigma _{11}} + {\sigma _{22}} + {\sigma _{33}}} \right)} 3}} \right. } 3}$和偏应力${\boldsymbol{s}} = {\boldsymbol{\sigma}} - p{\boldsymbol{I}}$表示，其中${\boldsymbol{{I}}}$是二阶单位张量，具体如式（8）的增量关系所示：

式中：

${\mathrm{d}}\varepsilon _q^e$ ——弹性偏应变增量；

${\mathrm{d}}\varepsilon _v^e$ ——弹性体应变增量；

G和K ——弹性剪切和体积模量，均为当前平均主应力p和孔隙比e的函数，如式（9）所示：

式中：

G₀ ——剪切模量参数；

e和ν ——土的孔隙比和泊松比。

如图3所示，三维空间中弹性区域以锥形屈服面限制，屈服面函数如式（10）所示：

Figure 3.  Yield surface of the SANISAND model^[27]

式中：

${{\boldsymbol{\alpha}}}$ ——偏背应力比张量；

m ——圆锥的张开半径，表示弹性区域的大小，一般m取值较小（小于0.05）。

当土的应力状态超过了屈服面所规定的弹性区域，土体将进入塑性应变阶段。土的塑性体应变增量$ {\mathrm{d}} \varepsilon _v^p$的计算与剪胀参数D相关，如式（11）所示：

式中：

A₀ ——无量纲的剪胀尺度缩放系数；

L ——塑性应变增量的大小，$\left\langle {} \right\rangle $为Macaulay括号，当括号中内容为正数时取它本身，为非正数时整个括号内容取0；

$ \left( {{{{\boldsymbol{\beta}}}_{{{\mathrm{d}}}}}{{: \boldsymbol n}}} \right) $——张量双点积，三维空间中的当前应力状态与相变面之间的距离；

${\boldsymbol{z}}$ ——组构张量，该模型中组构效应的影响，用于模拟循环加载，具体计算如式（12）所示。由式（12）可知式（11）中的$\left\langle {{\boldsymbol{z}}:{\boldsymbol{n}}} \right\rangle $在加载过程中保持为零，在卸载时发挥作用。

塑性总应变增量$ {\mathrm{d}}{\varepsilon ^p} $计算如式（13）所示：

式中：

L ——塑性应变增量的大小；

R ——塑性应变增量的方向，${\boldsymbol{R}} = \dfrac{1}{3}D{\boldsymbol{I}} + {\boldsymbol{R'}}$，其中含D项表示塑性体应变增量方向，$ {\boldsymbol{R}}' $表示塑性偏应变增量方向。

背应力比张量α的发展也由当前应力状态和边界面之间的距离决定，如式（14）所示：

式中：

h ——塑性硬化系数；

L ——塑性应变增量的大小。两个变量分别由式（15）和式（16）计算得来：

式中：

c_h和h₀ ——本构参数；

α_in ——初始背应力张量，以三轴应力状态下的调动应力比η与初始应力比η_in的差值表示。

式中：

K_p ——塑性模量，表示为$ {K_p} = \left( {{2 \mathord{\left/ {\vphantom {2 3}} \right. } 3}} \right)ph\left( {{{{\boldsymbol{\beta}}}_{\text{b}}}:{\boldsymbol{n}}} \right) $。

SANISAND模型的本构参数共计16个，可由常规的土工试验进行标定。本研究使用Blaker和Andersen于2015年发表的一组高质量密实砂土（Dr=80%）的土工试验数据^[28]对模型参数进行标定，最终选用的本构参数数值与具体物理意义在表1中总结。

本研究选取Blaker和Andersen于2015年发表的Dogger Bank sand土工试验数据^[28]标定SANISAND本构模型参数。Dogger Bank砂是一种典型的海洋石英砂，石英含量92%左右，颗粒形状呈次棱角状至次圆状，d₁₀和d₆₀分别为0.087 mm与0.174 mm。

土单元试验数值模拟中土样的竖向固结应力为200 kPa，侧向固结应力（三轴试验）为90 kPa（模拟K₀=0.45），与实际土工试验保持一致。文章仅考虑相对密实度为80%且细颗粒含量小于1%的试验结果。

本小节展示固结排水三轴压缩单元试验的模拟结果。从图4中可以发现，剪应力先增加并在轴向应变为0.8%时达到峰值（即峰值抗剪强度），随后剪应力逐渐减小并最终达到稳定值（即残余抗剪强度）。剪切模量参数G₀控制τ-ε_a曲线的初始阶段斜率，G₀增大时τ-ε_a曲线的初始阶段斜率加大；塑性模量系数h₀控制达到峰值抗剪强度时的应变，h₀增大时土在更小的应变水平下达到峰值强度；边界面参数n_b控制峰值抗剪强度和残余抗剪强度的比值，n_b增大时峰值抗剪强度和残余抗剪强度的比值增大。

Figure 4.  Shear stress-axial strain response from consolidated-drained triaxial compression test

图5中平台阶段的孔隙比取值即为土样达到极限状态时平均主应力p对应的极限状态孔隙比e_c，而e-ε_a曲线的形状和斜率主要由剪胀尺度缩放系数A₀控制。

Figure 5.  Void ratio-axial strain response from consolidated-drained triaxial compression test

本小节展示了固结不排水三轴压缩试验的模拟结果。在固结不排水三轴压缩试验中，调动应力比η=q/p是需要关注的变量，图6展示该变量在试验中的发展过程。从模拟结果中可以看出：随着平均主应力p的增加，调动应力比η先增大，继而在超越边界面时达到峰值；之后调动应力比减小，最终达到极限状态应力比M。此时，相变面、边界面和极限状态面重合（M_d=M_b=M=η），土的孔隙比达到极限状态孔隙比。

Figure 6.  Mobilized stress ratio-mean stress response from consolidated-undrained triaxial compression test

如图7所示，对固结不排水三轴压缩试验的应力-应变曲线中，塑性模量系数h₀也起到控制作用。如图7蓝色虚线所示，该值设置较小时计算得到较小的塑性模量H，调动应力比η下降的速度也较慢：当平均主应力p达到极限状态孔隙比对应的临界值时，调动应力比η大于极限状态应力比M和边界面应力比M_d，不排水应力-应变曲线出现峰值应力后的软化阶段，这一现象与实际情况是不相符的。

Figure 7.  Shear stress-axial strain response from consolidated-undrained triaxial compression test

图8为固结不排水三轴压缩试验中孔压随轴向应变的变化，泊松比ν的大小依据不排水试验的孔压预测结果进行修正：ν增大时体变刚度K增大，孔压上升。

Figure 8.  Pore pressure-axial strain response from consolidated-undrained triaxial compression test

以上展示了SANISAND本构模型参数标定的部分结果。在单元模拟中整理了一套依据实际土工试验结果标定各个SANISAND本构参数的方法，展示了SANISAND模型所选用的各个本构参数对数值模拟的控制作用。

除此以外，本研究还利用三轴拉伸与单剪试验的结果对SANISAND模型的模拟能力进行检验。对比发现由三轴压缩试验标定的模型参数对三轴压缩路径的试验结果吻合程度较高，对三轴拉伸和单剪路径的单元试验，需调整部分本构参数的取值以获得更佳的模拟效果。在实际工程问题中，各种应力路径可能同时存在。因此，在标定本构参数时应关注所研究问题中最重要的应力路径，调整参数取值以提升关键应力路径的模拟效果。针对大直径基础水平桩-土相互作用问题，水平加载的三轴拉伸路径最为重要，故在确定边值问题模型参数时将表1中（p=0时的极限状态孔隙比）e₀调整为0.7，其余参数不变。

数值模拟依托商用有限元软件ABAQUS开展，大直径单桩有限元网格如图9所示。模型由两部分组成，其中大直径单桩（直径D=10 m，长度3D=30 m）划分为6144个单元，单元类型C3D8；土体（直径为6D=60 m，深度4D=40 m）划分为22 080个单元，单元类型C3D8P。由于考虑的问题具有对称性，取模型的1/2（沿直径所在平面划分）作为研究对象。开展完全排水和完全不排水（仅土体的上表面设为自由排水边界）两种排水条件的数值模拟。桩顶高度在泥面处（Z=0 m），施加的荷载包括沿X轴正方向的水平荷载H=55 MN和绕Y轴正方向的弯矩M=1 650 MN·m，荷载施加在桩顶中心位置（图中点RP-1）。模型底面限制Z方向位移，侧面（弧面）限制X和Y方向位移，Y=0面设置为对Y方向的对称面，桩土接触面设置摩擦接触（μ=tanφ_c=0.7）并允许脱开。

Figure 9.  Finite element mesh of large diameter monopile foundation

模拟中将壁厚t=100 mm，直径D=10 m的钢管桩等效为同等直径的实心桩，以钢管桩的抗弯刚度（EI）等效计算实心桩的弹性模量，计算方法如式（17）所示：

式中：

E_eq ——等效弹性模量；

E ——钢材的弹性模量，取201 GPa；

d ——钢管桩内径（m）。

计算获得实心桩的等效弹性模量E_eq=16 GPa。桩身材料泊松比ν取0.3。

土的有效重度$\gamma ' = 10\;{\mathrm{kN}}/{{\mathrm{m}}^3}$，渗透系数$k = 5 \times {10^{ - 4}}\;{\mathrm{m}}/{\mathrm{s}}$。为便于地应力平衡，桩重度设置为$\gamma = 10\;{\mathrm{kN}}/{{\mathrm{m}}^3}$，与土有效重度相同。考虑该边值问题中最重要应力路径为三轴拉伸路径，选用的SANISAND本构模型参数中e₀取值为0.7，其余参数按表1取值。

本小节对大直径单桩在水平和弯矩荷载作用下桩-土相互作用开展分析。故本小节首先对大直径单桩在不同排水条件下水平受荷响应进行对比分析，然后将完全排水条件下的模拟结果与工程设计中常用的API方法计算得到的结果进行了对比。

图10为桩顶水平位移（参考点RP-1）随桩头荷载变化曲线。对比可以发现：桩顶水平荷载小于14.47 MN时，不排水加载与排水加载的力-位移曲线几乎重合；桩顶水平荷载大于14.47 MN而小于46.20 MN时，不排水加载模拟结果刚度更低；桩顶水平荷载超过46.20 MN后，不排水加载刚度明显大于排水加载；在桩头水平荷载H=55 MN时，不排水加载桩头水平位移为0.34 m，排水加载桩头水平位移为0.40 m。

Figure 10.  Variation of pile head displacement with applied load

图11为桩顶水平荷载小于14.47 MN时桩顶水平位移随桩头荷载的变化曲线，可以发现排水条件对桩-土相互作用几乎没有影响。

Figure 11.  Variation of pile head displacement with applied load under relatively low load levels

为解释图10和图11中展示的排水条件对桩-土相互作用刚度的影响，选取桩身右侧被动区潜在滑裂面上的特征单元进行研究。

图12～图13分别展示了该特征单元的剪应力（在此定义为q=σ_z－σ_x）以及孔压（完全不排水）随X方向正应变的发展过程，图14展示了两种排水条件下的有效应力路径。可以发现该特征单元在两种排水条件下的应力-应变曲线对比关系与泥面处荷载-位移曲线对比关系十分相似：在轴向应变较小时，排水条件对应力-应变曲线几乎没影响；随着轴向应变增大不排水加载应力在一段范围内始终低于排水加载；轴向应变继续增大，不排水加载应力很快提升并超过排水加载。图13和图14解释了不排水加载应力变化的原因：加载初期不排水加载单元发生剪缩，产生正孔压，降低土骨架的有效应力，导致不排水加载单元的刚度低于排水加载；随着轴向应变的继续增大，不排水加载单元进入剪胀阶段，正孔压消散随后转为负孔压，不排水加载单元的刚度很快提升且超过排水加载单元。

Figure 12.  Comparison of stress-strain responses at the characteristic element under different drainage conditions

Figure 13.  Development of pore pressure with strain in the X direction of the characteristic element under undrained conditions

Figure 14.  Comparison of stress paths of the characteristic element under different drainage conditions

图15为桩顶水平荷载为H=55 MN，弯矩M=1 650 MN·m时的桩身位移、桩身截面剪力和弯矩以及土抗力随深度的变化曲线。可见，在该荷载水平下，不排水加载条件下桩身位移小于排水加载条件，整体刚度更高。在较浅土层中，不排水加载条件桩身截面弯矩低于排水加载条件，且桩身剪力的衰减速度更快，这表明在浅层土中不排水加载下土抗力更大，这一点与图15(d)的土抗力曲线结果吻合。

Figure 15.  The variation of lateral deflection, cross-section bending moment, cross-section shear force and soil resistance varying with depth

综合所述，排水条件对大直径单桩桩-土相互作用的刚度具有显著影响，但影响的程度随荷载水平变化。在相对较小的荷载水平下，排水条件的影响几乎可忽略不计；随着荷载水平提高，完全排水条件下桩-土相互作用刚度超过不排水加载；然而，随着密实砂土剪胀特性的发挥，不排水条件下的桩-土相互作用刚度反超完全排水加载。从实际工程设计的角度而言，在文章所考虑的密实砂土地基中，设计荷载的水平位于图10中力和位移曲线的初始阶段，排水条件的影响可忽略不计，按照完全排水条件进行单桩桩-土相互作用设计是可行的。

API规范推荐的p-y曲线法（下简称API模型）假定大直径单桩水平受荷过程是完全排水的，故本节将对比完全排水条件下有限元数值模型预测结果与API模型预测结果。

图16为有限元模型与API模型预测的桩顶（泥面处）水平位移随桩头荷载变化的对比。可以发现API模型预测的桩-土相互作用刚度远高于有限元模型。当桩顶水平荷载为55 MN时，有限元模型预测泥面处水平位移为0.396 m，而API模型预测的泥面处水平位移仅为0.145 m。

Figure 16.  Comparison of pile head load-displacement responses predicted by finite element analysis and API p-y curves

从图16的对比中可以发现在荷载较低时，有限元模型与API模型的预测结果吻合度较高。图17为荷载低于5.8 MN时的预测结果对比：当水平荷载低于3.2 MN时，API模型和有限元模型的预测结果非常接近，然而水平荷载超过3.2 MN后，API模型预测的刚度迅速超过有限元模型结果。

Figure 17.  Comparison of pile head load-displacement responses predicted by finite element analysis and API p-y curves at relative low load levels

图18为各深度有限元模型和API模型预测的土反力-水平位移（p-y）曲线对比，图19为对应的小位移水平下的两种模型预测结果对比。图中散点为有限元模型预测结果，实线为API模型预测结果，相同颜色表示同一深度。对比图中结果可以发现：（1）随着深度的增加，两模型预测的p-y曲线刚度均上升；（2）在较小的位移水平下（小于0.005 m），两模型预测的曲线刚度很接近，但随着位移的增大，API模型预测的p-y弹簧刚度迅速超过有限元模型预测结果，证明API方法预测的p-y弹簧初始刚度过高，与现有研究结果相符^{[3, 11-15]}；（3）在深度不大于5 m的土层中API模型预测的极限土反力低于有限元模型结果。

Figure 18.  Comparison of p-y curves predicted by the finite element model and the API p-y curves

Figure 19.  Comparison of p-y curves predicted by the finite element model and the API model at low load levels

针对砂土中大直径单桩桩-土相互作用研究中的两个关键问题：（1）桩土相互作用的初始刚度以及与API方法的对比关系；（2）排水条件对桩-土相互作用的影响，本研究采用SANISAND高级岩土本构模型开展有限元数值模拟研究，得到了以下结论：

1）SANISAND模型作为一个成熟的高级本构模型，理论上具备使用同样一套参数对不同应力水平、孔隙比的土体和不同加载应力路径进行模拟的能力。然而与Blaker和Andersen公布的高质量土工试验数据^[28]对比表明：使用同一套参数模拟不同排水条件和不同应力路径下的试验存在一定难度。在实际应用SANISAND模型模拟具体的工程问题时，需针对所研究问题中最重要的应力路径有所侧重，合理选择模型参数。

2）排水条件对大直径单桩水平桩-土相互作用存在影响，且影响程度与荷载水平有关。本研究数值模拟表明，在较低荷载水平下，排水条件造成的影响几乎可以忽略不计：桩顶水平荷载小于14.47 MN时，排水条件不同的桩头荷载-位移曲线几乎重合；随着荷载提高，土开始剪胀，排水条件的影响变大：桩顶水平荷载大于14.47 MN而小于46.20 MN时，不排水加载刚度更低；桩顶水平荷载超过46.20 MN后，不排水加载刚度明显大于排水加载；在桩头水平荷载H=55 MN时，不排水加载桩头水平位移为0.34 m，排水加载桩头水平位移为0.40 m。对本研究所考虑的密实砂土地基，在实际工程设计关心的荷载水平下排水条件对桩-土相互作用刚度造成的影响可忽略不计，当前实际工作中按照完全排水条件进行大直径单桩的桩-土相互作用设计是可行的。

3）在较小的位移水平下（mm量级），有限元模型与API模型预测的p-y曲线刚度吻合度较好，但随着荷载水平增加，API p-y曲线高估桩-土相互作用刚度。