0.   引言

近年来，随着中国加大对海洋资源的开发和利用，桩基础被大量用于大型跨海桥梁、海洋油气平台以及海上风电场等工程。循环打入是桩基础的主要成桩方式，打入桩的承载特性是工程分析和设计的关键。揭示桩打入过程中桩土相互作用机理，深入认识桩打入引起的桩周土体应力场分布，为准确预测其承载特性提供依据，并对实际工程分析和设计具有重要意义。

过去几十年中，为了探究桩打入过程中的桩土作用机理，许多学者进行了大量的现场试验和室内标准模型槽试验，为建立更可靠的桩基础设计方法提供了关键依据。Lehane等^[1]在闭口桩桩身布置传感器，开展了系列的打入桩现场试验，发现桩身应力与距桩尖的相对位置有关（即h/R效应，其中h为测点与桩尖的竖向距离，向上为正，R为桩半径）。Chow^[2]利用相同的试验桩进行了现场试验，并利用预先打入土中的试验桩桩身上的传感器，测量了砂土中应力的变化规律，发现其同样具有显著的h/R效应。张明义等^[3]将预制桩静压入层状地基中，发现压桩力与锥尖阻力密切相关。Gavin等^[4]进行了开口桩的室内模型槽试验，并利用砂土中预埋的2个应力传感器，在室内试验中验证了应力的h/R效应。王浩等^[5]和周健等^[6]在室内模型槽中研究了静压桩的沉桩过程，并对桩端周围土体的位移场、应力场和孔隙率变化等进行了分析研究。Yang等^[7]在室内模型槽试验中对桩打入之后桩端与桩侧土的性状进行了研究，在桩侧发现了非常明显的颗粒破碎现象。Jardine等^[8-9]同样进行了标准模型槽试验，通过在桩身与砂土中安置大量的应力传感器，获得了闭口模型桩打入过程中桩身的应力变化规律与砂土中的三向（径向、环向和竖向）应力场，再次验证了h/R效应，同时提供了非常宝贵的桩侧砂土应力场数据。

然而实际的场地试验和室内标准试验涉及到土体的大变形以及桩土接触的界面滑移和颗粒破碎等众多复杂的物理现象，难以全面且准确地测量数据以及构建精准的理论，因此采用数值模拟研究砂土中桩打入过程十分必要，合理的数值模拟可以呈现土体变形的全过程，揭示现象发生的机理，为建立更合理的桩基础设计方法提供参考依据。许多学者已成功采用有限元法、离散元法和物质点法等数值模拟技术研究了桩的打入过程。Henke等^[10]使用CEL（Coupled Eulerian- Lagrangian，耦合欧拉-拉格朗日）有限元方法模拟了闭口桩的打入过程，得到了桩侧应力与砂土孔隙比的变化趋势。王腾等^[11]同样使用CEL方法模拟了钢管桩的大变形沉桩过程，分析了土塞的演变和形成机制。Zhang等^{[12, 13]}使用ALE（Arbitrary Lagrangian-Eulerian，任意拉格朗日-欧拉）有限元方法与颗粒破碎砂土本构模型模拟了Jardine等^[8-9]的室内模型槽试验，得到了砂土的颗粒破碎分布与桩侧的应力场。Yang等^[14]使用ALE方法与基于砂土平均正应力、塑性剪切应变和孔隙比的修正莫尔-库仑模型模拟了Jardine等^[8-9]的室内模型槽试验，得到了更为合理的桩侧应力场。周健等^[15]和詹永祥等^[16]基于颗粒流理论，采用离散元方法从微观角度研究了砂土中土塞的形成和发展。Wang等^[17]和Liu等^[18]使用二维离散元方法分别研究了桩打入过程中砂土的颗粒破碎分布与桩侧应力分布情况，模拟结果与试验结果较为一致。Lorenzo等^[19]使用物质点方法和改进的修正剑桥模型模拟了Jardine等^[8-9]的模型槽试验，获得了合理的应力模拟结果。

