CN112417544A - 一种横向载荷作用下超长桩在成层地基中的水平位移求解方法 - Google Patents

Abstract

本发明公开了一种横向载荷作用下超长桩在成层地基中的水平位移求解方法，忽略了土体沿深度方向的连续性。随着水平载荷的逐渐增大，土体的连续性对成层地基中的水平位移影响很大，对p‑y曲线法的结果进行修正，弥补其在考虑土体的纵向连续性时的不足。Poulos方法对水平位移影响系数的求解是建立在匀质土Mindlin解的基础上，且仅对应半空间体作用集中荷载情况。实际工程中土体，各土层间特性可能差异很大，本发明将各层土体视为各向同性弹性体，假设桩截面土体水平抗力平均分布于桩周，根据弹性地基理论及层状弹性体系理论，提出了以传递矩阵法为基础，同时适用于有限杆单元法的桩基水平位移影响系数矩阵的求解方法，并给出了桩周土体水平位移的矩阵表达式。

Description

一种横向载荷作用下超长桩在成层地基中的水平位移求解 方法

技术领域

本发明属于岩土工程技术领域，具体涉及一种横向载荷作用下超长桩在成层地基中的水 平位移求解方法。

背景技术

由于p-y曲线法假设桩周不同深度的p-y曲线互不相关，忽略了土体沿深度方向的连续性。 随着水平载荷的逐渐增大，必须对p-y曲线法的结果进行修正，弥补其在考虑土体的纵向连 续性时的不足。Poulos弹性分析理论提出了考虑土体纵向连续性的方法。但Poulos方法对水 平位移影响系数的求解是建立在匀质土Mindlin解的基础上，且仅对应半空间体作用集中荷载 的情况。实际工程中的土体，各土层间特性可能差异很大，故本发明将各层土体视为各向同 性弹性体，假设桩截面土体水平抗力平均分布于桩周，根据弹性地基理论及层状弹性体系理 论，提出了以传递矩阵法为基础，同时适用于有限杆单元法的桩基水平位移影响系数矩阵的 求解方法，给出了桩周土体水平位移的矩阵表达式。

发明内容

为了解决现有技术问题，本发明的目的在于克服已有技术存在的不足，提供一种横向载 荷作用下超长桩在成层地基中的水平位移求解方法，将各层土体视为各向同性弹性体，假设 桩截面土体水平抗力平均分布于桩周，根据弹性地基理论及层状弹性体系理论，提出了以传 递矩阵法为基础，同时适用于有限杆单元法的桩基水平位移影响系数矩阵的求解方法，并给 出了桩周土体水平位移的矩阵表达式。

为达到上述发明创造目的，本发明采用如下技术方案：

一种横向载荷作用下超长桩在成层地基中的水平位移求解方法，包括如下操作步骤：

基于应力及位移的三角级数表达式，桩周抗力的分布形式与Mindlin解中所采用的集中力 形式，或路面工程中所采用的圆形或双曲线分布模式差异较大，假设地基为无桩状态，桩身 水平抗力平均分布于桩周，为构建表达该分布形式的函数，需利用狄拉克函数；

a)狄拉克函数定义如下：

b)矩阵传递法传递式的推导：

水平载荷作用于层状土体表面的情况如层状土体表面作用水平载荷图所示，根据层间完 全接触条件，其中地基土总深度Z趋向无穷远，在Z＝Zm，根据层间接触条件有：

其中z_sn+，z_sn-分别表示相应土体深度z_sn处的上下表面，其余下标方式类同；

令

并取m₀＝1，则在第s_n+1层：

在第i层

上式中，

则有，

以上方程进行递推即为水平荷载作用于表面的多层体系矩阵传递法表达式：

c)定解条件的选定

由于在实际工程中，定义足够远处的应力、形变或位移等于零；对于超长桩问题，由于 桩端距桩顶较远，水平载荷影响极其微小，因而定义桩底端位置水平位移为0；此外，先进 行摩阻力的分布运算；对超长桩问题，主要是摩擦桩或摩擦端承桩，桩底端沉降值取为0或 固定常数，求解过程中，桩端的定解条件简化为：

U_k(ξ,Z₁)＝V_k(ξ,Z₁)＝ω_k(ξ,Z₁)＝0

即

X(ξ，Z₁)＝[σ_ck(ξ,Z₁) τ₁(ξ,Z₁) τ₂(ξ,Z₁) 0 0 0]^T

上式中z_l为桩的入土深度；

此外，表面Z＝0处：

根据三角函数变换关系，则可知任一非周堆成载荷仅取K＝1的一项，其他项为零，则上 式转化为如下关系：

其中，

