-
海上风能资源是极具发展潜力的可再生能源之一，我国海上风能资源十分丰富，且具有风速稳定、受土地、环境噪声制约小等优点，根据“十三五”可再生能源规划，海上风电资源将继续获得大力开发。海上风机桩基础形式有单桩基础、三桩基础、四桩导管架基础和高桩承台基础等，其中单桩基础制造简单、结构受力明确、施工速度快，因此应用最为广泛。单桩作为海上风机基础时，需要承受上部结构风压、波浪和海流等水平荷载的耦合作用，其在水平荷载作用下的桩土相互作用机理是研究的关键[1]。
目前对于水平荷载作用下，桩-土相互作用分析的方法主要有m值法、 p-y曲线法、有限单元法等，其中m值法是一种线弹性地基反力法，主要适用于桩侧土体处于弹性变形阶段的求解，其计算结果往往和实际情况有一定差距[2,3]。 p-y曲线法是一种弹塑性地基反力法，能够真实地反应桩周土的非线性特性及在此基础上桩的内力和变形特性，可用于桩基大变形非线性分析，有效克服了m值法的缺陷，被公认为是分析桩土水平非线性变形最有效的方法。Reese等[4]通过砂土中直径0.61 m桩基的水平加载现场试验提出砂土中的分段 p-y曲线，这一方法经过改进后被美国石油协会(API)[5]采用。随着计算机技术的广泛应用，有限元法被越来越多地应用于桩-土相互作用分析中：Yang[6]分别以弹塑性材料和纯弹性材料模拟土和桩，建立三维数值分析模型，研究了土体分层特性对桩基水平承载力的影响；Abdel-Rahman[7]采用有限元法，针对近海风电超大直径单桩基础，研究了API规范中砂土 p-y曲线对超大直径单桩基础的适用性。
基于有限元软件进行 p-y曲线分析可以结合 p-y曲线法和有限元法的优点，具有巨大应用前景。然而现有有限元软件中都没有直接提取 p-y曲线的模块，目前 p-y曲线是基于桩身弯矩和位移同深度的关系通过多项式拟合和二次求导来提取[8]，根据桩体性质和土体性质的不同，需要求解者依据经验确定相应的边界条件和多项式次数，不同拟合求导后的曲线往往存在较大差异。鉴于此，本文基于编程开发了一套新的 p-y曲线自动提取方法，可直接获得桩身剪力的特点，采用三次样条曲线拟合桩身截面剪力后一次求导得到桩身反力p，减小了边界条件对拟合结果的影响，降低了使用者的经验要求，可以准确方便地自动提取 p-y曲线。结合桂山海上风电场的现场试桩数据进行对比分析，证明了本文提出的 p-y提取方法合理可行。
HTML
-
在有限元软件划分桩体网格时采用结构化的网格剖分技术，可以保证桩体单元成层堆叠，同一层单元的纵坐标数值是相等的。为此，首先根据桩单元的纵坐标大小，将桩身单元纵坐标相同的单元创建为一个集合，然后以两层单元交界面处的单元为参考，沿桩身不同高度建立一系列目标截面。获得目标截面后，通过添加关键字即可让软件在计算过程中自动产生相应截面的弯矩、位移、Mises应力等一系列数据。
总结发现，p-y曲线的提取思路主要有以下四种：
1)首先从计算结果文件中提取桩身各深度的挠度，插值形成挠度曲线，根据其曲率及桩截面刚度获得桩身弯矩，进一步推算出桩周土反力。桩身位移、弯矩及桩周土反力三者间存在如下关系：
((1)) ((2)) 式中：M为桩身弯矩(kN·m)；EI为抗弯刚度(kN·m2)；p为桩周土反力(kN·m-1)；y为桩身挠度(m)。
2)同样根据以上关系，直接从计算结果中提取出桩身各深度的弯矩值，通过最小二乘法插值形成弯矩随深度变化的曲线，进一步推导出各深度的p值。
3)根据有限元结果，对不同深度桩身截面的应力张量通过面积分计算出各截面的剪力V，对剪力采用三次样条曲线沿桩长进行插值得到相应表达式，p值由V的表达式沿桩长度方向求导得到，表达式为：
((3)) 4)根据有限元计算结果，沿不同深度处桩身各截面对应力张量沿圆周进行线积分，直接得到p值，其表达式为：
((4)) 式中：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/m2的均布荷载，分别采用第3、4种思路编写相应p-y曲线提取程序。结果发现采用第4种方法编写的程序计算，通过多次加密网格后距桩端不同距离的截面p值相对误差始终保持在7.5%以上，采用第三种思路编程计算的结果如表1所示：
距固定端距离/m 3.0 7.0 11.0 15.0 19.0 p相对误差/% 3.5 2.98 2.74 2.74 2.46 Table 1.
Comparison of p value 由表1可看出p值计算结果距固定端近处相对误差较大，距固定端远处相对误差则可控制在3%以内。因此，虽然第4种方法能够更直接地得到p值，然而实际操作中其准确程度可能远不及第三种方法，原因可能在于第4种方法圆周积分涉及的点全部位于单元表面，应力由外推得到，而第3种方法面积分涉及的点位于梁内部，应力由内插得到。因此本文选择第3种方法来编写提取模块，该方法实现的流程如图1所示。
Figure 1. Technology roadmap for p-y curve extraction
1.1 p-y曲线提取思路
1.2 p-y曲线提取思路比较与选取
-
采用有限元方法进行桩基水平加载的数值分析。桩基采用线弹性模型，地基土选用非线性的本构模型。桩和土采用独立的实体单元建模并采用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]：
((5)) 式中：E0为土体在大气压强下的杨氏模量(MPa)；p为平均有效主应力(kPa)；Pa为大气压强(kPa)；a是和土体孔隙比有关的一个常量。
E0可以通过三轴试验获得，根据杨氏模量的定义可知，E0中应只包含弹性应变，因此选取土体应力应变曲线中的初始切线模量Eini为对应围压下的E0。
土的泊松比μ影响土体的变形特征，对计算结果影响较大[12]。饱和软粘土的泊松比接近0.49。参考Vermeer等[13]提出的三轴固结排水压缩试验双折线模型，由体应变与轴向应变关系曲线计算得到砂土试样泊松比为0.3。
2.1 有限元模型
2.2 本构及关键参数取值规则
-
为了检验本文所提出的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。
Figure 2. Cross-section of GK04
有效重度γ′/(kN·m-3) 比重Gs 液限LL/% 塑限PL/% 渗透系数k/(cm·s-1) 压缩系数av/ MPa-1 压缩模量Es/MPa 6.7 2.68 43.8 25.5 2E-7 1.2 2.1 Table 2.
Parameters of soil 钢桩弹性模量E为210 GPa，外径2.2 m，内径2.14 m。根据等效刚度替换原则计算得到有限元模型桩的弹性模量21.99 GPa，加粗部分模型桩弹性模量38.95 GPa。根据质量等效计算得到正常部分实心桩的密度为419.65 kg/m3，加粗部分模型桩等效密度为760.5 kg/m3。有限元模型中的桩参数如表3所示。
桩径/m 埋深/m 桩长/m 密度/(kg·m-3) 泊松比 弹性量/GPa 正常段 加粗段 2.2 57.4 70 419.65 0.3 21.99 38.95 Table 3.
Parameters of pile -
通过桂山桩基现场实测数据和有限元分析结果对比，对计算结果的可靠性进行验证。为了方便对比，对有限元模型施加同现场试桩一样的水平荷载，分别提取出加载点的荷载-位移曲线、桩身弯矩曲线、桩身位移曲线等结果，并通过自行开发的p-y曲线提取程序计算出各深度的p-y曲线，结果如下。
-
API规范计算结果、有限元模型计算结果以及现场试验的荷载-位移曲线如图3所示。比较可知，有限元模型的计算结果与现场试验的桩顶水平荷载-位移曲线实测结果吻合较好，相同荷载下二者的位移值都明显小于API规范计算值，说明API规范偏保守；与现场试验相比，有限元模型在荷载相同的情况下，计算的位移值始终略大于实测的位移值，这是因为现场测试中，桩基周围土体存在重固结的现象[14]。土体重固结现象随着加载的进行越来越明显，因此当荷载较小时，计算位移值和实测位移值差别较小；荷载随时间的推移逐渐增大，土体的重固结加重，二者差别也变大。
Figure 3. Load-displacement curve of the mount point
-
不同荷载下数值模拟和现场试验得到的桩身弯矩对比如图4所示。从图中可以看出，不同等级加载条件下，数值模拟和试验结果均吻合较好，二者的结果都表明：随着荷载的增加，桩身最大弯矩所在位置逐渐向更深处移动。
Figure 4. Comparison of bending moment for pile
-
计算不同荷载下桩身位移随深度的变化较为简单，由于在提取弯矩时已经沿桩身取了一系列的截面，只需在前处理模块作相应的设置，软件即可自动计算出这一系列截面中点的位移值，如图5数值模拟结果与实测结果的对比所示，不同荷载作用下二者的桩身变形吻合较好，考虑到土体在加载过程中的重固结效应，数值模拟计算的桩身位移总是略大于实测桩身位移。
Figure 5. Comparison between experimental and numerical simulation results
-
如图6所示为根据现场实测数据分析所得p-y曲线、数值计算p-y曲线、以及API规范法计算p-y曲线的对比。可以看出数值计算结果和现场实测结果比较接近，二者的初始刚度和极限土反力都明显大于API规范法的计算结果。
Figure 6. Comparison between field test and centrifuge model test
数值计算和现场实测结果在加载范围内得到的p-y曲线吻合较好。随着荷载增大，数值方法所得的承载系数Np(Np＝Pu/suD，Pu为极限承载力，su为不排水抗剪强度，D为桩基直径)已经趋向稳定，Np值最终稳定在9.5附近。Randolph[15]基于理论分析得到Np的塑性解最大值为11.94，Murff[16]通过开展离心模型试验得到Np＝7.0～15.8，Jeanjean[17]通过开展离心模型试验得到Np＝9.3～16.4，本文的数值计算结果位于上述各区间内，较为合理。而现场实测结果仍有增大的趋势，可以预见最终的Np会大于数值计算结果。笔者认为出现这种现象的一个重要的原因是打桩所导致的土体重固结。
基于以上分析，采用本文所提出的数值计算模型以及提取p-y曲线的后处理程序所得到的分析结果是合理可信的。
3.1 现场试桩工程参数
3.2 计算结果对比验证
3.2.1 荷载-位移曲线
3.2.2 桩身弯矩
3.2.3 桩身变形
3.2.4 p-y曲线
-
1)本文总结了4种p-y曲线提取思路，针对弹性梁对比了各种思路的计算结果和理论值的差异，明确了通过桩身剪力获得p-y曲线的方法的优越性，并基于此二次开发了自动化的提取程序。
2)将二次开发程序运用于桂山海上风电场项目的钢管桩水平承载力计算，进行了荷载-位移曲线、桩身弯矩、桩身变形以及p-y曲线等内容的计算，与现场实验结果进行对比，二者吻合较好，验证了本文提出的方法在钢管桩水平承载力计算时的适用性。
DownLoad:
