在有限元软件划分桩体网格时采用结构化的网格剖分技术，可以保证桩体单元成层堆叠，同一层单元的纵坐标数值是相等的。为此，首先根据桩单元的纵坐标大小，将桩身单元纵坐标相同的单元创建为一个集合，然后以两层单元交界面处的单元为参考，沿桩身不同高度建立一系列目标截面。获得目标截面后，通过添加关键字即可让软件在计算过程中自动产生相应截面的弯矩、位移、Mises应力等一系列数据。

总结发现，p-y曲线的提取思路主要有以下四种：

1)首先从计算结果文件中提取桩身各深度的挠度，插值形成挠度曲线，根据其曲率及桩截面刚度获得桩身弯矩，进一步推算出桩周土反力。桩身位移、弯矩及桩周土反力三者间存在如下关系：

式中：M为桩身弯矩(kN·m)；EI为抗弯刚度(kN·m²)；p为桩周土反力(kN·m^-1)；y为桩身挠度(m)。

2)同样根据以上关系，直接从计算结果中提取出桩身各深度的弯矩值，通过最小二乘法插值形成弯矩随深度变化的曲线，进一步推导出各深度的p值。

3)根据有限元结果，对不同深度桩身截面的应力张量通过面积分计算出各截面的剪力V，对剪力采用三次样条曲线沿桩长进行插值得到相应表达式，p值由V的表达式沿桩长度方向求导得到，表达式为：

4)根据有限元计算结果，沿不同深度处桩身各截面对应力张量沿圆周进行线积分，直接得到p值，其表达式为：

式中：S为应力张量；v为所计算p值的方向；n为桩身表面法线方向。

针对上述四种方法，笔者对各种方法的计算精度进行分析比较如下：

第1种方法相比于第2种方法，要先从桩基挠度间接获得桩身弯矩，显然没有直接提取桩身弯矩的精度高，因此不予考虑。第2种方法需要对桩身弯矩求两次导数，至少需要三次以上方程对桩身弯矩-深度关系进行拟合，求导和拟合都伴随着精度的降低。桩体根据刚度、埋深等的不同存在刚性桩和柔性桩的区别，不同实际工况又会导致边界条件存在区别。以上各种不确定因素导致拟合弯矩曲线最终获得p-y曲线的过程对数据分析人员的经验要求较高，很难得到准确p值。

上述第3、4种方法获得p值，降低了拟合阶数和求导次数，理论上精度明显优于第2种方法。有限元数值计算可以提取出截面剪力，对第3、4种方法有良好适用性。另外用三次样条曲线插值法代替最小二乘法，可以更好地利用已有数据点，减小边界条件不确定性所带来的计算结果不准确。下面将依据数值试验结果，对第3、4种方法进行比较，并相应给出选用建议。

对简单且易于得到理论值的悬臂梁进行研究，采用弹性材料模拟悬臂梁，尺寸为5 m×5 m×20 m，悬臂梁右端固定，在上表面施加1 MN/m²的均布荷载，分别采用第3、4种思路编写相应p-y曲线提取程序。结果发现采用第4种方法编写的程序计算，通过多次加密网格后距桩端不同距离的截面p值相对误差始终保持在7.5%以上，采用第三种思路编程计算的结果如表1所示：

由表1可看出p值计算结果距固定端近处相对误差较大，距固定端远处相对误差则可控制在3%以内。因此，虽然第4种方法能够更直接地得到p值，然而实际操作中其准确程度可能远不及第三种方法，原因可能在于第4种方法圆周积分涉及的点全部位于单元表面，应力由外推得到，而第3种方法面积分涉及的点位于梁内部，应力由内插得到。因此本文选择第3种方法来编写提取模块，该方法实现的流程如图1所示。

图  1  p-y曲线提取技术路线

Figure 1.  Technology roadmap for p-y curve extraction

采用有限元方法进行桩基水平加载的数值分析。桩基采用线弹性模型，地基土选用非线性的本构模型。桩和土采用独立的实体单元建模并采用6面体8节点线性减缩积分单元(C3D8R)划分网格。采用接触属性采用主-从(Master-Slave)算法，来定义桩土间接触属性模拟桩-土之间的相对位移和剪力传递，其中刚度大的桩基定义为主控面，桩周土表面定义为从属面。桩与土之间的接触面属性采用莫尔-库伦摩擦罚函数，界面相对滑动摩擦系数u＝tan(0.75φ)，其中φ为土体内摩擦角^[9]。整个地基模型为圆柱体形，直径20D(D为桩基直径)，高L＋10D(L为桩打入土深度，10D为下卧层厚度)，以上模型尺寸可忽略边界的影响^[10]。桩土相互作用问题中，距桩身5D范围以内土体对计算结果有较大影响，因此建模时对该部分网格进行加密，提高计算精度。

考虑到现场试验所提供的土体参数有限，在数值计算中采用相对简化的本构模型来模拟土体性质。表层粘土采用Von Mises屈服准则进行模拟^[11]，其性质完全由不排水抗剪强度控制；深层砂土采用非相关联的Drucker-Prager本构进行模拟。

不论是砂土还是粘土，它们的杨氏模量都和土体围压有关，如式(5)所示^[6]：

式中：E₀为土体在大气压强下的杨氏模量(MPa)；p为平均有效主应力(kPa)；P_a为大气压强(kPa)；a是和土体孔隙比有关的一个常量。

E₀可以通过三轴试验获得，根据杨氏模量的定义可知，E₀中应只包含弹性应变，因此选取土体应力应变曲线中的初始切线模量E_ini为对应围压下的E₀。