根据以上说明及推导，则土体表面边界有，

由以上各式即求得任意点应力位移情况；

d)水平荷载作用于层状地基内的传递矩阵解法：

d-1)矩阵传递法传递式的推导:

水平载荷作用于层状土体内部的情况如层状土体表面作用水平载荷图所示；

当水平载荷作用点为z＝z_v，将整个层状地基体系分为两部分考虑；其中，第一部分为从 z_inf至z_v+，第二部分为z_v-至Z＝0，z_v+表示z_v深度向下计算的起始点，z_v-表示z_v深度向上计算 的起始点；

分段写出从Z＝0处传递至底部的矩阵传递式，

上式中，

为水平载荷作用点应力位移矩阵增向量；

d-2)定解条件

对矩阵传递式进行求解的过程中，对桩底的边界条件仍按上式取用；

表面Z＝0处，由于地表处不存在载荷，故边界条件为：

根据两个边界条件，并不能完全反应层状地基体系中的实际受力情况；在z＝z_v处，有：

非轴对称载荷仍仅取k＝1的一项，其他项为零，则上式转化为如下关系：

则

表示为，

e)位移影响系数矩阵求解

e-1)矩阵传递法求解中的数值方法

通过上式并利用相应的定解条件，即求得任意深度z处经Hankel积分变换后的位移及应 力向量表达式

为求得实际位移及应力向量X(ξ,0)需进行相应的Hankel积分逆变换， 其表达式如下：

由此可见，层状弹性体系下的应力，位移公式均为无穷积分公式，难以利用手工完成计 算，必须采用数值积分方法，通过计算机编程完成相应的求解工作；

积分上限：

以上各式计算涉及无穷积分问题，由于计算机无法精确计算0～∞的积分值，采用一选 定的有限积分上限x_s取代无穷极限，则积分区间变更为[0，x_s]；通过对应力及位移分量的被 积函数特性进行分析来确定恰当的积分上限x_s；采用桩周水平抗力平均分布的假设进行运算， 所得被积函数为两部分乘积；若令x＝ξr₀，则两部分分别表示为与指数函数有关的部分E_x(x)， 及与Bessel函数有关的部分J(x)，则：

根据已知的传递矩阵，其指数函数为

的形式；k为指数函数系数，|k|≥1；H为计 算深度，其值不小于桩长；x_s～∞积分值为余项值，由于Bessel函数为波动收敛函数，当指数 函数中的幂次不小于15的情况下能保证积分余项中指数函数部分数值不超过10^-8，并保证积 分余项对数值解的精度不产生影响，写成不等式形式，则为，

相应的积分上限取 为:

则

根据精度的具体要求，适当放大积分上限值；

Newton-Cotes公式及自适应步长Simpson法

对于积分，在[a,b]区间任取n+1个结点a≤x₀<x₁<…<x_n≤b，构造f(x)的Lagrange插值多项 式L_n(x)，则插值型求积公式为：

若求积结点为等距，即x_i＝a+ih(i＝0,1,…,n)，步长

此时求积公式为，

称为n阶Newton-Cotes求积公式，其中，系数

为Cotes系数；

当n＝1时，求积公式为梯形公式:

当n＝2时，求积公式为Simpson求积公式：

e-2)为求得上式的近似数值解，利用自适应步长的Simpson法；若以

为 积分的基本表达式，具体过程如下，

①从梯形公式出发，计算：

②利用把区间逐次二分的办法，将区间[a,b]二分，令区间长度

计 算

则Simpson求积公式为

该递推公式仅需计算分割点上的函数值，不再重复计算原有各点上的函数值；

③重复步骤②，直至相邻两次结果

和

满足给定精度ε：

由于传递矩阵法主要为矩阵运算，利用Matlab中与自适应步长Simpson法对应的矢量化 函数，即在满足精度的前提下保证程序运行速度；

e-3)位移系数矩阵形式及求解流程

假设地基为层状横观各向同性各向同性弹性体，桩周土体位移参考水平荷载作用下的 Mindlin解给出；有限杆单元法的计算过程中，将桩身划分为若干单元，桩身结点编号为1～n， 则i结点处由于桩土相互作用产生的土体水平位移U_i表示为：

U_i＝u₁+u₂

其中u₁为i结点等效水平抗力在该节点处产生的水平位移，u₂为桩身其它节点处等效水 平抗力引起的i节点处土体水平位移；该式又表示为：

上式中

表示作用于j点的单位等效水平抗力在i点引起的位移，即水平载荷影响系数矩 阵[u^e]_n×n中i行j列的元素，

桩周土体水平位移矢量{U}为：

{U}＝[u^e]{p}

