0.   引　言

随着社会的不断发展，工程建设范围不断拓展，工程中经常会遇到各种复杂的地质问题。岩溶地质作为一种危害较大的工程地质条件，经常会出现在不同的工程建设中^[1−3]。桩基是岩溶区应用最为广泛的基础形式，因此，岩溶区桩基的研究对保证桩基及工程安全具有重要意义。

针对桩基下伏溶洞顶板稳定性问题，诸多科研人员应用理论、物理模型及数值模拟试验的手段开展了系统地研究。理论方面，赵明华等^[4−6]通过深入分析桩端岩层的变形特性及其破坏机理，综合考虑了桩端持力岩层的抗冲切、抗剪切和抗弯拉三种破坏形式，通过建立符合其工程特点的力学模型，进行应力分析计算，提出了一系列桩端岩层安全厚度的确定方法。蒋冲等^[7]根据岩溶区桩端溶洞顶板工程特点，引入突变理论，结合溶洞顶板稳定性因素的模糊性，建立岩溶区桩端溶洞顶板稳定性分析的突变评判新方法。物理模型试验方面，张慧乐等^{[8- 9]}基于室内相似模型试验，对岩溶区嵌岩桩的承载能力与破坏模式进行了研究，以顶板厚度、溶洞直径、溶洞位置、溶洞极半径、溶洞赤道半径等为变化因素，通过分析不同条件下岩溶区嵌岩桩的破坏模式，提出一系列安全可靠，便于工程应用的计算公式，对破坏模式进行了量化分析。李建东等^[10]通过分析桩顶荷载作用下荷载－沉降、桩身轴力、桩侧阻力、桩端土体塑性区发展变化规律。对超长单桩和群桩的承载特性进行研究，为大厚度黄土地区超长桩基的承载特性研究提供参考。数值模拟方面，黎斌等^[11]利用三维有限单元法对桩基下伏溶洞顶板进行应力分析计算，并采用多元线性回归方法得到相应的验算顶板安全厚度的关系式。鲁贤成等^[12]采用FLAC3D软件，根据丹霞枢纽互通区地质资料，建立概化地质模型，对该岩溶区大直径嵌岩状荷载传递特性进行数值计算，探讨了不同桩径，嵌岩深径比，桩土刚度比，桩岩刚度比条件下桩基荷载的传递规律。

以上研究主要考虑单溶洞情况，针对桩基含多溶洞的地基稳定性研究较少。黄明等^[13]通过ABAQUS建立串珠状溶洞地层的桩基模型，对桩基穿过串珠状溶洞时的荷载传递特性与岩层的破化型式进行了研究。梁冠亭等^[14]将下限分析法嵌入有限元模拟中，根据桩与溶洞间的相对关系，详细探讨了嵌岩深度，上覆土自重和溶洞直径对桩基承载力的影响程度。张乾青等^[15]基于现场物探和室内试验结果明确了溶洞的发育情况和岩土体的物理力学参数，采用有限元软件建立了14组数值计算模型，分别对不同溶洞半径、不同溶洞位置及不同桩径的工况进行了模拟，确定BoxLucas荷载传递模型的适用性，提出了穿越无充填溶洞单桩的承载特性计算方法。Liang等^[16]在抗滑桩嵌岩段下方和前方建立了两种不同形状的溶洞，通过水平加载试验，研究了椭圆桩桩底和桩前溶洞对抗滑桩嵌固端稳定性的影响。

工程中，岩溶区地基中往往含有多个溶洞，地基承载力影响因素较多，桩基承载机理复杂。显然，当前针对含多溶洞桩基－地基承载机理的研究落后于工程实践。因此，针对岩溶地基含多个溶洞的桩基承载机理，本文拟应用ABAQUS有限元数值分析软件，依托某岩溶桩基工程，建立含多溶洞的桩基－地基模型，分别探讨含单一溶洞、双溶洞、溶洞间距以及侧溶洞几何特征对桩基极限承载力的影响，以期为相似工程实践提供计算参考，进一步丰富岩溶区桩基承载机理的研究。

