岩溶区桥梁双桩基础有限元极限分析*

肖 尧 赵明华 张 锐 徐卓君 赵 衡

(①湖南大学岩土工程研究所 长沙 410082)(②中南大学土木工程学院 长沙 410075)

0 引 言

在岩溶区进行桥梁设计与施工时，桩基的极限承载力计算关系到工程设计的安全与经济。目前，众多学者对该问题进行了研究，研究手段主要包括试验研究、理论研究和数值分析等。试验研究方面，王革力(2002)、刘铁雄(2003)将溶洞视为梁板模型，探讨了溶洞顶板厚度、跨度对基桩极限承载力的影响； 而张慧乐等(2013)则进一步研究了溶洞位置、大小和形状等因素对基桩极限承载力的影响，江松等(2017)对岩溶桩基动力效应进行了试验研究。理论研究方面，刘之葵等(2003)基于弹性力学平面应变假定，通过判断圆形溶洞边界上关键点是否发生破坏，进而反算得到基桩极限承载力； 赵明华等(2004， 2016a)提出了溶洞顶板抗冲切、抗剪切、抗弯拉的验算方法，并以三者最小值作为基桩极限承载力； 韩红艳等(2012)对岩溶路基溶洞顶板稳定性进行了研究。此外，雷勇等(2014)、赵明华等(2017， 2018)采用极限分析的方法得到了基桩承载力计算公式。数值分析方面，黎斌等(2002)、阳军生等(2005)和黄明等(2017)采用有限元的方法，对岩溶区基桩极限承载力进行了系统研究，得出了一些有益的结论。以上研究取得了一系列的成果，但都是针对单桩方面的研究，对岩溶区桥梁双桩基础极限承载力的研究似尚未见报道。在实际工程中，双桩基础比单桩的情况更为常见，且受力更为复杂，其与单桩承载机理主要区别在于，双桩基础可能存在荷载叠加的效应，从而使整体承载力降低。因此，有必要针对该问题作进一步研究。

本文的目的在于，结合有限元和极限分析的优点，通过MATLAB编程求解岩溶区桥梁双桩基础极限承载力。该方法相较传统有限元和传统极限分析的优势在于：(1)不需要人为构造可静应力场(下限解)及机动许可的运动场(上限解)； (2)将直接给出极限承载力的下限和上限，不需要通过荷载-位移曲线来确定极限承载力，有效地克服了读数时的误差； (3)结果容易收敛，计算效率较高。为此，下文将首先对有限元极限分析的基本原理作简要介绍； 其次，为了描述岩体非线性特点，采用Hoek-Brown准则，并简要介绍将Hoek-Brown准则嵌入计算程序的方法； 然后，详细讨论溶洞大小、位置、桩距等因素对极限承载力的影响； 最后，计算无溶洞条件下单桩极限承载力，并与已有成果进行对比，验证本文方法的正确性。

1 有限元极限分析方法简介

有限元极限分析方法被广泛应用于岩土工程稳定性分析中(Sloan， 2013)，该方法根据极限分析基本理论，结合有限元的方法，利用计算机求解得到严格的上、下限解。该方法不需要人为构造可静应力场、机动许可的速度场，对岩溶区桥梁双桩基础承载力的研究尤为适合。为此，下文将对其基本原理作简要介绍。

1.1 单元离散

如图 1所示，采用线性三角形单元对计算域进行离散，其中，下限单元每个节点i有3个未知应力分量(σxi，σyi，τxyi)，每个单元共2个未知体积力分量，即he=(hx.hy)； 上限单元每个节点j有2个未知速度分量(uj.vj)，每个单元共3个未知应力分量，即σe=(σx，σy，τxy)。

图 1 三角形单元示意图Fig. 1 Sketch of triangular element

图 2 基于退化单元的速度间断线示意图Fig. 2 Sketch of velocity discontinuity based on degenerated triangular elements

为了克服线性三角形单元难以模拟复杂速度场分布的缺点，本文采用Krabbenhøft et al.(2005)提出的基于“退化”单元的速度间断线设置方法。如图 2所示，两个厚度为0的单元构成了“退化”单元的速度间断线，其有一对节点坐标重合。

1.2 数学规划模型建立

单元离散后，根据上、下限定理，建立节点应力和节点速度的约束方程，以总的内能耗散或外力荷载作为目标函数，并得到相应的数学规划模型。该过程在文献(Lyamin et al.， 2002, 2005)中已有详细介绍，本文不再赘述，下面将直接给出上、下限分析数学规划模型的具体形式。

上限分析数学规划模型具体形式为(Lyamin et al.， 2002)：

MinimizeQ=σTBu-cTu，

(1)

Subject toAu=b，

(2)

(3)

(4)

fi(σ)≤0，j∈Jσ，

(5)

(6)

u∈Rnu，σ∈Rnσ，λ∈RE

下限分析数学规划模型具体形式为(Lyamin et al.， 2005)：

MaximizeQ=cTx，

(7)

Subject toAx=b，

(8)

fi(x)≤0，j∈J

(9)

x∈Rn

(10)

式中，x为全局节点应力列向量；fj(x)为屈服准则或其他条件产生的不等式约束，其他符号意义同前。

1.3 有限元极限分析的计算机实现

有限元上、下限分析的求解过程如图 3所示，本文以MATLAB为平台编制了相关计算机程序，计算网格的划分，对优化模型建构和求解，并调用Tecplot360软件实现数据信息可视化。具体过程的论述可参考文献(张锐， 2015； 赵明华等， 2015， 2016b)。

图 3 有限元上、下限分析的计算机实现Fig. 3 Computer implementation of upper and low bound finite element method

1.4 基于Hoek-Brown准则的有限元极限分析

Hoek-Brown准则能较好地描述岩体非线性特征，在实际工程中得到了广泛应用，其最新的表达形式如下(Hoek et al.， 2007)：

σ′1=σ′3+σc(mbσ′3/σc+s)a

(11)

式中，σ′1和σ′3分别为最大、最小主应力；σc为岩石单轴饱和抗压强度；mb，s和a为与地质强度指标GSI有关的参数，其表达式为，

mb=miexp[(GSI-100)/(28-14D)]

(12)

s=exp[ (GSI-100)/(9-3D)]

(13)

a=1/2+(e-GSI152212/e-20/3)/6

(14)

式中，mi为与岩石完整程度有关的参数；D为扰动系数，没有扰动取0，完全扰动取1，本文中D=0。

针对本文所研究的问题，也采用Hoek-Brown准则，由于Hoek-Brown准则不同于其他线性屈服准则(如Mohr-Coulomb准则、Tresca准则等)，因此有必要对屈服函数迭代计算过程中所用的方法进行介绍。

在优化模型的求解过程中，除了需要计算当前迭代点的屈服函数值，还必须获得屈服函数对应力变量的一阶、二阶导数，即屈服函数的梯度向量和Hessian矩阵，具体过程如下：

为了便于推导，将Hoek-Brown准则转化为应力不变量表达的形式，

(15)

其中，I1为第一应力不变量，J2为第二偏应力不变量，θ为第三偏应力不变量J3相关的应力Lode角，参数β、χ及函数g(θ)、h(θ)的表达式如下：

g(θ)=-2cos(θ)

(16)

(17)

(18)

(19)

由于Hoek-Brown准则的应力空间包络面上存在奇异点，导致不能被求导，为了让Hoek-Brown准则适用于非线性规划算法，需对其进行“光滑化”处理。采用一个微小项ε对屈服函数进行“双曲型近似”处理，即采用式(20)对J2进行代替。

(20)

ε是一个与材料抗拉强度相关的常数，其值可由式(21)确定。

ε=min(δ，μρ|ρg(0)+(ρh(0)+χ)α=0)

(21)

其中，μ=10-1，δ=10-6，ρ为屈服面与纵轴对应的交点。

则双曲型近似Hoek-Brown准则可表示为：

(22)

根据式(22)求一阶、二阶偏导便可求得屈服函数的梯度向量和Hessian矩阵，由于篇幅所限，具体过程请参考文献(张锐， 2015)，本文不再赘述。

2 问题描述

2.1 计算模型及假定

Serrano et al.(2014)对嵌岩桩桩端极限承载力的研究表明：三维模型比二维模型承载力大约高1.3倍。此外，廖丽萍等(2010)研究表明：椭球形溶洞比椭圆形溶洞更稳定。以上研究表明，将岩溶区桩基承载力问题简化为平面模型是偏于安全，是有利于工程的。为了建立模型的同时突出关键因素的影响，本文将问题简化为二维平面模型，如图 4所示，并假定：(1)上部荷载全部由桩端持力岩层承担； (2)双桩基础由横系梁连接，上部荷载由双桩基础共同承担，且不考虑桩体自身变形的影响； (3)溶洞截面形状为圆形，周边岩层为均质材料，且符合修正的Hoek-Brown准则。

图 4 计算模型Fig. 4 Calculation model

图 4中：q为上部荷载，d为桩径，h1、hk、…hn分别为第1、k、…n层土的高度，γ为岩石的重度，γ1、γk、…γn分别为第1、k、…n层土的重度，hr为嵌岩深度，X、Y分别为溶洞中心与左桩桩端中心的水平距离和垂直距离，R为溶洞半径，L为左桩与右桩中轴线的水平距离。

图 5 网格划分示意图Fig. 5 Sketch of meshing

2.2 网格划分

数值模型与网格划分如图 5所示，分析域宽为40id，高为25id，模型左右边界及桩侧进行法向约束，底部进行完全约束，网格采用自适应划分技术(张锐， 2015； 赵明华等， 2016)，单元总数为6000个，进行5次网格优化，初始单元数取1000，均布极限荷载qu作用在横梁上，上覆土层以均布荷载qs的形式作用于岩层上，其表达式如式(23)所示。

(23)

桩侧与岩层的接触面光滑，桩端与岩层接触面完全粗糙。桩与横梁作为整体，视为无重度的刚体，岩层为均质材料，符合修正的Hoek-Brown准则，岩体的弹性模量E和泊松比ν分别根据式(24)、式(25)确定(Hoek et al.， 2006； Vasrhelyi， 2009)。

(24)

ν=-0.002GSI-0.003mi+0.457

(25)

2.3 评价指标

为了评价溶洞对双桩基础极限承载力的影响，定义一个无量纲参数Nσ，表达式如下：

Nσ=qu/σc

(26)

表 1 误差分析Table 1 Error analysis

由表 1 可知，上限解与下限解的误差在10%以内，以上限解与下限解的平均值作为设计指标，则其与真实解的相对误差在5%以内，为了后续便于分析，本文参数Nσ取上、下限解的平均值。

3 结果与讨论

岩溶区桥梁双桩基础整体的极限承载力影响因素主要包括：(1)上覆土层荷载qs； (2)嵌岩深度hr； (3)溶洞半径R； (4)左桩桩端中心与溶洞中心的水平距离X； (5)左桩桩端中心与溶洞中心的垂直距离Y； (6)桩的间距L； (7)岩石物理力学参数，GSI，mi，γ。下面将分别讨论各因素对嵌岩桩桩端极限承载力的影响

3.1 上覆土层的影响

图 6 qs对Nσ的影响Fig. 6 Effect of qson Nσ

由图 6可知，Nσ随着qs的增大而增大，且X越大，增长的幅度越大。当X=0时，qs对Nσ的影响基本可忽略不计。

3.2 嵌岩深度的影响

图 7 hr对Nσ的影响Fig. 7 Effect of hron Nσ

由图 7可知，嵌岩深度越大，桥梁双桩基础的承载力会得到一定提高，文献(赵明华等， 2016)也得出了相类似的结论。

3.3 溶洞半径R的影响

由图 8可知，Nσ随着溶洞半径R的增大而不断减小，且减小的幅度变缓，直至趋于某一稳定值。在文献(张慧乐等， 2013)中也得出了相类似的结论。

图 8 R对Nσ的影响Fig. 8 Effect of R on Nσ

3.4 水平距离X的影响

图 9 X对Nσ的影响Fig. 9 Effect of X on Nσ

图 10 不同X条件下双桩与单桩承载力系数对比Fig. 10 Comparison of bearing capacity factor for double-piles and single pile with different X

为了分析双桩与单桩承载性能的区别，将双桩所得结果平均分配给左桩和右桩，保持所有条件不变，取Y=2id，计算得到仅有左桩、右桩时的极限承载力系数，对比结果如图 10所示。根据图 10，当X≥0时，双桩基础的承载力大致等于仅有左桩、右桩时承载力的平均值； 当X=-1.5id时，双桩基础承载力大于仅有单桩时的承载力，主要原因在于发生了图 16a所示的整体剪切破坏。

3.5 垂直距离Y的影响

与图 10类似的，取X=0，得到不同Y条件下双桩与单桩承载力系数对比结果，如图 12所示。从图 12中可以看出，Y/d≥5时，双桩基础承载力系数大约为仅有单桩的90%。产生该现象的原因可能是双桩基础存在荷载叠加的效应，因此双桩基础相对于单桩的承载力降低了。综合图 10、图12的结果，建议工程设计时，双桩基础的承载力计算由离溶洞最近的单桩承载力控制，并乘以0.9的折减系数。

图 11 Y对Nσ的影响Fig. 11 Effect of Y on Nσ

图 12 不同Y条件下双桩与单桩承载力系数对比Fig. 12 Comparison of bearing capacity factor for double-piles and single pile with different Y

3.6 桩间距L的影响

在X确定的情况下，桩间距L的影响主要体现在右桩与溶洞的位置关系，当右桩离溶洞较远时，双桩基础的承载力也会较高。双桩基础与溶洞的位置关系在前面已论述，在此不再赘述。

3.7 地质强度指标GSI的影响

图 13 GSI对Nσ的影响Fig. 13 Effect of GSI on Nσ

表 2 γ对Nσ的影响Table 2 Effect of γ on Nσ

3.8 参数mi的影响

图 14 qs对Nσ的影响Fig. 14 Effect of qs on Nσ

3.9 岩体自重γ的影响

图 15 不同R的溶洞极限破坏模式Fig. 15 Failure patterns of the cave with different R

图 16 不同X的溶洞极限破坏模式Fig. 16 Failure patterns of the cave with different X

图 17 不同Y的溶洞极限破坏模式Fig. 17 Failure patterns of the cave with different Y

4 破坏模式

根据计算的结果，溶洞大小和位置是影响极限破坏模式的关键因素，将不同R、X、Y条件下溶洞的极限破坏模式总结出来，如图 15～图 17所示。图 15中，X=0、Y=2id，当R=1id时，左桩下方的溶洞顶板发生冲切破坏，同时右桩出现地基破坏模式； 随着R的增大，破坏模式由左桩控制，右桩对破坏模式影响不大，图 15c、图15d所示； 当R≥4时，发生整体剪切破坏，破坏面贯穿左、右桩的外侧。

图 16给出了不同X条件下溶洞的破坏模式，其中，R=2id、Y=2id，当-1.5d≤X≤0，发生整体剪切破坏； 随着X的增大，破坏模式逐渐变为由左桩控制的冲切破坏，同时右桩底部发生地基破坏模式。

图 17给出了不同Y条件下溶洞的破坏模式，其中，R=2id、X=0，溶洞与桩端竖直距离似乎对破坏模式影响不大，以整体剪切破坏为主，破坏面由左、右两桩桩侧延伸至溶洞表面，这可以用来解释Nσ与Y大致呈线性关系。此外，值得注意的是，当Y=1id时，溶洞顶部会出现冒落区。

5 结果验证

表 3 无溶洞条件下单桩桩端承载力系数NσTable 3 Bearing capacity factor Nσof pile tip without void

6 结 论

(1)考虑实际工程中岩溶区桥梁双桩基础的承载特点，根据极限分析上、下限定理，利用MATLAB平台编制了有限元极限分析计算程序，用Hoek-Brown准则来描述岩体的非线性特征，并通过“双曲线近似”处理，将其嵌入计算程序中。

(2)上覆土层荷载和嵌岩深度越大，双桩基础的承载力越大； 溶洞半径越大，极限承载力越小，逐渐趋向于0； 极限承载力随GSI的增大非线性增大，与mi大致呈线性关系； 当桩与溶洞距离较远时，岩石重度越大，承载力越高，但当桩与溶洞距离较近时，岩石重度对承载力影响可忽略不计。

(3)左桩与溶洞的水平距离X增大，承载力先增大后减小，在X=0时取最小值； 承载力随垂直距离Y的增大而大致线性增大。通过与单桩承载力的比较，在工程设计时建议，双桩基础的承载力计算由离溶洞最近的单桩承载力控制，并乘以0.9的折减系数。

(4)岩溶区桥梁双桩基础的极限破坏模式主要有3种： ①左桩下方的溶洞顶板发生冲切破坏，同时右桩出现地基破坏模式； ②由左桩控制的冲切破坏模式； ③整体剪切破坏，破坏面由左、右桩的外侧贯穿至溶洞表面。

(5)本文所有计算结果都进行了无量纲化处理，并通过计算无溶洞时桩端极限承载力，与理论方法、FLAC计算结果对比，验证了本文所提方法的正确性，可为同类工程提供参考。