摘要： 基于桩侧土体和桩端土体的变形与破坏机理不同，以及土力学 理论 、有关室内和现场试验结果，并为简化起见，分别采用双曲线模型、理想弹塑性模型来描述桩侧土体静阻力和桩端土体静阻力特性。考虑桩身自重及桩周土体阻力作用建立了一维动力打桩波动方程，采用有限差分法进行求解，编制了基于上述土体模型的打桩 分析 软件ADP，并用该软件对某海洋桩基平台的打桩工程进行了数值分析。结果表明，本文建立的模型更符合实际情况，有助于提高分析精度。
关键词： 双曲线模型 理想弹塑性 波动方程
1 概述
在运用波动方程法预测桩的可打入性及单桩极限承载力中，桩周土体静力模型的合理选择是个极其重要的 问题 。土体的静力特性远非线弹性、理想弹塑性能简单描述，而非线性、非弹性、弹塑性等模型可较好地描述。因此，改进土体静力模型及其 计算 参数的确定 方法 ，是进一步完善波动方程分析法的一个非常重要方面。
桩侧摩阻力的发挥一般是桩体和土体之间的剪切破坏，也可能是桩体带着部分土体，土体间的剪切破坏，而桩端阻力的发挥有的是 “ 刺入 ” 破坏，有的是 “ 压剪 ” 破坏 [1] 。由此可知，桩侧土主要承受剪切变形，而桩端土体变形主要是压缩，而且不能承受拉应力，桩侧土体和桩端土体的变形和破坏机理是截然不同的。
文献 [2] 通过室内剪切试验，测得不同法向压力下，钢和混凝土材料分别与土之间的摩阻力与剪切位移的关系曲线，用以描述桩、土间的荷载传递特性。结果表明摩阻力和剪切位移呈非线性关系，而且符合双曲线方程。汉森 (Hansen) 、瑞典桩基委员会和 ISSMFE 提案也都曾假定压载试验的荷载 - 位移 (P-S) 曲线为双曲线 [3] 。曹汉志 [4] 通过试桩发现实测到的荷载传递曲线可近似用双曲线来描述。王幼青、张克绪 [5] 等人通过分析 71 根桩的压载试验的荷载 - 位移 (P-S) 曲线，得到 S/P-S 的线性回归的相关系数的平均值为 0.9976 ，这表明桩的荷载 - 位移 (P-S) 曲线完全可近似用双曲线关系来拟和。但该文中不分桩侧土体、桩端土体，均采用双曲线模型来模拟，模型中参数完全基于桩的静载荷试验值，不易推广。
由上述土力学理论及室内、室外试验结果，都表明在静荷载作用下桩周土体表现出非线性特性，并可用双曲线来描述荷载与位移的关系。但基于桩侧土体和桩底土体的变形及破坏机制不一样，而且桩端土不能承受拉力的特点，因此，桩侧与桩端土体静力模型 应用 不同的模型来描述。为简化起见，文中桩侧土体静摩阻力与剪切位移的模型采用双曲线关系，桩端土体仍采用理想弹塑性模型来描述 ( 即同 Smith 法 [6] ) 。
2 桩周土体模型
2.1 改进的桩侧土体模型 在动力打桩过程中，桩侧土体单元 i 在时刻 t 时所发挥的静阻力和动阻力分别由非线性弹簧 ( 双曲线 ) 和缓冲壶组成的模型来模拟 ( 如图 1) 。
2.1.1 桩侧土体静摩阻力双曲线模型 桩侧土体单元 i 在时刻 t 时所发挥的静摩阻力 R s (i,t)( 下文简写为 R s ) 与剪切位移 S(i,t)( 下文简写为 S) 假定用康纳 (Kondner) 双曲线拟和，即 R s -S 的关系曲线，如图 2 中 (a) 可用下式表示： R s =S/(a+bS) (1) 式中： a 、 b 为待定系数。
式 (1) 可改写为 S/R s =a+bS (2) 设以 S/R s 为纵轴， S 为横轴，可构成 S-S/R s 的新坐标系，则双曲线可转换 成直线，如图 2 中 (b) 其截矩为 a, 斜率为 b 。由式 (2) 可得，当 S→∞ 时， (3) 即用极限摩阻阻力 R ult 表示 S→∞ 时的值，也即为 R s 的渐近值。实际上工程中 S 不允许趋向无穷大，当达到一定值后，土体单元就被认为已破坏了，破坏荷载为 R f ，它总是小于 R ult 。可令： F r =R f /R ult (4) 称 F r 为桩侧土体单元的破坏比，则有 F r ＜ 1.0 。文中的 F r 假定可参照土体 Duncan Chang 模型 [7] 的破坏比取值。 同样，由式 (2) 可得，当 S→0 时， (5) 而 (SR s ) S→0 是曲线 R s -S 的初始切线刚度，也是该曲线上最大的切线刚度，用 K max 表示 ( 见图 2) ，因此 a=1/K max (6) 可将 (S f ,R f ) 代入式 (1) ，得 R f =S f /a+bS f (7) 由式 (7) 得到 a=1/K max =S f (1-F r )/R ult F r (8) 关于 S f 取值，可参照文献 [8] ：黏性土剪切破坏的最大位移一般为 5 ～ 8mm ，无黏性土剪切破坏的最大位移一般为 8 ～ 12mm 。
本文根据 Masing 准则和 Pyke 对 Masing 模型的修正方法 [9] 来确定卸载与加载的荷载 R S 与位移 S 的关系曲线。从第一象限状态点 (S 1 ,R s1 ) 处卸载阶段 ( 如图 3 中曲线 2) ，土体单元的静摩阻力为 R s =R s1 +(nR ult K max (S-S1))/nR ult +K max (S-S1) (9) 式中： n=1-R s1 /R ult 。
从状态点 (S 1 ,R s1 ) 处加载阶段 ( 如图 3 中曲线 3) ，土体单元的静摩阻力为 R s =R s1 +nR ult K max (S-S1)/nR ult +K max (S-S1) (10) 式中： n=-1-R s1 /R ult 。
采用 Pyke 法能有效地避免了卸载和反向加载时土体承受的荷载超过其最终强度的现象。
2.1.2 桩侧土体的动阻力计算 假定桩侧土体单元 i 在时刻 t 时的动阻力 R d (i,t) 与上节中计算出的静阻力 R s (i,t) 和单元速度 V(i,t) 成正比，比例常数采用桩侧土体单元的经验阻尼系数 J(i) ，即 R d (i,t)=J(i)R s (i,t)V(i,t) (11) 则在所发挥的总阻力 R(i,t) 为 R(i,t)=R s (i,t)+R d (i,t) (12) 2.2 桩端土体模型 桩端土体单元用理想弹塑性模型来描述土体力学与变形特性，即弹性、塑性与动阻力分别由弹簧、摩擦键及缓冲壶组成的土体流变模型来模拟 ( 图 4) 。土体单元的加载弹性变形由线段 OC( 弹簧 ) 模拟，塑性变形由线段 CD( 摩擦键 ) 来模拟， DE 段表示桩 端土卸载弹性变形阶段。详细的模型描述参见 文献 [6] 。
3 一维动力打桩波动方程
1931 年， D.V.Isaacs 指出能量从桩锤传递到桩底不是简单的刚体撞击动力 问题 ，而是撞击应力波在桩身内的传播问题 [10] 。他将桩周土体阻力 R′ 引入古典的一维波动方程得到 (13) 式中： x 为桩截面的位置坐标； u 为 x 处桩截面的质点位移； t 为时间； R′ 为桩身土阻力； C 为弹性应力波波速，单位： m/s ； ； E 、 ρ 分别为桩材料弹性模量、 密度。
然而式 (13) 仅仅在古典的一维波动方程中机械而简单地引入桩周土体阻力 R′ ，并没有真正反映出打桩过程中波动响应。根据方程各个分量的量纲 分析 ，桩周土体阻力 R′ 的量纲应为 1/m ，而力的量纲为 N ，因此， R′ 的物理意义不明确。
考虑桩身自重和桩周土阻力的情况下来建立一维动力打桩波动方程如下： (14) 式中： R 为桩单元受到的土阻力，单位： N ； G 为桩单元的重量，单位： N ； g 为重力加速度，单位： m/s 2 。其它量同式 (13) 。方程右边的第二项为桩的自 重对应力波传播的 影响 ；方程右边的第三项为桩周土体阻力对应力波传播的影响。
显然，该方程与方程 (13) 明显不同，而且有： (15) 因此，严格地讲 R′ 不是桩周土阻力，而是反映桩周土体阻力的一个参数，其量纲为 1/m 。
目前 求解波动方程 (14) 常采用特征线法、有限元法及差分法。由于差分法简单且实用，故文中仍然采用该法，以一维应力波动 理论 为基础，采用上述桩周土体模型编制了动力沉桩分析 软件 ADP(Analysis of Driving Pile) 。
4 工程实例
运用 ADP 程序对某海洋石油平台动力打桩工程进行 计算 分析。在实际试桩过程中进行了沉桩过程中贯入度、应力等数据监测，并在打桩完成后一段时间进行了复打及静载荷试验。
4.1 工程介绍 某海洋石油平台地区水深 24.5m ，潮差 1.0 ～ 1.5m 。从地质资料由上而下基本可分 5 层 ( 表 1) ，土质不均匀，夹层较多，表面有较厚的淤泥层在覆盖。平台桩基试桩为钢管桩，贯入深度 0 ～ 16.1m 范围内用 M40 型柴油锤打桩，在贯入深度 16.1 ～ 26.48m 范围内采用 MB70 型柴油锤施工。该桩总长为 69.11m ，整个桩打入到设计贯入深度 26.48m ，即桩端持力层位于第 ⑤ 层 —— 粉细砂层，如图 5 所示。桩身壁厚变化，即 δ=14 ～ 36mm 。 表 1 某海洋平台试桩处的各土层的物理力学指标 土体名称 标高
/m 含水量
ω 饱和容重
γ/(t·m -2 ) 比重 液限
W L 塑限
W P 塑性指数
I P 孔隙比
e 粘聚力 C/kPa 内摩擦角
φ/° ① 淤泥
② 淤泥质粘土
③ 粉细砂
④ 亚黏土
⑤ 粉细砂 0.0--4.5
-4.5--5.5
-5.5--10.5
-10.5--16.5
-16.5--29.5 41.5-54.4
21.5-25.8
15.7-27.1 1.7-1.75
2.0-2.05
1.93-2.12 2.65-2.67
2.68-2.72
2.68-2.70 44.5-53.8
27.1-33.3
- 20.0-21.4
14.9-18.1
- 24.5-32.4
11.7-15.2
- 1.16-1.41
0.599-0.72
0.5-0.773
5.0-8.0
-
23.5-27.0
36.0-40.0 表 2 各单元的重量 G(i) 和横截面积 A(i) 单元编号 i 3-7 8-11 12-15 16-19 20-23 24-27 28-31 32-35 36-39 A(i)/cm 2
G(i)/kN 不包括土塞 853.
6.77 958.6
7.66 1017.5
8.14 1155.0
9.07 1155.0
9.07 1155.0
9.07 1155.0
9.07 1155.0
9.07 1060.8
8.46 单元编号 i 40-43 44-47 48-51 52-55 56-59 60-63 64-67 68-71 A(i)/cm 2
G(i)/kN 不包括土塞 595.6
4.75 568.0
4.46 568.0
4.46 568.0
4.46 568.0
4.46 568.0
4.46 568.0
4.46 568.0
4.46 注：桩身外露单元共 45 个，其中 46 ～ 49 号单元桩周土为表层淤泥，程序中不考虑该土层对桩阻力的贡献。 4.2 计算参数选取 分别采用 Smith 法模型 ( 桩侧土体和桩端土体均采用理想弹塑性模型 ) 和本文改进模型 ( 桩侧土体用双曲线模型，桩端土体用理想弹塑性模型 ) 来进行动力打桩数值计算分析。为简化起见，计算中仅考虑用 MB-70 型柴油锤进行分析。具体参数选取如下： (1)MB-70 型柴油锤：锤心重 72kN ，锤心落高为 2.7m ，锤心刚度为 2.1671×10 6 MPa ，锤击效率 η=0.95 ； (2) 锤垫为白棕绳，其弹簧常数为 5.0MN/cm ，恢复系数 e=0.5 ； (3) 冲击块和桩帽：冲击块重 30kN ，桩帽重 21kN ； (4) 无桩垫； (5) 钢管桩的几何参数详见图 5 、表 2( 不考虑土塞效应 ) ，桩材弹性模量 E=2.1×10. 5 MPa 。为提高分析精度及更好地考虑土体的性质，每个桩单元长度设为 1m ； (6) 据时间步长选取原则，计算时间步长 Δt=0.00018s ； (7) 桩周阻力分布假定：桩端土的阻力分配比 PER=40% ，且桩侧土体侧摩阻力均匀分布； (8) 土质参数：因桩是大直径管桩，据文献 [11] 可适当加大桩端土最大弹性量的取值。 Smith 法模型：桩侧土 Q S =0.254cm ，桩端土 Q P =0.3048cm ；改进模型：桩侧土体参数包括：土体破坏变形值 S f ，土体的破坏比 F r ( 见表 3) ，桩端土体最大弹性变形 Q P =0.3048cm 。两模型的土体阻尼系数均取：桩侧土 J S =0.22s/m ，桩端土 J P =0.66s/m 。 表 3 桩侧土体 计算 模型参数 土体名称 土体破坏变形 S f /mm 破坏比 F r ① 淤泥
② 淤泥质粘土
③ 粉细砂
④ 亚黏土
⑤ 粉细砂 5.0
5.0
8.0
5.0
8.0 0.8
0.85
0.85
0.9
0.85 4.3 桩的可打入性预测 分析 运用 ADP 程序对试桩进行可打入性分析。图 6 、图 7 及图 8 分别为两模型预测桩在打入过程中所需锤击数、所发挥的静阻力及桩身最大拉、压应力值；图 9 、图 10 分别为在桩打入到 26.48m 处在一次锤击过程中桩顶单元的应力响应曲线、桩身所有单元的最大应力值。由图 6 分析，在整个贯入过程中，改进模型得到的所需锤击数与实测贯入曲线非常接近。从泥面开始到贯入深度 10m 范围内，两模型得到的所需锤击数基本一致，但从 10 ～ 26.48m 范围内，改进模型的计算值小于 Smith 模型的相应结果。在最大贯入深度 26.48m 处，本文模型、 Smith 模型得到的贯入 30cm 锤击数分别为： 357 击、 681 击，前者约为后者的 50% ，而实测值为贯入 30cm303 击。因此，改进模型的预测结果较准确，相对误差为 18% 。从计算和实测结果发现，沉桩到最大贯入深度处时 ( 位于 ⑤ 粉细砂层 ) ，非常难打。 由图 7 分析，随着贯入深度增加，发挥的桩周静阻力基本线性增加，这与计算中假定的桩侧土体均匀分布有关。本文改进模型计算值略小于 Smith 法模型的相应值，在最大贯入深度处所发挥的桩周静阻力分别为： 9.6MN 和 8.699MN 。 由图 8 分析，随着贯入深度增加，桩身的最大拉应力反而减小。在泥面附近贯入深度处，拉应力值最大，两模型的计算值均约为 139.0MPa 。故在沉桩初期 ( 特别对于混凝土桩 ) ，应尽量采取措施减少锤击能量 ( 如轻锤，或重锤小落距 ) ，以减小锤击拉应力值，防止桩身拉裂。在整个沉桩过程中，两模型得到的最大压应力值变化趋势非常一致，即先增加后减小，但本文改进模型所得到的最大压应力值略小于 Smith 法模型的相应值，在最大贯入深度一半处，两模型的最大压应力值都达到最大，分别为： 183.2MPa 、 197.6MPa 。钢管桩的屈服强度 F y 为 360MPa ，显然在沉桩过程中桩不会发生柱状屈曲，即满足强度要求。
由图 9 分析，在 0 ～ 15ms 范围内两模型分析得到的在一次锤击过程中桩顶单元应力响应的时程曲线基本一致，在 t15ms 时，本文改进模型所得的应力值略小于 Smith 法。两模型得到 的最大值均为 120.1MPa ，略小于实测最大值 133.3MPa ，误差约为 10% ，满足精度要求。
由图 10 分析，在最大贯入深度时一次锤击过程中，两模型得到的桩身单元最大压应力包络图非常接近，且与实测最大压应力包络图也较接近，可用于打桩工程应力控制。其中 Smith 法模型、改进模型计算得到的最大值分别为 149.68MPa 、 146.65MPa ，发生的位置自桩顶 44m 处。而实测到的最大的最大值为 141.2MPa ，发生的位置自桩顶 7.19m 处。显然，最大应力值满足钢材允许应力值。
总之，本文改进模型所预测锤击数精度高于 Smith 法模型，更接近于实测结果，且得到的应力响应结果和实测结果吻合较好，可用于实际打桩工程可打入性预测分析。
4.4 桩的极限静承载力预测分析 用 ADP 程序对桩在设计贯入深度 26.48m 处进行了单桩极限承载力预测分析。图 11 为两模型预测桩的极限承载力曲线。由图 11 分析，在贯入 30cm 时击数在 0 ～ 20 击范围内，两模型预测的单桩极限静承载力基本一 致；当贯入 30cm 大于 20 击时，本文改进模型的预测值逐渐大于 Smith 法模型的相应值。复打贯入 30cm 实测锤击数为 300 击， Smith 法模型、本文改进模型的预测值分别为 8.861MN 、 9.442MN ，而根据静载试验得到单桩的极限承载力为 9.6MN 。两模型计算结果与试验值的相对误差分别为 7.7% 、 1.6% 。二者的计算结果均满足工程精度要求，但改进模型的预测值的精度优于 Smith 法模型，更接近于实测结果，可用于实际打桩工程中单桩极限静承载力值预测分析。
5 结 论
(1) 本文根据土性 理论 的非线性和非弹性的特性以及室内与现场试验结果，而且基于桩侧土体和桩端土体在竖向荷载作用情况下的变形和破坏机理不同，假定用双曲线模型模拟桩侧土体静阻力，桩端土静阻力采用理想弹塑性模型来模拟是比较合理的。 (2) 原波动方程中土阻力 R′ 的物理意义不明确，只是反映桩周土体阻力的一个参数。本文以所推导的应力波动方程为基础，开发了可用于打桩分析的 软件 ADP 。 (3) 用程序 ADP 对某海洋平台打桩工程进行了 Smith 法模型和本文改进模型的对比计算分析。与实测结果比较表明：采用改进模型得到的预测单桩极限静承载力值、贯入过程中的所需的锤击数的计算精度更高，更接近于实测结果，而两模型得到的应力响应结果相近，和实测结果吻合较好。这也说明，双曲线模型更接近于实际情况，可用于实际打桩工程的可打入性及单桩极限承载力预测分析。