1.   模型设置分析

1.1   工程背景

本文以湖北省武汉市洪山区某桩基工程为参考案例，该场地岩溶属于武汉地区典型的覆盖性溶洞。其中下伏基岩为中风化石灰岩，较坚固且完整性良好。基岩面之上覆盖有20~40 m厚的土层，桩径1 m，桩长43.8 m，穿过上覆土层，嵌入岩层2 m^[17]。基桩参数及工程地质条件见图1。

图  1  桩基参数及工程地质条件

Figure  1.  Pile foundation parameters and engineering geological conditions

下载: 全尺寸图片 幻灯片

1.2   单桩极限承载力特征分析^[18−20]

工程实践表明，单桩可能出现屈曲破坏，整体剪切破坏或刺入破坏，不同的破坏形式，监测到荷载－沉降（Q-S）曲线特征不同。当桩周土体软弱，桩身侧向抵抗力弱，桩身将呈现屈曲破坏；当桩穿过抗剪强度较低的土层，到达强度较高的土层，且桩身长度不大时，桩底土体形成滑动面而出现整体剪切破坏。以上两种情况，Q-S曲线多为“急剧破坏”的陡降型。对于陡降型Q-S曲线，《规范》规定，取曲线发生明显陡降的起始点所对应的荷载值作为桩基极限承载力。当桩身入土深度较大，桩周土层抗剪强度均匀，桩在轴向荷载作用下出现刺入破坏，Q-S曲线多为缓变型。对于缓变型Q-S曲线，《规范》规定，宜根据桩顶总沉降量，取S=40 mm时对应的荷载值作为桩的极限承载力；对于桩径d不小于800 mm的桩，可取S =0.05 d对应的荷载值。轴向荷载作用下单桩破坏模式及Q-S如图2。

图  2  轴向荷载下桩基破坏模式

Figure  2.  Failure mode of pile foundation under axial load

下载: 全尺寸图片 幻灯片

本文模拟工况中，桩身入土深度较大且桩周土体抗剪强度较均匀，Q-S曲线多为“渐进破坏”的缓变型。从安全、经济两方面综合考虑，本文将取S=0.05，d=50 mm时对应的荷载值作为桩基极限承载力。

1.3   Mohr-Coulomb塑性本构模型

Mohr-Coulomb本构模型可以很好的反映岩土类材料抗压和抗拉强度的不对称性，所需岩土体参数少且其中c，ϕ值可以通过各种不同的常规试验测定。因此，较其他本构模型具有较好的可比性，在工程实践中有着重要的作用和地位，得到广泛的应用^{[21- 22]}。

摩尔–库伦本构模型属于理想的弹塑性模型，摩尔–库伦模型应符合摩尔–库伦屈服准则。

摩尔–库伦屈服准则表达式：

式中：c为岩土体黏聚力；φ为内摩擦角；τ_n为极限抗剪强度；σ_n为剪切面上法向应力。

Mohr-Coulomb屈服准则较Tresca屈服准则、Miese屈服准则及双剪屈服准则的屈服应力是最小的，具有较大的安全性。岩土材料由剪应力与垂直应力共同作用使粒子间克服摩擦产生相对滑移破坏，因此Mohr-Coulomb屈服准则对岩土体较为实用。

2.   模型介绍

2.1   假设条件

由于实际工程的复杂性，在模拟中很难将各种影响因素全部考虑在内，因此，为使建立的模型更具针对性，模型建立及计算满足以下假设条件^{[11, 23]}：

（1）溶洞周围岩体以连续介质考虑，视作Mohr-Coulomd材料，为各向同性的地下空间半无限体，且溶洞在天然状态下稳定；

（2）不考虑地下水的影响；

（3）溶洞内部无填充或填充物无承载特性；

（4）桩为圆柱形，分析时不考虑桩身破坏。

2.2   几何模型