桩身受水平载荷作用时，桩周水平抗力实际为连续分布；对分布载荷，有限单元法采用 下式进行计算结点的等效载荷，

其中单元内的抗力计算采用了连续分布弹簧的假设进行，在结点i的水平抗力系数计算 过程中已考虑了相邻上下两结点i-1，i+1范围内土体抗力的连续分布，但忽略了相应抗力对 更远处土体水平位移的影响，因而上式写为：

其中u₁为i结点及相邻结点i-1及i+1上的等效水平抗力产生的水平位移，按p-y曲线法 计算结果取用；对水平载荷影响系数矩阵做相应改变，

根据p-y曲线法等方法计算所得的土体水平位移的计算结果为：

则较高的荷载水平下，桩周土体水平位移矢量{U}表示为：

{U}＝[u_i]+[u^e]{p}

e-4)位移系数矩阵数值求解流程

以层状弹性体系理论为基础，利用传递矩阵法求解水平载荷影响系数矩阵[u^e]_n×n数值解；

具体步骤说明如下：

①导入有限元计算过程中的桩身单元及结点划分信息，并增加结点与所属土层的信息， 取初始计算结点号为i＝1；

②于i结点作用水平单位等效载荷，根据传递矩阵计算式计算各土层传递矩阵 [G(ξ，H_i)]，并利用求出表层

表达式，并利用数值积分求出泥面位置位移，并注意 以下关系：

③求解水平荷载影响系数矩阵第i列元素，即利用

求结点j处

表达式和 数值解，并将水平位移数值解作为水平荷载影响系数矩阵[u^e]的j行i列元素，注意求解过程 中需对j结点位置进行判定：

当z_j≤z_i，根据下式进行计算；

当z_j>z_i，根据下式进行计算；

④令i＝i+1，并重复步骤②～③直至i＝n为止，n为结点数；

⑤对所得水平荷载影响系数矩阵[u^e]进行修改，以便用于修正较高荷载水平下p-y曲线 改进有限杆单元法求得的桩周土位移。

本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点：

1.本发明忽略了土体沿深度方向的连续性；随着水平载荷的逐渐增大，土体的连续性对 成层地基中的水平位移影响很大，由此必须对p-y曲线法的结果进行修正，弥补其在考虑土 体的纵向连续性时的不足；Poulos弹性分析理论提出了考虑土体纵向连续性的方法，但Poulos 方法对水平位移影响系数的求解是建立在匀质土Mindlin解的基础上，且仅对应半空间体作 用集中荷载的情况；实际工程中的土体，各土层间特性可能差异很大，故本发明方法将各层 土体视为各向同性弹性体，假设桩截面土体水平抗力平均分布于桩周，根据弹性地基理论及 层状弹性体系理论，提出了以传递矩阵法为基础，同时适用于有限杆单元法的桩基水平位移 影响系数矩阵的求解方法，并给出了桩周土体水平位移的矩阵表达式；

2.本发明方法更加精确获取横向载荷作用下超长桩在成层地基中的水平位移，对沿岩土 工程理论和实践具有重要意义。

附图说明

图1为本发明中单桩水平面边界示意简图。

图2为本发明中所述的狄拉克函数定义。

图3为本发明中桩周抗力分布假设示意图。

图4为本发明中层状地基示意图。

图5为本发明中所述层状土体表面作用水平荷载简图。

图6为本发明中所述层状土体内部作用水平荷载简图。

图7为本发明中所述桩基模型。

图8为本发明中所述程序流程图。

图9为本发明中传递矩阵法与Mindlin解计算结果比较图，图中5m处作用单位载荷。

图10为本发明中上下层土体模量不同时计算结果比较图，图中5m处作用单位载荷。

图11为本发明中不同深度单位荷载作用对土体表面产生的位移影响系数。

具体实施方式

以下结合具体的实施例子对上述方案做进一步说明，本发明的优选实施例详述如下：

实施例一：

在本实施例中，参见图1-11，一种横向载荷作用下超长桩在成层地基中的水平位移求解 方法，包括如下操作步骤：

基于应力及位移的三角级数表达式，桩周抗力的分布形式与Mindlin解中所采用的集中力 形式，或路面工程中所采用的圆形或双曲线分布模式差异较大，假设地基为无桩状态，桩身 水平抗力平均分布于桩周，参见图1，为构建表达该分布形式的函数，需利用狄拉克函数， 参见图3；

a)狄拉克函数定义如下：

b)矩阵传递法传递式的推导：

水平载荷作用于层状土体表面的情况如层状土体表面作用水平载荷图5所示，根据层间 完全接触条件，其中地基土总深度Z趋向无穷远，在Z＝Zm，根据层间接触条件有：