上述的数值模拟研究中，桩的打入过程被简化为单调压入，与实际的循环打入过程存在差异，因此具有一定的局限性。Ciantia等^[20]使用三维离散元模拟了模型桩的循环打入过程，采用将桩缓慢提升至桩顶力为零模拟稳态阶段，通过该方法成功获得了稳态阶段的砂土应力，模拟结果与试验结果较为一致。

然而离散元法计算量巨大，在实际工程应用中存在困难，适合小尺度问题的定性分析。有限元法等基于连续介质力学的方法依然是工程问题分析中的主流方法，其在分析计算时需要选取合适的土体本构模型。莫尔-库仑模型由于简单易用，被广泛用于砂土的有限元计算分析。如Yang等^[14]使用修正莫尔-库仑模型模拟了桩的单调打入过程，其中模型参数随打桩过程变化，需通过迭代确定不断更新，在模拟应力变化频繁的循环打入过程时较为不便。本文为了兼顾数值计算的稳定性与合理性，在莫尔-库仑本构模型的基础上开发了状态相关莫尔-库仑模型，形式简单且可自动反映应力状态及密实度对土体力学特性的影响。

本文考虑砂土中的桩循环打入过程，参照Jardine等^[8-9]的室内模型槽试验进行数值模拟研究。模拟采用大变形有限元ALE法，网格和材料相互分离且可通过频繁网格重划分在计算过程中始终保持高质量的网格，很好地解决了传统有限元计算中由于大变形而产生的网格畸变问题；同时使用状态相关莫尔-库仑模型，可较好地反映应力状态及密实度对砂土强度与变形特性的影响。模拟分析得到桩循环打入过程中全阶段的应力场分布规律，并与模型槽试验结果进行对比验证。

1.   ICL-INPG模型槽试验与有限元模型

ICL（英国帝国理工学院）-INPG（法国格勒诺布尔理工学院）模型桩试验槽布置如图 1所示，槽直径1.2 m，高1.5 m，模型桩尺寸与标准CPT（静力触探试验）仪保持一致，锥尖角度为60°，直径D为36 mm（R = 18 mm），采用不锈钢材质制造；模型槽中的砂土为中密实度（e₀ = 0.62，D_r = 72%）的枫丹白露干砂（Fontainebleau NE34 Sand）。为了测量桩打入过程中砂土内部的应力变化，落砂制样过程中及时布置大量微型土应力传感器，用于获得距桩轴线径向距离为2R~20R范围内砂土的径向应力${\sigma _r}^\prime$、环向应力${\sigma '_\theta }$和竖向应力${\sigma '_z}$。

图  1  模型槽试验布置图（改自Jardine等^[9]）

Figure  1.  Schematic diagram of calibration chamber tests (after Jardine et al^[9])

下载: 全尺寸图片 幻灯片

ICL-INPG模型槽试验中，砂土首先在150 kPa的超载下进行了一段时间的预压。为了提高数值计算效率，在有限元分析平台ABAQUS中采用二维轴对称的ALE方法，对桩的打入过程进行模拟。桩与砂土都采用对称边界条件进行模拟，其中桩模拟为拉格朗日刚体，尺寸与模型桩相同，锥尖60°，直径D = 2R = 36 mm，桩长L设为1.5 m，最终的打入深度L_p = 1m（L_p/R ≈ 56）；砂土为可变形体，尺寸与模型槽相同，高度H为1.5 m，宽度w = 0.6 m（w/R ≈ 33.3），土体重度γ = 16 kN/m³。砂土整体均为ALE的网格重划分区域，划分参数使用ABAQUS中ALE方法的默认值。桩-土间接触选用面-面接触，刚性桩外表面设为主面，可变形砂土设为从面，摩擦系数为0.49，与模型试验一致。

2.   状态相关莫尔-库仑模型与参数选取

由于砂土力学行为的状态相关性，在桩打入过程中，桩周砂土状态（应力、孔隙比等）经历了较大变化，位于桩尖、桩侧附近砂土的摩擦角和剪胀角等参数与初始状态相比有较大差异，采用相同的参数无法准确模拟桩打入过程中土体的力学响应。为了考虑这一特性，本文将Been等^[21]提出的砂土状态参数$\psi = e - {e_{\text{c}}}$引入莫尔-库仑模型中，采用改进后的状态相关莫尔-库仑模型模拟桩的循环打入过程。