依托上述工程背景，采用大尺寸模拟半无限空间体，选取土体半径为桩横截面半径的45倍，模型宽45 m，高60 m，其中上覆土层41.8 m，岩层18.2 m，包括2 m桩身嵌岩段建立几何模型，模型简化为图3。

图  3  模型工况基本信息示意图

Figure  3.  Schematic diagram of basic information of model

下载: 全尺寸图片 幻灯片

2.3   本构模型及物理力学参数的选取

本文模拟桩身采用线弹性模型，岩土体采用Mohr-Coulomd本构模型。

在数值模拟分析中，岩土体岩石力学参数的选取至关重要，参数选取的准确性直接影响计算结果的准确性及精度。对参考工况岩体参数进行简化，岩层及桩身力学参数见表1。

表  1  桩和岩土体物理力学参数

Table  1.  Mechanical parameters of pile and rock

岩土名称	层厚/m	密度ρ/kg·m−3	弹性模量E/MPa	泊松比µ	内聚c/kPa	内摩擦角 φ/°
粉质黏土	7.3	1 940	15	0.4	38	15
含碎石黏土	29.2	1 970	23	0.4	40	17
红黏土	5.3	1 960	27	0.4	35	12
中风化石灰岩	18.2	2 550	2e4	0.25	1e4	45
桩	43.8	2 500	3.55e5	0.2

下载: 导出CSV 

| 显示表格

2.4   边界条件及网格划分

模型左右边界进行水平约束，模型底部边界进行水平，竖直约束。采用ABAQUS中的主–从接触算法进行桩土接触模拟，选择主从面时应遵循以下原则：从面网格划分比主面更精细；当密度相近时，选择材料刚度较大的平面作为主面。ABAQUS中经常采用库仑摩擦模型、罚函数摩擦模型、Lagrange摩擦模型以及动力学摩擦模型等^[24]。本文采用适用大多数接触问题的罚函数摩擦模型，模型边界及约束信息详情如图4。

图  4  模型接触及约束示意图

Figure  4.  Schematic diagram of model contact and constraints

下载: 全尺寸图片 幻灯片

ABAQUS软件中的网格划分精度影响着计算结果的精度和有效性，本文采用四结点双线性平面应变的四边形单元进行网格划分，同时为了减小计算误差及缩短计算时间，桩土接触面与溶洞附近单元网格划分较细，远离部分相对稀疏。

2.5   计算分析步及模拟工况设置

应力计算大概分三步进行：（1）地应力平衡分析，ABAQUS有多种地应力平衡方式，本文采用自动地应力平衡法，生成初始应力场。由于接触的存在对地应力平衡的计算收敛会造成影响，因此地应力平衡分析中暂时移除桩体。（2）激活桩体，地应力平衡分析完成后，激活桩体，添加桩土，桩岩接触。（3）外荷载分析，完成接触定义后，对桩顶施加分级竖向荷载，模拟桩体受荷过程^[25]。

设置G1工况组和G2工况组，探讨溶洞形状对桩基承载力的影响。

G1工况组：桩端含单一溶洞，洞高h不变。在不同溶洞顶板厚度H的情况下，以溶洞形状及溶洞跨度l为变量，探讨溶洞顶板厚度H、溶洞形状及溶洞跨度l对桩基承载力的影响。模型基本工况如图5。G2工况组：桩端含双溶洞，保持桩端溶洞为椭圆溶洞，且洞高h、洞跨l不变，桩侧溶洞洞跨l'不变。在不同溶洞顶板厚度H工况下，以桩侧溶洞形状，洞高h'为变量，探讨溶洞顶板厚度H，侧溶洞形状，洞高h'对桩基承载力的影响。模型基本工况如图6，G1与G2工况详细参数见表2。

图  5  G1组模型基本信息示意图

Figure  5.  Schematic diagram of basic information on G1 model

下载: 全尺寸图片 幻灯片

图  6  G2组模型基本信息示意图

Figure  6.  Schematic diagram of basic information on G2 model

下载: 全尺寸图片 幻灯片

表  2  G1、G2组模拟工况参数表

Table  2.  Table of parameters for G1 and G2

溶洞数	分组	溶洞特征
		顶板厚度H	桩端溶洞					桩侧溶洞				溶洞间距D
			形状	洞高h	洞跨l	高跨比h/l		形状	洞高h'	洞跨l'	高跨比h'/l'
单	G1	d/2d/3d/4d/5d	矩形椭圆	3d	3d	1		/
					6d	1/2
					9d	1/3
					12d	1/4
					15d	1/5
双	G2	2d/3d/4d/5d	椭圆	3d	6d	1/2		矩形椭圆	3d	3d	1	d
									6d		2
									9d		3

下载: 导出CSV 

| 显示表格

3.   计算结果及分析

3.1   荷载-位移曲线分析

导出G1，G2组模拟工况下的荷载、位移数据，绘制成Q-S曲线（图7）。计算结果显示，G2组桩基极限承载力明显小于G1组，说明桩侧溶洞的出现对桩极限承载造成了一定的影响；2组工况下Q-S均可认定为缓变型曲线，与理论分析一致，结合单桩基极限承载力确定方法，选取S=0.05，d=50 mm时对应的荷载值作为桩基极限承载力。

图  7  G1、G2工况组Q-S曲线

Figure  7.  Q-S curves of G1 and G2

下载: 全尺寸图片 幻灯片

本文模拟工程地质参数来自于“武汉地区岩溶对桩基承载力影响数值模拟研究”，因此为验证本文模拟的正确性及可靠性，提取相同及相近的四种工况下文献[17]模拟结果与本文模拟结果进行对比分析。

根据文献[17]模拟结果与本文模拟结果对比图（图8），在相同工况下，两者的Q-S曲线基本重合；在相似工况下，两者模拟曲线的出入符合溶洞对桩基极限承载力影响的基本理论。

图  8  模拟结果与文献结果对比图

Figure  8.  Comparison between simulation results and literature

下载: 全尺寸图片 幻灯片

以沉降量S=50 mm为控制指标，提取出各种工况下极限承载力，并按工况组的划分进行对比分析。

3.2   溶洞形状对桩基极限承载力的影响

图9为G1组桩基极限承载力变化特征图。计算结果表明：相同工况下，圆形桩端溶洞的桩基承载力皆大于方形溶洞。图9（a）表明，保持桩端溶洞洞高h=3d不变，桩基极限承载力随着溶洞洞跨的增大而减小，变化幅度随着溶洞顶板厚度的增大而减小。图9（b）表明，桩基极限承载力随溶洞顶板厚度的增大而增大；随着溶洞顶板厚度的增大，溶洞几何要素的改变对极限承载力的影响逐渐减小。图9（c）表明，桩端下伏圆形溶洞与方形溶洞对桩基承载力影响的差别随着溶洞顶板厚度的增大而减小，当溶洞顶板厚度达到5d时，溶洞形状对桩基承载力的影响可忽略不计。

图  9  G1组桩基极限承载力

Figure  9.  Ultimate bearing capacity of G1 pile foundation

下载: 全尺寸图片 幻灯片

张慧乐等^[9]通过室内模型试验方法探讨了桩端单溶洞时桩基极限承载力随溶洞顶板，溶洞直径的变化规律，得出桩基极限承载力随溶洞顶板厚度的增大而增大，随溶洞直径的增大而减小，与本文模拟所得结论一致。

图10为G2组桩基极限承载力变化特征图。图10（a）表明，保持桩端下伏溶洞为椭圆，仅改变桩侧溶洞形状，椭圆形侧溶洞的桩基极限承载力大于矩形侧溶洞；保持侧溶洞洞跨不变，桩基极限承载力随着侧溶洞洞高的增大而减小，侧溶洞形状对桩基承载力的影响程度随着侧溶洞尺寸的增大而增大。图10（b）表明，不同溶洞顶板厚度下，相同工况下，侧溶洞形桩对桩基极限承载力不利影响的差值基本恒定。

图  10  G2组桩基极限承载力

Figure  10.  Ultimate bearing capacity of G2 pile foundation