其中z_sn+，z_sn-分别表示相应土体深度z_sn处的上下表面，其余下标方式类同； 令

并取m₀＝1，则在第s_n+1层：

在第i层

上式中，

则有，

以上方程进行递推即为水平荷载作用于表面的多层体系矩阵传递法表达式：

c)定解条件的选定

由于在实际工程中，定义足够远处的应力、形变或位移等于零；对于超长桩问题，由于 桩端距桩顶较远，水平载荷影响极其微小，因而定义桩底端位置水平位移为0；此外，先进 行摩阻力的分布运算；对超长桩问题，主要是摩擦桩或摩擦端承桩，桩底端沉降值取为0或 固定常数，求解过程中，桩端的定解条件简化为：

U_k(ξ,Z₁)＝V_k(ξ,Z₁)＝ω_k(ξ,Z₁)＝0

即

X(ξ，Z₁)＝[σ_ck(ξ,Z₁) τ₁(ξ,Z₁) τ₂(ξ,Z₁) 0 0 0]^T

上式中z_l为桩的入土深度；

此外，表面Z＝0处：

根据三角函数变换关系，则可知任一非周堆成载荷仅取K＝1的一项，其他项为零，则上 式转化为如下关系：

其中，

根据以上说明及推导，则土体表面边界有，

由以上各式即求得任意点应力位移情况；

d)水平荷载作用于层状地基内的传递矩阵解法：

d-1)矩阵传递法传递式的推导:

水平载荷作用于层状土体内部的情况如层状土体表面作用水平载荷图6所示；

当水平载荷作用点为z＝z_v，将整个层状地基体系分为两部分考虑；其中，第一部分为从 z_inf至z_v+，第二部分为z_v-至Z＝0，z_v+表示z_v深度向下计算的起始点，z_v-表示z_v深度向上计算 的起始点；

分段写出从Z＝0处传递至底部的矩阵传递式，

上式中，

为水平载荷作用点应力位移矩阵增向量；

d-2)定解条件

对矩阵传递式进行求解的过程中，对桩底的边界条件仍按上式取用；

表面Z＝0处，由于地表处不存在载荷，故边界条件为：

根据两个边界条件，并不能完全反应层状地基体系中的实际受力情况；在z＝z_v处，有：

非轴对称载荷仍仅取k＝1的一项，其他项为零，则上式转化为如下关系：

则

表示为，

e)位移影响系数矩阵求解

e-1)矩阵传递法求解中的数值方法

通过上式并利用相应的定解条件，即求得任意深度z处经Hankel积分变换后的位移及应 力向量表达式

为求得实际位移及应力向量X(ξ,0)需进行相应的Hankel积分逆变换， 其表达式如下：

由此可见，层状弹性体系下的应力，位移公式均为无穷积分公式，难以利用手工完成计 算，必须采用数值积分方法，通过计算机编程完成相应的求解工作；

积分上限：

以上各式计算涉及无穷积分问题，由于计算机无法精确计算0～∞的积分值，采用一选 定的有限积分上限x_s取代无穷极限，则积分区间变更为[0，x_s]；通过对应力及位移分量的被 积函数特性进行分析来确定恰当的积分上限x_s；采用桩周水平抗力平均分布的假设进行运算， 参见图3，所得被积函数为两部分乘积；若令x＝ξr₀，则两部分分别表示为与指数函数有关的 部分E_x(x)，及与Bessel函数有关的部分J(x)，则：

根据已知的传递矩阵，其指数函数为

的形式；k为指数函数系数，|k|≥1；H为计 算深度，其值不小于桩长；x_s～∞积分值为余项值，由于Bessel函数为波动收敛函数，当指数 函数中的幂次不小于15的情况下能保证积分余项中指数函数部分数值不超过10^-8，并保证积 分余项对数值解的精度不产生影响，写成不等式形式，则为，

相应的积分上限取 为:

则

根据精度的具体要求，适当放大积分上限值；

Newton-Cotes公式及自适应步长Simpson法

对于积分，在[a,b]区间任取n+1个结点a≤x₀<x₁<…<x_n≤b，构造f(x)的Lagrange插值多项 式L_n(x)，则插值型求积公式为：

若求积结点为等距，即x_i＝a+ih(i＝0,1,…,n)，步长

此时求积公式为，

称为n阶Newton-Cotes求积公式，其中，系数

为Cotes系数；

当n＝1时，求积公式为梯形公式:

当n＝2时，求积公式为Simpson求积公式：

e-2)为求得上式的近似数值解，利用自适应步长的Simpson法；若以

为 积分的基本表达式，具体过程如下，

①从梯形公式出发，计算：

②利用把区间逐次二分的办法，将区间[a,b]二分，令区间长度

计 算

则Simpson求积公式为

该递推公式仅需计算分割点上的函数值，不再重复计算原有各点上的函数值；

③重复步骤②，直至相邻两次结果

和

满足给定精度ε：

由于传递矩阵法主要为矩阵运算，利用Matlab中与自适应步长Simpson法对应的矢量化 函数，即在满足精度的前提下保证程序运行速度；

e-3)位移系数矩阵形式及求解流程

假设地基为层状横观各向同性各向同性弹性体，桩周土体位移参考水平荷载作用下的 Mindlin解给出；有限杆单元法的计算过程中，将桩身划分为若干单元；参见图7，桩身结点 编号为1～n，则i结点处由于桩土相互作用产生的土体水平位移U_i表示为：

U_i＝u₁+u₂

其中u₁为i结点等效水平抗力在该节点处产生的水平位移，u₂为桩身其它节点处等效水 平抗力引起的i节点处土体水平位移；该式又表示为：

上式中

表示作用于j点的单位等效水平抗力在i点引起的位移，即水平载荷影响系数矩 阵[u^e]_n×n中i行j列的元素，

桩周土体水平位移矢量{U}为：

{U}＝[u^e]{p}

桩身受水平载荷作用时，桩周水平抗力实际为连续分布；每两个连续单元受土体水平抗 力作用简图见图7，对分布载荷，有限单元法采用下式进行计算结点的等效载荷，

其中单元内的抗力计算采用了连续分布弹簧的假设进行，在结点i的水平抗力系数计算 过程中已考虑了相邻上下两结点i-1，i+1范围内土体抗力的连续分布，但忽略了相应抗力对 更远处土体水平位移的影响，因而上式写为：

其中u₁为i结点及相邻结点i-1及i+1上的等效水平抗力产生的水平位移，按p-y曲线法 计算结果取用；对水平载荷影响系数矩阵做相应改变，

根据p-y曲线法等方法计算所得的土体水平位移的计算结果为：

则较高的荷载水平下，桩周土体水平位移矢量{U}表示为：

{U}＝[u_i]+[u^e]{p}

e-4)位移系数矩阵数值求解流程

以层状弹性体系理论为基础，利用传递矩阵法求解水平载荷影响系数矩阵[u^e]_n×n数值解 的流程简图见图8；

具体步骤说明如下：

①导入有限元计算过程中的桩身单元及结点划分信息，并增加结点与所属土层的信息， 取初始计算结点号为i＝1；

②于i结点作用水平单位等效载荷，根据传递矩阵计算式计算各土层传递矩阵 [G(ξ，H_i)]，并利用求出表层

表达式，并利用数值积分求出泥面位置位移，并注意 以下关系：

③求解水平荷载影响系数矩阵第i列元素，即利用

求结点j处

表达式和 数值解，并将水平位移数值解作为水平荷载影响系数矩阵[u^e]的j行i列元素，注意求解过程 中需对j结点位置进行判定：

当z_j≤z_i，根据下式进行计算；

当z_j>z_i，根据下式进行计算；

④令i＝i+1，并重复步骤②～③直至i＝n为止，n为结点数；

⑤对所得水平荷载影响系数矩阵[u^e]进行修改，以便用于修正较高荷载水平下p-y曲线 改进有限杆单元法求得的桩周土位移。

本实施例方法将各层土体视为各向同性弹性体，假设桩截面土体水平抗力平均分布于桩 周，根据弹性地基理论及层状弹性体系理论，提出了以传递矩阵法为基础，同时适用于有限 杆单元法的桩基水平位移影响系数矩阵的求解方法，给出了桩周土体水平位移的矩阵表达式。

实施例二：

本实施例与实施例一基本相同，特别之处在于：

在本实施例中，一种横向载荷作用下超长桩在成层地基中的水平位移求解方法，该计算 方法包括以下步骤：

某码头单桩入土深度36m，直径D＝1m。桩周土为两层，表层厚度2m，第二层至无穷远处。两层土泊松比相同，为0.3。为使结果更为直观，模量取无量纲数E₁＝1且E₂≥E₁。分 别利用本发明方法及经典Mindlin公式，计算0-36m深度范围内的桩周土体水平位移影响系数，并进行比较。

式中，c为水平集中力与泥面距离，

为增加可比性，并避免水平位移影响系数出现趋向无穷大的情况，对传递矩阵法所得结 果按下式计算：

在Mindlin公式计算中，取利用下式计算的桩横截面周长上多点的水平位移影响系数的平 均值作为Mindlin解的数值计算结果。

当土体模量E₂＝E₁时，5m深度处作用单位载荷对应的水平位移影响系数计算结果见图9。 由于两层土体模量一致，简化成弹性半无限空间问题，可直接将传递矩阵法与Mindlin公式结 果进行比较。由图9可见，除荷载作用点及临近两点Mindlin公式的计算结果偏大外，其余各 点计算结果基本一致。根据下式形式，构建用于修正计算的水平位移影响系数矩阵时，将去 除以上三点数据，可见假设抗力沿桩横截面平均分布的矩阵传递法可用于求解弹性半无限体， 证明了本发明方法的正确性。

当土体模量E₂/E₁分别为1，3，5时，5m深度处作用单位载荷对应的水平位移影响系数 计算结果见图10。由图可见，除单位载荷作用点以外，水平位移影响系数求解结果是连续的， 进一步证明了本方法及所编制程序的正确性。此外，由图直观可见，由于E₂/E₁逐渐增大，其 计算所得水平位移影响系数结果差异明显。根据Mindlin公式的形式，当模量或等代模量变化 时，不同深度水平位移影响系数按固定比例变化，不能很好反映层状体系各层分布情况造成 的影响。随着模量的变化，分析E₂/E₁＝1及E₂/E₁＝3情况的差异，可见当模量比增大时，不同 深度的水平位移影响系数的变化并不相同，自土体表面至深度5m处，其变化率范围为 0.49～0.31，可见本方法能更好地体现层状土体的实际分布差异的影响。

当土体模量E₂/E₁分别为1，2，5时，沿深度不同位置作用单位载荷对应的泥面水平位移 影响系数见图11。当E₂/E₁＝1时，即两层土体性质完全相同时，相应曲线较为平滑，而E₂/E₁>1 时，接近土体分层处，曲线变化趋势的转变较为明显，随着土体模量比增大，转折区域进一 步扩大，这一变化恰好说明当土层间差异较大时，各土层力学及承载性能将受到临近土层影 响。当E₂/E₁＝1时，曲线较为光滑，由于其结果可近似看作Mindlin公式的计算结果，可见 Mindlin公式并不能很好反映这一影响。

在较大的水平载荷作用下，由于p-y曲线法不考虑土体沿深度方向的连续性，计算所得 桩周土体水平位移误差可能需利用水平位移影响系数矩阵进行修正。针对该矩阵的求法，本 发明所做的工作及得出的结论如下：

1.提出假设抗力沿桩横截面平均分布，按弹性地基及层状弹性体系理论，以传递矩阵法 求解桩基水平位移影响系数矩阵的方法，克服了Mindlin公式应用于层状地基时的缺陷与不 便。并给出了修正较高荷载水平下p-y曲线法计算结果的土体水平位移矩阵式，以及适用于 该式的水平位移影响矩阵[u^e]的形式。

2.基于本发明的桩基水平位移影响系数矩阵的求解方法，编制了相应的Matlab程序并对 一双层地基土算例取土体模量E₂＝E₁进行程序及计算模型验证，可见：本方法计算结果同 Mindlin公式的计算结果基本一致，本方法是正确的。

3.根据Mindlin公式的形式，模量或等代模量变化时，不同深度水平位移影响系数按固 定比例变化，不能很好反映层状体系各层分布情况造成的影响。当E₂/E₁逐渐增大，本方法计 算所得水平位移影响系数结果差异明显，且不同深度的水平位移影响系数的变化并不服从固 定比例，可见本方法能更好地体现层状土体实际分布差异的影响。

4.当土体模量E₂/E₁分别为1，2，5时，沿深度不同位置作用单位载荷对泥面水平位移 的影响随深度的增大而减小。根据E₂/E₁＝1或Mindlin公式求解的结果绘制的曲线较为平滑。 当E₂/E₁>1，在接近土体分层处，曲线将出现明显的转折区，随着土体模量比增大，转折区域 进一步扩大，说明当土层存在差异时，本方法求解的结果能够反映临近土层间的相互影响。

本实施例横向载荷作用下超长桩在成层地基中的水平位移求解方法忽略了土体沿深度方 向的连续性。随着水平载荷的逐渐增大，土体的连续性对成层地基中的水平位移影响很大， 由此必须对p-y曲线法的结果进行修正，弥补其在考虑土体的纵向连续性时的不足。Poulos 弹性分析理论提出了考虑土体纵向连续性的方法。但Poulos方法对水平位移影响系数的求解 是建立在匀质土Mindlin解的基础上，且仅对应半空间体作用集中荷载的情况。实际工程中 的土体，各土层间特性可能差异很大，故本实施例方法将各层土体视为各向同性弹性体，假 设桩截面土体水平抗力平均分布于桩周，根据弹性地基理论及层状弹性体系理论，提出了以 传递矩阵法为基础，同时适用于有限杆单元法的桩基水平位移影响系数矩阵的求解方法，并 给出了桩周土体水平位移的矩阵表达式。

上面对本发明实施例结合附图进行了说明，但本发明不限于上述实施例，还可以根据本 发明的发明创造的目的做出多种变化，凡依据本发明技术方案的精神实质和原理下做的改变、 修饰、替代、组合或简化，均应为等效的置换方式，只要符合本发明的发明目的，只要不背 离本发明的技术原理和发明构思，都属于本发明的保护范围。

Claims (1)

1.一种横向载荷作用下超长桩在成层地基中的水平位移求解方法，其特征在于，包括如下操作步骤：

基于应力及位移的三角级数表达式，桩周抗力的分布形式与Mindlin解中所采用的集中力形式，或路面工程中所采用的圆形或双曲线分布模式差异较大，假设地基为无桩状态，桩身水平抗力平均分布于桩周，为构建表达该分布形式的函数，需利用狄拉克函数；

a)狄拉克函数定义如下：

b)矩阵传递法传递式的推导：

水平载荷作用于层状土体表面的情况如层状土体表面作用水平载荷图所示，根据层间完全接触条件，其中地基土总深度Z趋向无穷远，在Z＝Zm，根据层间接触条件有：

其中

分别表示相应土体深度z_sn处的上下表面，其余下标方式类同；

令

并取m₀＝1，则在第s_n+1层：

在第i层

上式中，

则有，

以上方程进行递推即为水平荷载作用于表面的多层体系矩阵传递法表达式：

c)定解条件的选定

由于在实际工程中，定义足够远处的应力、形变或位移等于零；对于超长桩问题，由于桩端距桩顶较远，水平载荷影响极其微小，因而定义桩底端位置水平位移为0；此外，先进行摩阻力的分布运算；对超长桩问题，主要是摩擦桩或摩擦端承桩，桩底端沉降值取为0或固定常数，求解过程中，桩端的定解条件简化为：

U_k(ξ,Z₁)＝V_k(ξ,Z₁)＝ω_k(ξ,Z₁)＝0

即

X(ξ，Z₁)＝[σ_ck(ξ,Z₁) τ₁(ξ,Z₁) τ₂(ξ,Z₁) 0 0 0]^T

上式中z_l为桩的入土深度；

此外，表面Z＝0处：

根据三角函数变换关系，则可知任一非周堆成载荷仅取K＝1的一项，其他项为零，则上式转化为如下关系：

其中，

根据以上说明及推导，则土体表面边界有，

由以上各式即求得任意点应力位移情况；

d)水平荷载作用于层状地基内的传递矩阵解法：

d-1)矩阵传递法传递式的推导:

水平载荷作用于层状土体内部的情况如层状土体表面作用水平载荷图所示；

当水平载荷作用点为z＝z_v，将整个层状地基体系分为两部分考虑；其中，第一部分为从z_inf至z_v+，第二部分为z_v-至Z＝0，z_v+表示z_v深度向下计算的起始点，z_v-表示z_v深度向上计算的起始点；

分段写出从Z＝0处传递至底部的矩阵传递式，

上式中，

为水平载荷作用点应力位移矩阵增向量；

d-2)定解条件

对矩阵传递式进行求解的过程中，对桩底的边界条件仍按上式取用；

表面Z＝0处，由于地表处不存在载荷，故边界条件为：

根据两个边界条件，并不能完全反应层状地基体系中的实际受力情况；在z＝z_v处，有：

非轴对称载荷仍仅取k＝1的一项，其他项为零，则上式转化为如下关系：

则

表示为，

e)位移影响系数矩阵求解

e-1)矩阵传递法求解中的数值方法

通过上式并利用相应的定解条件，即求得任意深度z处经Hankel积分变换后的位移及应力向量表达式

为求得实际位移及应力向量X(ξ,0)需进行相应的Hankel积分逆变换，其表达式如下：

由此可见，层状弹性体系下的应力，位移公式均为无穷积分公式，难以利用手工完成计算，必须采用数值积分方法，通过计算机编程完成相应的求解工作；

积分上限：

以上各式计算涉及无穷积分问题，由于计算机无法精确计算0～∞的积分值，采用一选定的有限积分上限x_s取代无穷极限，则积分区间变更为[0，x_s]；通过对应力及位移分量的被积函数特性进行分析来确定恰当的积分上限x_s；采用桩周水平抗力平均分布的假设进行运算，所得被积函数为两部分乘积；若令x＝ξr₀，则两部分分别表示为与指数函数有关的部分E_x(x)，及与Bessel函数有关的部分J(x)，则：