砂土的临界状态孔隙比计算公式采用Li等^[22]提出的线性形式：

式中，e_Γ，λ_c，ξ为砂土的临界状态相关参数，可通过常规三轴试验获得，${p_{\text{a}}}$为大气压，取101.33 kPa。

建立状态相关莫尔-库仑模型的关键在于确认砂土的状态相关剪胀角以及状态相关内摩擦角，为了合理反映土体的状态相关剪胀性，Li等^[23]提出了以下剪胀方程：

式中：d₀和m为土体的模型参数；η和M分别为应力比和临界状态应力比。在莫尔-库仑模型中，土体只有在应力比达到临界状态应力比并进入屈服状态后才会发生塑性剪胀。因此可将η = M代入式（2），将其退化为

因此，考虑状态相关的剪胀角${\varphi _{\text{d}}}$为

同时密砂在排水状态下会出现显著的应变软化现象。在剪切过程中，密砂的应力-应变曲线会在达到峰值后下降并逐渐稳定，砂土的内摩擦角可能会超过临界状态摩擦角${\varphi _{\text{c}}}$，而常规的莫尔-库仑模型则无法模拟上述现象。因此采用类似的方法考虑状态相关内摩擦角${\varphi _{\text{f}}}$：

式中，${M_{\text{b}}}$为土的峰值应力比（Li等^[23]），而${n_{\text{b}}}$和${n_{\text{d}}}$同样是土的模型参数，它们分别反映了内摩擦角以及剪胀角与状态参数之间的关系。采用式（4），（6）对剪胀角和内摩擦角进行修正，可以采用同一套参数模拟砂土的硬化/软化，剪胀/剪缩现象。

根据Aghakouchak^[24]的室内试验标定结果，以及Dano等^[25]、Gaudin等^[26]、Andria-Ntoanina等^[27]和Altuhafi等^[28]的枫丹白露干砂三轴排水试验结果拟合，可以确认状态相关莫尔-库仑本构模型参数如表 1所示，其中E为杨氏模量，$\nu$为泊松比。

表  1  状态相关莫尔-库仑模型参数

Table  1.  Sand parameters adopted for state-dependent Mohr- Coulomb model

弹性参数	临界状态参数	状态相关参数
E = 60 MPa	φc = 30°	nd = 1.70
$\nu$ = 0.3	eΓ = 0.9	nb = 2.30
	λc = 0.019
	ξ = 0.23

下载: 导出CSV 

| 显示表格

3.   结果与分析

3.1   锥尖阻力q_c

在打入桩问题研究中，锥尖阻力${q_{\text{c}}}$与桩的承载力密切相关，是重要的研究指标，通过对比模拟得到的${q_{\text{c}}}$值与试验值，可验证模型的可靠性与所用本构模型及其参数选取的合理性。图 2为本文模拟所得到的锥尖阻力与试验数据的对比。

图  2  锥尖阻力q_c结果对比

Figure  2.  Comparison of q_c

下载: 全尺寸图片 幻灯片

可以发现数值模拟结果与试验数据吻合较好，在归一化打入深度L_p/R > 20时达到稳定值19 MPa，略低于模型槽试验值21 MPa。从图 2中还可以发现数值模拟得到的${q_{\text{c}}}$数据波动整体控制在±0.5 MPa范围内，表明了该方法在模拟桩打入过程大变形问题上有较好的数值稳定性。为了减少模型结果与试验结果的差异，本文的应力分析均采用锥尖阻力${q_{\text{c}}}$进行归一化处理。

3.2   桩周砂土状态

在桩打入过程中，随着砂土状态（应力、孔隙比）的改变，砂土的力学行为也会相应变化。图 3~5分别给出了桩完成打入后砂土的平均有效正应力、体积应变和状态参数等值线图，其中砂土在桩打入前由于自重与超载导致其初始状态参数约为−0.16。

图  3  桩打入完成后砂土平均有效正应力等值线图

Figure  3.  Contours of mean normal stress of sand after pile penetration

下载: 全尺寸图片 幻灯片

图  4  桩打入完成后砂土体积应变等值线图

Figure  4.  Contours of volumetric strain of sand after pile penetration

下载: 全尺寸图片 幻灯片

图  5  桩打入完成后砂土状态参数等值线图

Figure  5.  Contours of state parameter of sand after pile penetration

下载: 全尺寸图片 幻灯片

由图 3可知，桩在打入过程中，桩端附近砂土的平均有效正应力非常大（> 3 MPa），但其桩端砂土的体积应变很小，甚至在一定区域内接近零，呈现剪缩趋势。如图 4所示，由于桩端附近砂土处于高应力状态，导致该区域内的砂土达到了临界状态（$\psi = 0$）。当随着桩尖继续刺入，桩端附近砂土会发生向下流动，并停留在桩侧，导致桩侧砂土出现非常显著的应力松弛现象（≈ 0.2 MPa）。由于受到持续剪切作用，桩侧砂土发生了比较大的体积膨胀而仍然处于临界状态。同时由于体积应变在桩侧的变化梯度较大，从而导致2 < r/R < 3范围内的桩侧砂土在平均有效正应力几乎不变的情况下，状态参数也发生了显著变化。而在r/R > 8的区域，砂土的平均有效正应力与体积应变都非常低，砂土的状态参数接近于初始值，表明该区域的砂土力学行为几乎不受桩打入影响。

本文模拟得到的平均有效正应力与体积应变的分布规律与Yang等^[14]和Arshad^[29]等的结果基本一致。但本文的模拟中桩侧砂土始终处于剪胀状态，无法模拟Yang等^[7]在试验中发现的桩侧砂土颗粒破碎形成剪切带导致砂土变密的现象，采用考虑砂土颗粒破碎本构模型可以得到更符合试验的模拟结果。