下载: 全尺寸图片 幻灯片

综上所述，溶洞顶板厚度小于5d时，桩下溶洞对桩基极限承载力有着明显不利影响，当溶洞顶板厚度达到5d时，溶洞对桩基极限承载力的影响可忽略不计；相同工况下，矩形溶洞对桩基极限承载力的不利影响大于圆形溶洞。

在含溶洞地基桩基的施工建设中，当溶洞顶板厚度较小，需充分考虑溶洞形状对桩基极限承载力的不利影响，以保证桩基施工及使用的安全性。

3.3   侧溶洞洞高、洞跨对桩基极限承载力的影响

为探讨桩侧溶洞洞高、洞跨对桩基承载力的影响，根据G1、G2组计算结果，优化计算工况，设置G3工作组：桩下双溶洞，桩端溶洞为椭圆溶洞，溶洞间距为D=1d，不同溶洞顶板厚度H工况下，侧溶洞形状不变，以桩侧溶洞洞高h'，洞跨l'为变量，探讨溶洞顶板厚度H，侧溶洞洞高h'，洞跨l'对桩基承载力的影响。G3工况详细参数见表3，模型基本工况如图11。

表  3  G3组模拟工况参数表

Table  3.  Parameters for G3

溶洞数	分组	溶洞特征
		顶板 厚度 H	桩端溶洞					桩侧溶洞				溶洞 间距 D
			形状	洞高h	洞跨l	高跨比h/l		形状	洞高h'	洞跨l'	高跨比h'/l'
双溶洞	G3	2d/3d/4d/5d	椭圆	3d	6d	1/2		椭圆	3d	3d	1	d
										6d	1/2
										9d	1/3
									6d	3d	2
									9d		3

下载: 导出CSV 

| 显示表格

图  11  G3组模型基本信息示意图

Figure  11.  Schematic diagram of basic information on G3 model

下载: 全尺寸图片 幻灯片

图12为G3组Q-S曲线。结果表明，G3工况组桩极限承载力明显小于G1工况组，出现这种状况的原因有：第一，桩侧溶洞的存在会进一步破坏地基完整性，地基与溶洞相接点的应力集中现象更加突出；第二，桩基下伏溶洞与桩侧溶洞之间岩层的破坏也会导致桩极限承载力的降低。G3组Q-S曲线为缓变型曲线，仍选取选取S=0.05，d=50 mm时对应的荷载值作为桩基极限承载力。

图  12  G3工况组Q-S曲线

Figure  12.  Q-S curves of G3

下载: 全尺寸图片 幻灯片

图13为G3组桩基承载力变化特征图，图中横坐标0表示无侧溶洞。计算结果表明：桩端存在侧溶洞时，桩基承载力明显下降；保持侧溶洞洞高、洞跨其中一个不变，极限承载力都随着另一个的增大而减小；当溶洞顶板厚度为H=2d时，侧溶洞几何尺寸的改变的改变对极限承载力的影响较大，但随着溶洞顶板的增大，这种影响在逐渐减小。图13（c）表明，当溶洞顶板厚度小于等于3d时，侧溶洞洞高的改变对桩基极限承载力的影响远大于洞跨的改变对桩基极限承载力的影响；当溶洞H=5d时，侧溶洞洞跨改变对桩基承载力的影响超过洞高改变对桩基承载力的影响。

图  13  G3组桩基极限承载力

Figure  13.  Ultimate bearing capacity of G3 pile foundation

下载: 全尺寸图片 幻灯片

3.4   溶洞间距对桩基承载力的影响

为使溶洞间距对桩基承载力的影响区分更明显，选取最不利溶洞顶板工况探讨溶洞间距对桩基承载力的影响。计算结果表明，当溶洞顶板厚度H=2d时，溶洞对桩基极限承载力的影响最大。设计G4工况组为桩下双溶洞，保持溶洞顶板厚度H=2d，桩端溶洞与桩侧溶洞工况相同，以溶洞间距D，溶洞洞跨为变量，探讨双溶洞工况，溶洞间距D对桩基承载力的影响。G4工况详细参数见表4，模型基本工况如图14。

表  4  G4组模拟工况参数表

Table  4.  Parameters for G4

溶洞数

分组

溶洞特征

顶板 厚度 H

桩端溶洞

桩侧溶洞

溶洞 间距 D

形状

洞高h

洞跨l

高跨比h/l

形状

洞高h'

洞跨l'

高跨比h'/l'

G4

2d

椭圆

3d

3d

1

椭圆

3d

3d

1

d 2d 3d 4d 5d

6d

2

6d

2

9d

3

9d

3

矩形

3d

3d

1

矩形

3d

3d

1

6d

2

6d

2

9d

3

9d

3

下载: 导出CSV 

| 显示表格

图  14  G4组模型基本信息示意图

Figure  14.  Schematic diagram of basic information on G4 model

下载: 全尺寸图片 幻灯片

图15为G4组桩基极限承载力变化特征图。图15（a）计算结果显示，桩下溶洞皆为椭圆时的桩基极限承载力大于矩形溶洞；保持溶洞高跨比一定，桩基承载力随着溶洞间距的增大而增大。图15（b）计算结果显示，随着溶洞间距的增大，侧溶洞对桩基承载力的影响逐渐变小，当溶洞间距达到5d后，桩基极限承载力趋近于稳定，即侧溶洞对桩基极限承载力的影响可忽略；对比溶洞形状的影响可以发现，矩形溶洞下变量对极限承载力的影响较大。

图  15  G4组极限承载力变化特征

Figure  15.  Ultimate bearing capacity of G4 pile foundation

下载: 全尺寸图片 幻灯片

赵明华等^[26]根据相似原理设计室内模型试验，得到了顶板底部的应变曲线。其结果表明：当溶洞顶板厚度为桩径2倍时，距桩中心距为5d位置的顶板极限应变趋紧为0，即溶洞顶板对桩基荷载扩散传递的响应范围为桩径的5倍（图16），与本文模拟结果相符：当溶洞间距达到5d后，桩基极限承载力趋近于稳定，即侧溶洞对桩基极限承载力的影响可忽略。

图  16  溶洞顶板极限应变曲线

Figure  16.  Curve of ultimate strain of the roof

下载: 全尺寸图片 幻灯片

综上所述，桩端侧溶洞的出现极大地降低了桩基极限承载力。在相同工况下，随着溶洞顶板厚度的增加，桩基极限承载力受溶洞的影响越来越小；圆形溶洞参数的改变对桩基承载力的影响较小；随着溶洞间距的改变，侧溶洞对桩基极限承载力的影响逐渐减小，当溶洞间距达到5d时，侧溶洞对桩基极限承载力的影响可忽略。

4.   结论与探讨

本文通过ABAQUS数值模拟对多溶洞桩基承载力特征进行分析，探讨了含单一溶洞、双溶洞以及侧溶洞几何特征对桩基极限承载力的影响，得到如下主要结论：

（1）含溶洞地基–桩基中，当溶洞顶板厚度小于5倍桩径时，溶洞的存在对桩基极限承载力有着不利影响；相同工况下，方形溶洞对桩基极限承载力的不利影响始终大于圆形溶洞。

（2）含多溶洞地基–桩基中，矩形侧溶洞对桩基极限承载力的不利影响大于圆形侧溶洞，侧溶洞单一尺寸增大，侧溶洞形状对桩基极限承载力的影响降低。

（3）含多溶洞地基–桩基中，当两溶洞边界间距达到5倍桩径时，可仅考虑桩端溶洞对桩基极限承载力的影响，桩侧溶洞对桩基极限承载力的影响可以忽略。

本文采用数值模型计算，仅考虑了桩下存在双溶洞，且双溶洞形心位于同一高程面上的情形。因此，地基中含有多溶洞，且溶洞形态及溶洞空间相对关系更为复杂时，多溶洞共同作用对地基的影响是值得后续更深入探讨的问题。