根据已知的传递矩阵，其指数函数为

的形式；k为指数函数系数，|k|≥1；H为计算深度，其值不小于桩长；x_s～∞积分值为余项值，由于Bessel函数为波动收敛函数，当指数函数中的幂次不小于15的情况下能保证积分余项中指数函数部分数值不超过10^-8，并保证积分余项对数值解的精度不产生影响，写成不等式形式，则为，

相应的积分上限取为:

则

根据精度的具体要求，适当放大积分上限值；

Newton-Cotes公式及自适应步长Simpson法

对于积分，在[a,b]区间任取n+1个结点a≤x₀<x₁<…<x_n≤b，构造f(x)的Lagrange插值多项式L_n(x)，则插值型求积公式为：

若求积结点为等距，即x_i＝a+ih(i＝0,1,…,n)，步长

此时求积公式为，

称为n阶Newton-Cotes求积公式，其中，系数

为Cotes系数；

当n＝1时，求积公式为梯形公式:

当n＝2时，求积公式为Simpson求积公式：

e-2)为求得上式的近似数值解，利用自适应步长的Simpson法；若以

为积分的基本表达式，具体过程如下，

①从梯形公式出发，计算：

②利用把区间逐次二分的办法，将区间[a,b]二分，令区间长度

计算

则Simpson求积公式为

该递推公式仅需计算分割点上的函数值，不再重复计算原有各点上的函数值；

③重复步骤②，直至相邻两次结果

和

满足给定精度ε：

由于传递矩阵法主要为矩阵运算，利用Matlab中与自适应步长Simpson法对应的矢量化函数，即在满足精度的前提下保证程序运行速度；

e-3)位移系数矩阵形式及求解流程

假设地基为层状横观各向同性各向同性弹性体，桩周土体位移参考水平荷载作用下的Mindlin解给出；有限杆单元法的计算过程中，将桩身划分为若干单元，桩身结点编号为1～n，则i结点处由于桩土相互作用产生的土体水平位移U_i表示为：

U_i＝u₁+u₂

其中u₁为i结点等效水平抗力在该节点处产生的水平位移，u₂为桩身其它节点处等效水平抗力引起的i节点处土体水平位移；该式又表示为：

上式中

表示作用于j点的单位等效水平抗力在i点引起的位移，即水平载荷影响系数矩阵[u^e]_n×n中i行j列的元素，

桩周土体水平位移矢量{U}为：

{U}＝[u^e]{p}

桩身受水平载荷作用时，桩周水平抗力实际为连续分布；对分布载荷，有限单元法采用下式进行计算结点的等效载荷，

其中单元内的抗力计算采用了连续分布弹簧的假设进行，在结点i的水平抗力系数计算过程中已考虑了相邻上下两结点i-1，i+1范围内土体抗力的连续分布，但忽略了相应抗力对更远处土体水平位移的影响，因而上式写为：

其中u₁为i结点及相邻结点i-1及i+1上的等效水平抗力产生的水平位移，按p-y曲线法计算结果取用；对水平载荷影响系数矩阵做相应改变，

根据p-y曲线法等方法计算所得的土体水平位移的计算结果为：

则较高的荷载水平下，桩周土体水平位移矢量{U}表示为：

{U}＝[u_i]+[u^e]{p}

e-4)位移系数矩阵数值求解流程

以层状弹性体系理论为基础，利用传递矩阵法求解水平载荷影响系数矩阵[u^e]_n×n数值解；

具体步骤说明如下：

①导入有限元计算过程中的桩身单元及结点划分信息，并增加结点与所属土层的信息，取初始计算结点号为i＝1；

②于i结点作用水平单位等效载荷，根据传递矩阵计算式计算各土层传递矩阵[G(ξ，H_i)]，并利用求出表层

表达式，并利用数值积分求出泥面位置位移，并注意以下关系：

③求解水平荷载影响系数矩阵第i列元素，即利用

求结点j处

表达式和数值解，并将水平位移数值解作为水平荷载影响系数矩阵[u^e]的j行i列元素，注意求解过程中需对j结点位置进行判定：

当z_j≤z_i，根据下式进行计算；

当z_j>z_i，根据下式进行计算；

④令i＝i+1，并重复步骤②～③直至i＝n为止，n为结点数；

⑤对所得水平荷载影响系数矩阵[u^e]进行修改，以便用于修正较高荷载水平下p-y曲线改进有限杆单元法求得的桩周土位移。

Images