3.3   径向应力${\sigma '_r}$

Jardine等^[8-9]的统计结果显示砂土径向应力的测量结果要明显优于另外两个方向，因此本文重点分析砂土中的径向应力分布。与以往研究不同，本文模拟考虑了桩的循环打入过程，在每个打入段后将桩顶荷载撤为0，使桩达到稳态后才继续下一步打入操作，因此本文分别考虑打入过程中的径向应力${\sigma '_{{\text{rm}}}}$及稳态阶段的径向应力${\sigma '_{{\text{rs}}}}$。图 6给出了桩打入过程中距离桩中心轴径向距离为r/R = 2位置处，3个不同的深度（顶层：z/R = 10.6，中层：z/R = 30.6和底层：z/R = 46.1）的归一化${\sigma '_{{\text{rm}}}}$结果对比。3个深度处的结果均显示出与打入深度L_p明显的高度依赖性，可发现${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$随L_p变化而剧烈变化。取某一固定深度分析，在桩刚打入时，桩尖离该深度的距离还较远，土体的应力水平较低；随着桩尖逐渐靠近，${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$急剧增大，并在桩尖抵达该深度时达到峰值；随着桩的继续打入，土体的应力水平迅速下降，最终稳定在某一略高于土体初始应力的值。从峰值结果看，有限元的模拟结果略低于试验结果，但总体上两者曲线的变化趋势基本一致。

图  6  桩打入过程中距桩轴r/R = 2处归一化径向应力${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$结果对比

Figure  6.  Comparison of ${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$ at r/R = 2 at driving stage

下载: 全尺寸图片 幻灯片

图 7展示了3个位置处砂土归一化径向稳态应力${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$的结果。可以发现稳态应力的变化趋势与打入应力的变化趋势基本相同，${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$随桩尖的靠近而逐渐增大到最大值，在桩尖经过之后衰减到某一较小的平台值。曲线的变化趋势也基本相符，但有限元的模拟结果略高于试验结果。

图  7  桩打入过程中距桩轴r/R = 2处归一化径向稳态应力${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$结果对比

Figure  7.  Comparison of ${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$ at r/R=2 at pause stage

下载: 全尺寸图片 幻灯片

图 8绘制了砂土的径向应力随着h/R变化的趋势曲线。其中图 8（a），（c）分别给出离桩轴径向距离r/R = 2处${\sigma '_{{\text{rm}}}}/{q_{\text{c}}}$与稳态应力${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$的模拟结果随h/R变化的曲线，图 8（b），（d）则分别展示了对应的试验结果。图 8（e），（f）分别给出了离桩轴径向距离r/R = 3位置处，砂土归一化稳态应力${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$的模拟与试验结果。从上述组图中可以发现，应力的变化趋势均较为相似：当桩打入深度在测点之上时（h/R为负值），砂土的打入应力与稳态应力均维持在比较低的初始水平（0.3%q_c~0.4%q_c）；随着桩的继续打入，砂土的应力水平不断提高，并在h/R = 0（即桩尖抵达测点深度）附近达到峰值；当桩尖超过测点深度时（h/R变为正数），此时应力又会迅速减小，直到稳定到某一平台值（1.0%q_c~2.0%q_c）。结果对比可以发现径向打入应力与稳态应力存在显著差别，且稳态应力模拟结果的数值稳定性比打入应力差，与Jardine等^[9]统计的试验中稳态应力的变异系数显著大于打入应力的趋势一致。图 8结果同时证明了径向稳态应力${\sigma '_{{\text{rs}}}}/{q_{\text{c}}}$在空间上也是h/R以及r/R的函数。

图  8  桩打入过程中距桩轴r/R = 2和3处归一化径向应力随h/R变化结果对比

Figure  8.  Comparison of normalized radial stresses at h/R= 2 or 3

下载: 全尺寸图片 幻灯片

图 9为桩打入完成之后，桩尖上方区域（5 < h/R < 41）砂土的径向稳态应力随r/R变化的模拟与试验结果。结果显示，径向稳态应力的模拟结果虽然偏大，但其整体变化趋势与试验结果基本一致，应力的最大值并不是出现在桩侧（r/R = 1），而是在土体内部。以图 9（a）模拟结果为例，径向稳态应力在桩侧（r/R = 1）位置处还相对较小，仅为1.8%q_c~2.8%q_c；而在r/R =3~4范围内，径向稳态应力达到峰值（2.8%q_c~4.4%q_c），约为桩侧值的1.5~2.0倍；之后随着r/R的继续增大，径向稳态应力逐渐减小，直到r/R > 30后，应力水平稳定在0.8%q_c左右。同时从模拟与试验结果对比中可以发现，距桩轴相同水平距离（r/R）的砂土，其径向稳态应力会随h/R的增大而减小，验证了径向稳态应力的h/R效应。

图  9  桩打入完成后桩尖上方区域归一化径向稳态应力随r/R变化结果对比

Figure  9.  Comparison of normalized radial stresses at pause stage with r/R after pile penetration

下载: 全尺寸图片 幻灯片

4.   结论

本文使用大变形有限元ALE方法，结合状态相关莫尔-库仑模型，模拟了砂土中桩的循环打入过程，与ICL-INPG室内模型桩试验结果对比，获得了一致的桩打入全过程应力演化与分布结果。

（1）所采用的大变形有限元ALE方法，可以较好地模拟桩在打入过程中的土体大变形，并获得相应的桩身和桩侧砂土中应力随桩打入过程的变化规律。

（2）桩打入过程中打入应力与稳态应力的数值模拟结果变化趋势均与模型槽试验结果较为一致，证明了ALE方法结合状态相关莫尔-库仑模型模拟桩打入过程的可靠性和准确性。

（3）数值模拟结果验证了当桩处于不同的循环阶段（打入阶段与稳态阶段）时，砂土的应力水平会有较大的差异，且稳态应力与打入应力一样，在空间上均与桩尖的相对位置（h/R，r/R）密切相关。

（4）由于高应力、大位移等原因，试验中距桩较近的区域（r/R < 2）应力无法准确测量，本文数值模拟可以较好地补充靠近桩身区域的应力结果，建立完整的打入桩桩周应力场分布，为打入桩的承载力计算分析和设计提供依据。