基于有限元极限分析的下伏溶洞路堤承载力计算

2021-06-23陈言章赵明华

地基处理 2021年2期

陈言章，赵明华

（1.广东省建筑设计研究院有限公司，广东广州 510010；2.湖南大学岩土工程研究所，湖南长沙 410082）

0 引言

我国岩溶地貌分布广泛，大量在建的山区公路不可避免地会穿越岩溶发育区，由于岩溶地质的复杂多样性，公路路堤容易发生溶洞的塌陷破坏。因此，如何确定下伏溶洞路堤承载力具有重要的工程价值。

目前，众多学者对该课题开展了相关研究，可归结为试验研究、理论研究和数值分析3个方面。试验研究方面，Al-Tabbaa等[1]基于小比例模型试验对下伏空洞地基的稳定性进行了研究；Kiyosumi等[2]对坚硬地层中含多个矩形空洞的条形基础承载力开展了离心机试验；刘庭金等[3]基于室内模型试验研究了含矩形空洞地基的破坏过程。当通过试验手段获得数据后，有必要对模型做适当的假定和简化，进一步对数据进行分析，提出符合试验规律的理论方法，并结合数值分析方法进行验证。理论研究方面，Wang等[4]假定作用在圆形空洞上条形基础的破坏模式，采用极限分析的方法得到了基础的极限承载力公式；刘辉等[5]利用极限分析上限法，建立了与假定破坏模式对应的速度场，求得了空洞上方条形基础极限承载力；刘之葵等[6]基于弹性理论求得溶洞周围岩体的应力状态，并结合格里菲斯强度准则，分析了含溶洞岩石地基的稳定性。数值分析方面，Azam等[7]利用二维有限元软件对含空洞时的地基承载力进行了数值计算，讨论了空洞顶板厚度和空洞位置等因素对地基承载力的影响；彭芳乐等[8]利用PLAXIS分析了单溶洞的存在对基础承载力和沉降的影响，对其发生机制进行了研究；在此基础上，Kiyosumi等[9]进一步考虑了条形基础作用在含多个溶洞地层时的极限承载力。阳军生等[10]采用ABAQUS软件对岩溶地基圆形基础作用下溶洞顶板的稳定性进行了分析，探讨了极限承载力与顶板跨度、顶板厚度等影响因素的关系。以上有限元方法虽取得了一定的成果，但采用传统有限元方法研究承载力问题尚存在计算结果容易发散、效率较低等问题，有限元极限分析法通过降低强度或增加荷载使岩体最终达到极限破坏状态，自动生成破坏面，同时可得到极限荷载，克服了传统有限元方法需根据位移-荷载曲线确定极限承载力的不足。其基本思想是利用有限元将应力场离散化，然后在离散的应力场内按照上、下限定理的相关要求构建相应的数学规划模型，最后选用合适的数学规划算法求解该模型，搜索出极限应力场和上、下限荷载[11-12]。该方法不需要人为构造静力容许的应力场、机动许可的速度场，对岩溶区路堤极限承载力的研究较为适合。

鉴于此，本文采用修正的Hoek-Brown准则对下伏溶洞路堤承载力进行计算，重点分析溶洞形状、分布形态、溶洞位置、岩石物理力学参数对路堤承载力的影响，并给出极限破坏模式，最后，将作用在岩体上的条形基础承载力与本文计算结果进行对比，以验证程序的正确性及路堤承载力结果的可靠性。

1 计算原理与方法

本文将采用有限元上、下限分析方法研究下伏溶洞路堤承载力问题，下面简单介绍其基本原理和计算程序求解过程。

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

1.2 非线性规划模型的建立

单元离散后，根据上、下限定理，建立节点应力和节点速度的约束方程，确定合理的优化目标函数。下面将分别给出上、下限分析数学规划模型的具体形式。

上限分析数学规划模型具体形式为[11]：

下限分析数学规划模型具体形式为[12]：

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

1.3 有限元上、下限分析的计算机实现

有限元上、下限分析的求解过程如图2所示，本文以 MATLAB为平台编制相关计算机程序，计算网格的划分，对优化模型建构和求解，并调用Tecplot360软件实现数据信息的可视化。具体过程的论述可参考文献[13-15]。

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

2 问题的描述

2.1 基本假定

廖丽萍等[16]对地基中椭球型空洞稳定性进行了分析，分析结果表明：在远场应力状态相同的条件下，椭球洞比椭圆孔更稳定。戴自航等[17]基于ABAQUS对室内试验进行了数值模拟，研究表明：相同条件下，三维椭球状溶洞的稳定性要大于二维椭圆形溶洞的稳定性，而后者的稳定性又高于同等条件下矩形溶洞的稳定性。鉴于此，本文将公路路堤下伏溶洞顶板的承载力问题简化为二维平面模型，计算模型如图3所示，并假定：（1）路堤荷载视为均布荷载，土层按层状分布；（2）溶洞截面形状为矩形，周边岩层为均质材料，且满足修正的Hoek-Brown准则。

图3 计算模型Fig.3 Calculation model

图3中：q为路堤荷载集度；d为路堤荷载的宽度；h1、hk、… hn分别为第1、k、…n层土的高度；γ为岩石的重度；γ1、γk、…γn分别为第1、k、…n层土的重度；X为路堤荷载中心与溶洞中心的水平距离；Y为岩层顶面与溶洞中心的垂直距离；W、H分别为溶洞跨度和高度；θ为溶洞的旋转角度，以逆时针方向为正。

2.2 数值分析模型

本文将路堤荷载作为极限荷载qu，上覆土层以均布荷载qs的形式作用于岩层上，其表达式如式（11）所示。

数值模型与网格划分如图4所示，路堤荷载宽度d=10 m，为减小边界条件的影响，取分析域宽为4 d，高为3 d，约束模型左右边界的法向位移，模型底部则进行完全约束，网格采用自适应划分技术[14]，单元总数为14 000个，以剪切耗散作为控制变量，初始单元数取1 000个，对网格进行3次自适应迭代。由图4可见，在能量耗散剧烈的区域，单元划分的较多，而在能量耗散较小的区域单元分布较为稀疏，在单元总数一定的条件下，使得计算域内塑性变形越剧烈的区域可获得越大的变形自由度，从而有效地降低数值离散误差。

图4 数值模型及网格划分Fig.4 Numerical model and mesh

2.3 材料参数及评价指标

为了描述岩石固有的非线性破坏特点，采用修正的 Hoek-Brown准则[18]，其表达式如式（12）所示。

式中：mi为与岩石完整程度有关的常数；D为扰动参数，现场无扰动岩体为0，而完全扰动岩体为1，本文取D=0。

由于Hoek-Brown准则是非线性准则，在屈服函数迭代计算过程中，所用的方法不同于其它线性屈服准则，在优化模型的求解过程中，除了需要计算当前迭代点的屈服函数值，还必须获得屈服函数对应力变量的一阶、二阶导数，即屈服函数的梯度向量和Hessian矩阵，具体过程可参考文献[15]。

采用Hoek-Brown准则的关键在于确定岩体参数GSI和岩石参数mi的范围。对于岩体参数GSI，Sonmez等[19]考虑不连续面的分布率、风化程度、粗糙度和填充物性质等影响，提出了详细的GSI定量评价方法。当GSI取值在40～90之间时，能基本满足实际工程的要求。对于岩石参数 mi，Hoek等[20]结合室内试验和工程经验，提出了较为全面的mi取值方法。鉴于岩溶主要成分为碳酸盐，其岩石参数mi的取值见表1，根据表1可知岩溶区的岩石参数mi变化范围为6～15。

表1 碳酸盐类岩石参数mi的取值Table 1 Values of parameter mifor carbonate rock

Hoek 等[21]和 Vá sá rhelyi[22]基于参数 GSI和 mi分别提出了计算岩石弹性模量E和泊松比v的经验公式，其表达式为：

3 路堤极限承载力

影响路堤极限承载力的因素主要有：（1）上覆土层荷载qs；（2）溶洞宽度W；（3）溶洞高度H；（4）溶洞旋转角度 θ；（5）路堤荷载中心与溶洞中心的水平距离X；（6）岩层顶面与溶洞中心的垂直距离Y；（7）岩石物理力学参数：GSI，mi等。由于本文采用了网格自适应划分技术，所得下限解和上限解相对误差在10%以内，为了使图表便于分析，本文取下限解和上限解的平均值作为路堤的极限承载力。下面将分别讨论各因素对路堤极限承载力的影响。

3.1 上覆土层荷载qs对极限承载力qu的影响

取H=0.1 d、θ=0°、X=0、Y=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。qs与 qu的关系如图5所示。

由图5可知，当W≤0.1 d时，qu随着qs增大而增大，当qs＞100 kPa后，qu的增长幅度逐渐趋于缓慢；当W≥0.2 d时，qu基本保持不变。因此，当溶洞跨度较大时，上覆土层荷载qs对路堤极限承载力影响不大，在实际设计时可不考虑。值得注意的是，在无溶洞时，qu与qs大致呈线性关系。

图5 qs对qu的影响Fig.5 Effect of qs on qu

3.2 溶洞跨度W对极限承载力的影响

取 qs=100 kPa、H=0.1 d、θ=0°、X=0、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。W与 qu的关系如图6所示。

图6 W对qu的影响Fig.6 Effect of W on qu

由图6可知，当Y＜0.4 d时，qu随着W的增大而减小，并趋向于0，qu减小的幅度也趋于缓慢；特别地，当Y=0.1 d时，qu接近于0，为保证溶洞顶板的稳定性，应保证溶洞顶板的厚度不能过小；当Y≥0.4 d时，qu随着W的增大大致呈线性减小。

3.3 溶洞高度H对极限承载力的影响

取qs=100 kPa、Y-H/2（即溶洞顶板厚度）=0.1 d、θ=0°、X=0、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。H与qu的关系如图7所示。

由图7可知，当W≤0.1 d时，qu随着H的增大呈非线性减小；当W≥0.2 d时，qu随着H的增大基本保持不变，说明当溶洞跨度较大时，溶洞高度对路堤极限承载力结果影响不大，该结论与文献[10]所得结论一致。

图7 H对qu的影响Fig.7 Effect of H on qu

3.4 旋转角度θ对极限承载力的影响

取 qs=100 kPa、X=0、Y=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。θ与 qu的关系如图8 所示。

图8 θ对qu的影响Fig.8 Effect of θ on qu

由图8可知，当W≤0.1 d时，θ<30°，qu随θ的增大而减小；θ＞30°，qu随θ的增大，先减小后增大；θ=30°，qu取最小值。当W≥0.15 d时，θ≤60°，qu基本不变；θ＞60°，qu增大或减小，要根据溶洞跨度确定。由以上分析可知，当溶洞跨度较大时，且θ≤60°时，qu随θ的变化幅度并不大，在实际工程设计时可不考虑θ的影响。

3.5 水平距离X对极限承载力的影响

取qs=100 kPa、W=0.2 d、H=0.1 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa。X与qu的关系如图9所示。

由图9可知，由于路堤荷载宽度远大于溶洞跨度，因此，qu随着X的增大先保持不变而后逐渐增大，且Y值越小，增大的幅度越大。根据图5，没有溶洞存在时，qu=35 MPa，当X≥0.9 d时，qu可达到无溶洞条件下 qu的 80%以上。由以上分析可知，实际工程中，轴对称情况为最不利情况，但Y值较小时，X值的改变对qu的影响比较明显，在设计时应当考虑其有利的影响。

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

3.6 垂直距离Y对极限承载力的影响

取 qs=200 kPa、hr=1 d、W=2 d、H=2 d、θ=0°、GSI=50、mi=7、γ=27 kN/m3、σci=40 MPa。Y 与 qu的关系如图10所示。

图10 Y对qu的影响Fig.10 Effect of Y on qu

由图10可知，当W≤0.3 d时，qu随Y的增大而增加，Y＜0.35 d时增加速度较快，Y＞0.35 d时增加速度逐渐变缓。当W≥0.4 d时，qu随Y的增大而增加，Y＜0.25 d时增加速度较快，Y＞0.25 d时近似线性增加。由以上分析可知岩层顶面与溶洞中心的垂直距离对路堤承载力的影响较为明显，在实际工程中，应注意保证溶洞顶板具有足够的安全厚度。

3.7 岩石物理力学参数对极限承载力的影响

取 qs=100 kPa、W=0.2 d、H=0.1 d、θ=0°、X=0、Y=0.2 d。岩石各物理力学参数对qu的影响如表2、表3和图11所示。

由表2、表3可知，有溶洞与无溶洞条件下，γ对qu的影响均不大，qu随σci的增大而增大。因此，在设计时可不考虑 γ的影响。由图 11可知，qu随GSI的增大而增大，且增大的速度随GSI增大而增大，此外，mi越大，qu越大。

表2 不同 γ条件下 qu的大小（GSI=50、mi=10、σci=40 MPa）Table 2 Ultimate bearing capacity quwith different γ(GSI=50、mi=10、σci=40 MPa)

表3 不同σci条件下qu的大小（GSI=50、mi=10、γ=25 kN/m3）Table 3 Ultimate bearing capacity quwith different σci(GSI=50、mi=10、γ=25 kN/m3)

图11 GSI对 qu的影响 (γ=25 kN/m3、σci=40 MPa)Fig.11 Effect of GSI on qu(γ=25 kN/m3、σci=40 MPa)

4 极限破坏模式

溶洞的存在是影响岩层极限破坏模式的最主要因素，因此下面将主要从溶洞各参数对极限破坏模式的影响展开讨论。

4.1 溶洞宽度W对破坏模式的影响

取qs=100 kPa、H=0.1 d、θ=0°、X=0、Y=0.15 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同W条件下的极限破坏模式，如图12所示。当W≤0.1 d时，溶洞侧壁发生破坏，顶板发生冒落破坏；当W=0.2 d时，溶洞顶板发生冲切破坏，同时伴有冒落破坏；当W≥0.4 d时，溶洞顶板发生冲切破坏。可见，随着溶洞宽度的增加，溶洞的破坏模式从侧壁发生破坏向顶板发生冲切破坏转变，转变过程中伴有顶板的冒落破坏。

图12 不同W条件下极限破坏模式Fig.12 Failure mechanisms with different values of W

4.2 溶洞跨度H对破坏模式的影响

取qs=100 kPa、Y-H/2（即溶洞顶板厚度）=0.2 d、θ=0°、X=0、W=0.2 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同H条件下的极限破坏模式，如图13所示。溶洞顶板发生冲切破坏，H对破坏模式的类型并无影响，这也印证了H对路堤极限承载力结果影响不大的结论。

图13 不同H条件下极限破坏模式Fig.13 Failure mechanisms with different values of H

4.3 旋转角度θ对破坏模式的影响

取qs=100 kPa、X=0、Y=0.3 d、W=0.3 d、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同θ条件下的极限破坏模式，如图14所示。破坏面由侧壁延伸至岩层顶面，且随着θ的变化而旋转。

图14 不同θ条件下极限破坏模式Fig.14 Failure mechanisms with different values of θ

4.4 水平距离X对破坏模式的影响

取qs=100 kPa、Y=0.15 d、W=0.2 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同X条件下的极限破坏模式，如图15所示。当X≤6时，溶洞顶板发生冲切破坏；当X≥8时，由于路堤荷载远大于上覆土层荷载，溶洞顶板和靠近路堤荷载的侧壁将发生破坏。

图15 不同X条件下极限破坏模式Fig.15 Failure mechanisms with different values of X

4.5 垂直距离Y对破坏模式的影响

取 qs=100 kPa、X=0、W=0.2 d、θ=0°、GSI=50、mi=10、γ=25 kN/m3、σci=40 MPa，得到不同Y条件下的极限破坏模式，如图16所示。当Y≤1.5 d时，溶洞顶板发生冲切破坏；当2 d≤Y＜3 d时，溶洞侧壁发生破坏，破坏面延伸到岩层顶面；当Y≥3 d时，溶洞侧壁发生破坏。可见，随着岩层顶面至溶洞中心垂直距离的增加，溶洞的破坏模式从顶板发生冲切破坏向侧壁发生破坏转变。

图16 不同Y条件下极限破坏模式Fig.16 Failure mechanisms with different values of Y

5 结果验证

Serrano等[23]假定岩体的重度 γ=0，基于修正的Hoek-Brown准则提出了一种计算条形基础极限承载力的方法。在此基础上，Merifield等[24]利用极限理论对岩石地基的极限承载力进行了数值模拟。为了验证本文方法的正确性，考虑无溶洞条件下条形基础作用在岩层上的极限承载力，将本文计算结果与文献[23-24]的研究成果进行对比，对比结果如表4所示。

由表4可知，本文所得条形基础极限承载力与文献[23-24]的结果基本一致，误差在 5%以内，因此本文所提方法较为可靠。