土的泊松比μ影响土体的变形特征，对计算结果影响较大^[12]。饱和软粘土的泊松比接近0.49。参考Vermeer等^[13]提出的三轴固结排水压缩试验双折线模型，由体应变与轴向应变关系曲线计算得到砂土试样泊松比为0.3。

为了检验本文所提出的p-y曲线提取方法在桩土相互作用问题中的适用性，以下将对有限元模型的计算结果和广东桂山海上风电场项目的现场试验结果进行对比。

现场GK04试桩，是桂山海上风电场工程的一根试桩。如图2所示，试桩为开口钢管桩，材料为Q345B，外径2.2 m，壁厚0.03 m，桩身泥面以上3.6 m到泥面以下6 m部分加厚，壁厚0.055 m^[2]。GK04试桩桩长70 m，所在地质土层如图2所示。现场软粘土的其他基本参数如表2所示。水平推力加载装置采用油压千斤顶(最大出力为2.5 MN)，打桩结束58 d后进行加载。分15个荷载等级，最大荷载为2.025 MN。

图  2  GK04试桩剖面图

Figure 2.  Cross-section of GK04

钢桩弹性模量E为210 GPa，外径2.2 m，内径2.14 m。根据等效刚度替换原则计算得到有限元模型桩的弹性模量21.99 GPa，加粗部分模型桩弹性模量38.95 GPa。根据质量等效计算得到正常部分实心桩的密度为419.65 kg/m³，加粗部分模型桩等效密度为760.5 kg/m³。有限元模型中的桩参数如表3所示。

通过桂山桩基现场实测数据和有限元分析结果对比，对计算结果的可靠性进行验证。为了方便对比，对有限元模型施加同现场试桩一样的水平荷载，分别提取出加载点的荷载-位移曲线、桩身弯矩曲线、桩身位移曲线等结果，并通过自行开发的p-y曲线提取程序计算出各深度的p-y曲线，结果如下。

API规范计算结果、有限元模型计算结果以及现场试验的荷载-位移曲线如图3所示。比较可知，有限元模型的计算结果与现场试验的桩顶水平荷载-位移曲线实测结果吻合较好，相同荷载下二者的位移值都明显小于API规范计算值，说明API规范偏保守；与现场试验相比，有限元模型在荷载相同的情况下，计算的位移值始终略大于实测的位移值，这是因为现场测试中，桩基周围土体存在重固结的现象^[14]。土体重固结现象随着加载的进行越来越明显，因此当荷载较小时，计算位移值和实测位移值差别较小；荷载随时间的推移逐渐增大，土体的重固结加重，二者差别也变大。

图  3  加载点的荷载-位移曲线

Figure 3.  Load-displacement curve of the mount point

不同荷载下数值模拟和现场试验得到的桩身弯矩对比如图4所示。从图中可以看出，不同等级加载条件下，数值模拟和试验结果均吻合较好，二者的结果都表明：随着荷载的增加，桩身最大弯矩所在位置逐渐向更深处移动。

图  4  桩身弯矩对比

Figure 4.  Comparison of bending moment for pile

计算不同荷载下桩身位移随深度的变化较为简单，由于在提取弯矩时已经沿桩身取了一系列的截面，只需在前处理模块作相应的设置，软件即可自动计算出这一系列截面中点的位移值，如图5数值模拟结果与实测结果的对比所示，不同荷载作用下二者的桩身变形吻合较好，考虑到土体在加载过程中的重固结效应，数值模拟计算的桩身位移总是略大于实测桩身位移。

图  5  试验与数值模拟结果对比

Figure 5.  Comparison between experimental and numerical simulation results

如图6所示为根据现场实测数据分析所得p-y曲线、数值计算p-y曲线、以及API规范法计算p-y曲线的对比。可以看出数值计算结果和现场实测结果比较接近，二者的初始刚度和极限土反力都明显大于API规范法的计算结果。

图  6  现场试验与离心模型试验结果对比

Figure 6.  Comparison between field test and centrifuge model test

数值计算和现场实测结果在加载范围内得到的p-y曲线吻合较好。随着荷载增大，数值方法所得的承载系数N_p(N_p＝P_u/s_uD，P_u为极限承载力，s_u为不排水抗剪强度，D为桩基直径)已经趋向稳定，N_p值最终稳定在9.5附近。Randolph^[15]基于理论分析得到N_p的塑性解最大值为11.94，Murff^[16]通过开展离心模型试验得到N_p＝7.0～15.8，Jeanjean^[17]通过开展离心模型试验得到N_p＝9.3～16.4，本文的数值计算结果位于上述各区间内，较为合理。而现场实测结果仍有增大的趋势，可以预见最终的N_p会大于数值计算结果。笔者认为出现这种现象的一个重要的原因是打桩所导致的土体重固结。

基于以上分析，采用本文所提出的数值计算模型以及提取p-y曲线的后处理程序所得到的分析结果是合理可信的。

1)本文总结了4种p-y曲线提取思路，针对弹性梁对比了各种思路的计算结果和理论值的差异，明确了通过桩身剪力获得p-y曲线的方法的优越性，并基于此二次开发了自动化的提取程序。

2)将二次开发程序运用于桂山海上风电场项目的钢管桩水平承载力计算，进行了荷载-位移曲线、桩身弯矩、桩身变形以及p-y曲线等内容的计算，与现场实验结果进行对比，二者吻合较好，验证了本文提出的方法在钢管桩水平承载力计算时的适用性。