Bearing Capacity of Karst Cave Roof under Pile Foundation Load Using Limit Analysis

Abstract

Reducing the computational workload by simplifying the analysis of karst foundations into a plane strain problem can yield significant advantages. Yet, such an approach fails in reproducing the engineering situation in a rigorous manner. In this regard, this paper proposes an upper-bound method that can effectively analyze the bearing capacity of three-dimensional karst foundations. This method is utilized to investigate the impact of various pile diameters, the ratio of roof thickness to pile diameter, and the ratio of cave width to pile diameter on the stability of karst foundations. The validity of this method is established through an illustrative example. The outcomes illustrate that when subjected to both tensile and compressive horizontal stresses if the maximum horizontal stress surpasses the tensile strength of the rock mass, the roof rock mass may suffer damage. Increasing the ratio of roof thickness to pile diameter can bring down the horizontal stress value. The stability factor is largely influenced by the ratio of roof thickness to pile diameter. The most prominent growth trend of the stability factor is observed when the ratio is less than 3. If the ratio of the roof thickness to pile diameter exceeds 3, the prediction of the bearing capacity estimation for the karst foundation in three-dimensional circumstances is more conservative than that in two-dimensional circumstances.

1. Introduction

The stability of karst foundations is a crucial area of study within the field of geotechnical engineering. Research has revealed that the roof’s thickness, cave size, cave span, and location of the karst cave are the primary factors that influence the ultimate bearing capacity of the karst foundation [1,2]. Therefore, conducting an analysis of the correlation between these various factors and the bearing capacity of the karst cave roof is of utmost importance in guaranteeing the safety of karst foundation engineering projects [3].

Currently, scholars within this field are primarily concerned with addressing the issue of determining the ultimate bearing capacity of karst foundations. Drawing upon the fundamental principles of mechanical analysis, Xie, P. has produced an analytical solution to ascertain the ultimate bearing capacity of the karst cave roof, of which the calculated outcomes closely correspond with those produced by numerical simulation [4,5]. However, more complex strata distribution will usually be encountered in practical projects. This causes situations where karst foundations cannot be simplified as axisymmetric models, which inevitably brings great difficulties in deriving analytical solutions for the bearing capacity of karst foundations. The establishment of a geotechnical structure’s constitutive model through the finite element method (FEM) allows for predicting the karst foundation’s deformation and stability, rendering it highly applicable in resolving intricate load and boundary issues [6]. However, the accuracy of the FEM calculation’s outcomes can be influenced by multiple factors, such as meshing, type of element, and convergence criteria. Furthermore, the imprecise approximation of the displacement function poses a challenge in accurately assessing whether the FEM results are predisposed towards safety or hazard; precisely, it is difficult to determine whether it is an upper-bound solution or a lower-bound solution [7,8]. According to the upper-bound method (UBM), the real solution of the ultimate load is located below the upper-bound solution. Therefore, the UBM can provide an upper-bound prediction for the bearing capacity of karst foundations. At present, certain investigations employing UBM have been conducted on the ultimate bearing capacity of karst foundations. However, these inquiries have omitted the impacts of three-dimensional spatial factors and have streamlined karst foundations to a two-dimensional plane strain problem [9,10,11,12,13,14]. Although this strategy has successfully curbed computational expenses, it unavoidably leads to computational inaccuracies in complex three-dimensional circumstances [14].

The key to extending the UBM of finite element limit analysis from two-dimensional problems to three-dimensional problems lies in treating the yield function and selecting calculation methods. In two-dimensional cases, the Mohr–Coulomb yield function can be linearized through inscribed polygons, thereby, reducing the difficulty of solving and improving convergence speed [15]. However, the calculation workload will significantly increase if this method is applied to three-dimensional cases. Thus, the nonlinear programming solution of three-dimensional problems is one of the difficulties to be overcome by the UBM of finite element limit analysis. Many scholars have done much work in this regard. For example, the two-stage quasi-Newton algorithm and alternating direction method of multipliers are effective methods for solving high degree of freedom nonlinear programming problems [16,17].

Thus far, research on the ultimate bearing capacity of a three-dimensional karst foundation has predominantly relied on experimental observations [18,19,20]. Moreover, earlier investigation indicates that the primary contributor to cave failure is tensile stress [21]. Despite these efforts, there remains a scarcity of studies investigating the limit analysis of three-dimensional karst foundations and the stress experienced by karst cave roofs. To enable a thorough evaluation of karst foundation stability, an effective approach is required that utilizes the UBM of finite element limit analysis to examine the ultimate bearing capacity of the three-dimensional karst foundation.

Considering the aforementioned, the present paper amalgamates the concept of finite element discretization and uses tetrahedral elements to discretize the foundation. Kinematically admissible velocity fields have been established to comply with plastic flow constraints, velocity boundary conditions, and external power equal to internal power constraints. An ultimate bearing capacity analysis of the karst foundation is conducted through a three-dimensional nonlinear mathematical programming model, taking the ultimate load as the objective function. As the UBM fails to derive the structure’s stress field and the high degree of freedom nonlinear optimization model proves challenging to solve, the primary nonlinear mathematical programming model is transformed into a set of Kuhn–Tucker conditions and accomplished using Python.

This study directly investigates the impact of varying pile diameters, roof thickness, and cave width on the bearing capacity of a three-dimensional karst foundation. Rather than simplifying the foundation as a two-dimensional plane strain problem, the stability factor of the karst foundation is determined by reducing the shear strength parameters of the rock mass. This methodology is well suited for geotechnical structures exhibiting shear failure [22,23,24]. Additionally, the distribution of horizontal stress on the roof is calculated across different scenarios, including four pile diameters, the ratio of cave width to pile diameter, and the ratio of the roof to pile diameter.

2. Theory and Hypothesis

The UBM solely takes into account the movement state of the structure, utilizing the speed of the structure as the primary variable to resolve the ultimate bearing capacity of the structure during the critical state. The ultimate load consistently surpasses the actual load, implying that the solution is relatively secure for the structure. Through the incorporation of finite elements, the UBM has achieved recognition as an effective method for geotechnical structural stability analysis. The fundamental steps of this approach are

• (1). Discreting the geotechnical structure with finite elements, then establishing the kinematically admissible velocity field of the structure. Each velocity field corresponds to an external load;

• (2). Calculating the equivalent strain rate and equivalent strain of the velocity field using the geometric equation;

• (3). Calculating the upper-bound power generated by each part, which includes the internal power of the element, the power on the velocity discontinuity, and the power generated by external loads;

• (4). Calculating the external load that minimizes the total power consumption by optimization principles, and the minimum external load in all velocity fields is the ultimate load of the geotechnical structure.

The process outlined above can be succinctly expressed as a mathematical programming problem. Moreover, we aim to introduce a solution strategy centered around the concept of strength reduction. In order to render the computation process practicable, this article’s calculation model relies upon the ensuing hypothesis

• (1). Solely the bearing capacity of the pile end is taken into account, whilst disregarding the pile’s lateral friction resistance;

• (2). The rock mass is deemed to be of the ideal rigid plastic variety;

• (3). The external load conforms to the requirements of proportional loading;

• (4). The cave is devoid of any filling materials;

• (5). The impact of the roof’s self-weight is not deemed pertinent.

3. Mathematical Model

The karst cave roof exhibits primarily punching failure under the load of the pile foundation [11,12]. Consequently, the shear strength parameters of the materials may experience sustained reduction throughout the calculation process, to drive the karst foundation into a critical state. In this paper, the ratio of shear strength parameters before and after reduction is defined as a stability factor k, specifically

(1)$k=\frac{c}{c^{\prime }}=\frac{\tan \phi }{\tan {\phi }^{\prime }}$

where $c$ and $c^{\prime }$ are the cohesion before and after strength reduction, respectively. $\phi$ and ${\phi }^{\prime }$ are the internal friction angle before and after strength reduction, respectively.

3.1. Elements Discretization

As mentioned above, the first step is to discretize the karst foundation using tetrahedral elements, the element is shown in Figure 1a, and in the kinematically admissible velocity field, velocity discontinuity on the common surface of adjacent elements is shown in Figure 1b. Each element has six stress components and three velocity components. Assuming a linear change in velocity within each element, the velocity of each element can be expressed by Equation (2).

(2)${\mathit{u}}^e={\displaystyle \sum _{i=1}^4{N}_i{\mathit{u}}^i}$

where $i$ is the serial number of the node, ${\mathit{u}}^i$ is the velocity of the node $i$, and $N_i$ is a shape function related to the node coordinates $(x,y,z)$, which is given by

(3)$N_i=\frac{1}{6{\mathrm{V}}^e}(a_i+b_{i}x+c_{i}y+d_{i}z)(i=1,2,3,4)$

where ${\mathrm{V}}^e$ is the volume of the element and $a_i,b_i,c_i,d_i$ are the coefficients related to the geometric position of the node $i$.

3.2. Yield Condition of Rock Mass

When the combination of shear stress and normal stress in any plane reaches a critical condition, the rock mass begins to yield. The relationship between shear and normal stress can usually be described using the Mohr–Coulomb strength criterion, which can be expressed as

(4)$f^e(\mathit{\sigma })=\frac{I_1}{3\sin \phi }-\sqrt{J_2}\left[\cos {\theta }_{\sigma }+\frac{1}{\sqrt{3}}\sin \phi \sin {\theta }_{\sigma }\right]+c\cos \phi \le 0;e=(1,\cdots ,{\mathrm{N}}_1)$

where ${\theta }_{\sigma }$ is the stress Lode angle, $I_1$ is the first invariant of the stress deviator, $J_2$ is the second invariant of the stress deviator, $c$ is the cohesion of the material, and $\phi$ is the internal friction angle of the material.

3.3. Flow Law Constraints

For ideal rigid plastic materials, the strain component obtained from the deformation compatibility condition should be equal to the plastic strain rate component obtained from the associated flow rule and yield condition. Since the strain within each element is continuous, the flow law constraint can be expressed as

(5)$\left\{\begin{array}{l}{\mathbf{B}}^e{\mathit{u}}^e-{\lambda }^e\nabla f^e({\mathit{\sigma }}^e)=0;(e=1,\dots ,{\mathrm{N}}_1)\\ {\lambda }^e{f}^e({\mathit{\sigma }}^e)=0;(e=1,\dots ,{\mathrm{N}}_1)\\ {\lambda }^e\ge 0;(e=1,\dots ,{\mathrm{N}}_1)\end{array}\right.$

where ${\mathrm{N}}_1$ is the number of elements, ${\mathbf{B}}^e$ is the compatibility matrix of the element, ${\mathit{u}}^e$ is the velocity vector of the element, ${\lambda }^e$ is the non-negative plastic multiplier of the element, $f^e(\mathit{\sigma })$ is the yield function, and the specific expressions for ${\mathbf{B}}^e$ and ${\mathit{u}}^e$ are as follows

(6)$\left\{\begin{array}{l}{\mathbf{B}}^e=\left[{\mathbf{B}}_1,{\mathbf{B}}_2,{\mathbf{B}}_3,{\mathbf{B}}_4\right]\\ {\mathbf{B}}_i={\displaystyle \frac{1}{6{\mathrm{V}}^e}}{\left[\begin{array}{cccccc}{b}_i& 0& 0& c_i& 0& d_i\\ 0& c_i& 0& b_i& d_i& 0\\ 0& 0& d_i& 0& c_i& b_i\end{array}\right]}^{\mathrm{T}}(i=1,2,3,4)\\ {\mathit{u}}^e=\left[{\mathit{u}}_1^e,{\mathit{u}}_2^e,{\mathit{u}}_3^e,{\mathit{u}}_4^e\right]\\ \nabla ={\left\{\partial /\partial {\sigma }_{xx},\partial /\partial {\sigma }_{yy},\partial /\partial {\sigma }_{zz},\partial /\partial {\gamma }_{xy},\partial /\partial {\gamma }_{yz},\partial /\partial {\gamma }_{zx}\right\}}^{\mathrm{T}}\end{array}\right.$

where ${\mathbf{B}}_i(i=1,2,3,4)$ is the compatibility matrix of the node $i$ of the element.

On the global kinematically admissible velocity field, the flow rule can be expressed as

(7)$\left\{\begin{array}{l}\mathbf{B}\mathit{u}-{\displaystyle \sum _{e=1}^E{\lambda }^e\nabla f^e(\mathit{\sigma })=0};(e=1,\dots ,{\mathrm{N}}_1)\\ {\lambda }^e\nabla f^e(\mathit{\sigma })=0;(e=1,\dots ,{\mathrm{N}}_1)\\ {\lambda }^e\ge 0;(e=1,\dots ,{\mathrm{N}}_1)\\ \end{array}\right.$

where $\mathbf{B}$ is the global compatibility, $\mathit{u}$ is the global continuous velocity vector, and $\mathit{\sigma }$ is the global stress vector.

In order to satisfy the flow law for a velocity discontinuity occurring on a common surface, it is imperative to apply constraints on the nodes associated with the discontinuity. The introduction of a node velocity variable in the local coordinate system between any adjacent nodes allows for the description of the velocity discontinuity, thereby guaranteeing that both tangential and normal velocities traverse the discontinuity to satisfy the flow law. Within the context of the Mohr–Coulomb yield criterion, the correlation between the jump in normal velocity and the jump in tangential velocity is established.

(8)$\Delta u_n=|\Delta u_t|\tan \phi$

where $\Delta u_n$ and $\Delta u_t$ are the normal and tangential velocity jump, respectively, $\phi$ is the internal friction angle of the material.

The tangential velocity jump vector $\Delta {{\mathit{u}}^{\prime }}_t^{lm}$ of each pair of nodes $(l,m)$ on the discontinuity is expressed by the sum of two non-negative velocity vectors ${\mathit{v}}^{+lm}$ and ${\mathit{v}}^{-lm}$.

(9)$\Delta {{\mathit{u}}^{\prime }}_t^{lm}={\mathit{v}}^{+lm}-{\mathit{v}}^{-lm}$

where ${\mathit{v}}^{+lm}={\left\{v_{x^{\prime }}^{+lm},v_{y^{\prime }}^{+lm}\right\}}^{\mathrm{T}}$, ${\mathit{v}}^{-lm}={\left\{v_{x^{\prime }}^{-lm},v_{y^{\prime }}^{-lm}\right\}}^{\mathrm{T}}$, $\Delta {{\mathit{u}}^{\prime }}_t^{lm}={\left\{\Delta {u^{\prime }}_{x^{\prime }}^{lm},\Delta {u^{\prime }}_{y^{\prime }}^{lm}\right\}}^{\mathrm{T}}$.

Then the overall velocity jump $\Delta {{\mathit{u}}^{\prime }}^{lm}$ in the local coordinate system can be expressed as a vector related to the tangential velocity jump $\Delta {{\mathit{u}}^{\prime }}_t^{lm}$ and the normal velocity jump $\Delta {u^{\prime }}_n^{lm}$

(10)$\Delta {{\mathit{u}}^{\prime }}^{lm}={\left\{\left\{\Delta {{\mathit{u}}^{\prime }}_t^{lm}\right\},\Delta {u^{\prime }}_n^{lm}\right\}}^{\mathrm{T}}$

In different coordinate systems, the tangential velocity jumps have the following relationship

(11)$\Delta {\mathit{u}}^{\prime }=\mathsf{\beta }\Delta \mathit{u}$

where $\mathsf{\beta }$ is the direction cosine matrix between the local coordinate system and the global coordinate system.

For the purpose of removing the absolute value symbol in Equation (8), and avoiding non-differentiable constraints, the tangential velocity jump can be expressed as the sum of two pairs of non-negative velocity variables $v_x^{\pm }$ and $v_y^{\pm }$.

(12)$\Delta u_t=(v_x^{+}+v_x^{-})+(v_y^{+}+v_y^{-})$

Substituting the above equation into Equation (8)

(13)$\Delta u_n=\left[(v_x^{+}+v_x^{-})+(v_y^{+}+v_y^{-})\right]\tan \phi$

We represent the flow law constraints on the velocity discontinuity within the kinematically admissible velocity field in the form of matrix equations

(14)${\mathbf{A}}_u^d{\mathit{u}}^d+{\mathbf{A}}_v^d{\mathit{v}}^d=\mathbf{0}$

where

(15)${\mathit{v}}^d\ge 0$

and

${\mathbf{A}}_u^d=\left[\left.\begin{array}{cccccc}{\mathsf{\beta }}^{14}& -{\mathsf{\beta }}^{14}& & & & \\ & & {\mathsf{\beta }}^{25}& -{\mathsf{\beta }}^{25}& & \\ & & & & {\mathsf{\beta }}^{36}& -{\mathsf{\beta }}^{36}\end{array}\right]\right.$

${\mathbf{A}}_u^d=\left[\left.\begin{array}{cccccc}{\mathsf{\beta }}^{14}& -{\mathsf{\beta }}^{14}& & & & \\ & & {\mathsf{\beta }}^{25}& -{\mathsf{\beta }}^{25}& & \\ & & & & {\mathsf{\beta }}^{36}& -{\mathsf{\beta }}^{36}\end{array}\right]\right.$

${\mathbf{A}}_u^d=\left[\left.\begin{array}{cccccc}{\mathsf{\beta }}^{14}& -{\mathsf{\beta }}^{14}& & & & \\ & & {\mathsf{\beta }}^{25}& -{\mathsf{\beta }}^{25}& & \\ & & & & {\mathsf{\beta }}^{36}& -{\mathsf{\beta }}^{36}\end{array}\right]\right.$

${\mathit{u}}^d=\left\{\left\{{\mathit{u}}^1\right\},\left\{{\mathit{u}}^4\right\},\left\{{\mathit{u}}^2\right\},\left\{{\mathit{u}}^5\right\},\left\{{\mathit{u}}^3\right\},\left\{{\mathit{u}}^6\right\}\right\}$

${\mathit{v}}^d=\left\{{\left\{{\mathit{v}}^{+14}\right\}}^{\mathrm{T}},{\left\{{\mathit{v}}^{-14}\right\}}^{\mathrm{T}},{\left\{{\mathit{v}}^{+36}\right\}}^{\mathrm{T}},{\left\{{\mathit{v}}^{-36}\right\}}^{\mathrm{T}},{\left\{{\mathit{v}}^{+25}\right\}}^{\mathrm{T}},{\left\{{\mathit{v}}^{-25}\right\}}^{\mathrm{T}}\right\}$

where $d=(1,\cdots ,{\mathrm{N}}_2)$, ${\mathrm{N}}_2$ is the number of common faces, ${\mathsf{\beta }}^{l_i{m}_i}$ is the coordinate change matrix between the local coordinate system and the global coordinate system of the pair of nodes $(l_i,m_i),i(i=1,2,3)$ on the discontinuity, ${{\mathbf{A}}_u}^d$ is the matrix of equality constraint coefficients for the continuous velocity, ${{\mathbf{A}}_v}^d$ is the matrix of equality constraint coefficients for the discontinuous velocity, ${\mathit{u}}^d$ is the continuous velocity vector of common faces, and ${\mathit{v}}^d$ is the discontinuous velocity vector of common faces.

3.4. Velocity Boundary Condition

Assuming that the prescribed velocity distribution in the global coordinate system is ${\mathbf{B}}^b$, the velocity boundary constraint equation is

(16)${\mathbf{A}}^b{\mathit{u}}^b={\mathbf{B}}^b$

where ${\mathbf{B}}^b={\left[{\mathit{u}}_x^b,{\mathit{u}}_y^b,{\mathit{u}}_z^b\right]}^{\mathrm{T}}$, ${\mathit{u}}_x^b,{\mathit{u}}_y^b,{\mathit{u}}_z^b$ are the components of prescribed velocity vectors in different directions in a global coordinate system, ${\mathbf{A}}^{\mathrm{b}}$ is the coordinate transformation matrix, and ${\mathit{u}}^b={\left[{\mathit{u}}_{x^{\prime }},{\mathit{u}}_{y^{\prime }},{\mathit{u}}_{z^{\prime }}\right]}^{\mathrm{T}}$ is the components of known velocity vectors in a local coordinate system. $b=(1,\cdots ,{\mathrm{N}}_3)$, ${\mathrm{N}}_3$ is the number of boundary triangle elements with known velocities.

3.5. Constraints on Equality of Internal and External Power

Owing to the potential for plastic deformation to take place either within the components or at the velocity discontinuity, the internal dissipated power is comprised of two distinct parts. Concerning the dissipated power within the components, it may be articulated as follows

(17)$W_c^e={\displaystyle {\int }_{V^e}{\mathit{\sigma }}^{\mathrm{T}}\dot{\mathit{\epsilon }}}dV={\displaystyle {\int }_{V^e}{\mathit{\sigma }}^{\mathrm{T}}{\mathbf{B}}^e{\mathit{u}}^e}dV=V^e{\mathit{\sigma }}^{eT}{\mathbf{B}}^e{\mathit{u}}^e$

where $W_c^e$ is the dissipated power within the element, $V^e$ is the volume of the element, ${\mathit{\sigma }}^e$ is the stress vector of the element, and ${\mathbf{B}}^e{\mathit{u}}^e$ denotes the strain rate of the element and is given by Equation (6).

The dissipated power on the velocity discontinuity with an area of $S_d$ can be calculated using the following equation

(18)$W_d^s={\displaystyle {\int }_{S_d}c}|\Delta u_t|dS={\displaystyle {\int }_{S_d}c}\left[(v_x^{+}+v_x^{-})+(v_y^{+}+v_y^{-})\right]dS$

where $W_d^s$ is the dissipated power on the velocity discontinuity, $S_d$ is the area of the velocity discontinuity, $c$ is the cohesion of the material, and the tangential velocity jump $|\Delta u_t|$ is given by Equation (12). The external power applied by the self-weight on the element is expressed as

(19)$W_g^e={{\displaystyle {\int }_{V^e}\mathit{g}}}^{\mathrm{T}}\mathit{u}dV={{\displaystyle {\int }_{V^e}{\mathit{g}}^e}}^{\mathrm{T}}{\mathbf{N}}^e{\mathit{u}}^{e}dV={{\mathit{g}}^e}^{\mathrm{T}}({\displaystyle {\int }_{V_e}{\mathbf{N}}^e}dV){\mathit{u}}^e$

where ${\mathbf{N}}^e=\left[\mathbf{E}{\mathrm{N}}_1,\mathbf{E}{\mathrm{N}}_2,\mathbf{E}{\mathrm{N}}_3,\mathbf{E}{\mathrm{N}}_4\right]$, ${\mathbf{g}}^e={\left[g_1,g_2,g_3,g_3\right]}^{\mathrm{T}}$, $\mathbf{E}$ is an identity matrix with 4 × 4, and ${\mathbf{g}}^e$ is the column vector of the load act on the node by the self-weight of the element.

The external power exerted by the external load act on the element is

(20)$W_q^e={{\displaystyle {\int }_{A_q}\mathit{q}}}^{\mathrm{T}}{\mathit{u}}^{e}dA$

where $\mathrm{A}q$ is the area of load action and the load component applied by the external load on the node is represented as ${\mathit{q}}^{\mathrm{T}}={\left\{{\left\{{\mathit{q}}_1\right\}}^{\mathrm{T}},{\left\{{\mathit{q}}_2\right\}}^{\mathrm{T}},{\left\{{\mathit{q}}_3\right\}}^{\mathrm{T}},{\left\{{\mathit{q}}_4\right\}}^{\mathrm{T}}\right\}}^{\mathrm{T}}$.

The UBM indicates that the internal power dissipation and external power dissipation of the element are equal when it is destroyed. Specifically

(21)$W_c^e+W_s^d=W_g^e+W_q^e$

In order to obtain the upper-bound value of the external load, Equation (20) can be expressed as

(22)$Q={\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}\mathit{u}+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}$

where $Q$ is the external load, ${\mathit{\sigma }}^T$ is the transposition of stress vector, $\mathit{u}$ and $\mathit{v}$ are the continuous and discontinuous velocity vectors, respectively, ${\mathit{c}}_u$, ${\mathit{c}}_v$ are the coefficient vectors corresponding to the continuous and discontinuous velocity vectors, and $\mathbf{B}$ is the global compatibility matrix.

3.6. Objective Function

The primary purpose of this paper is to obtain the stability factor of the karst foundation and the roof’s horizontal stress. Since the UBM analyzes the failure mechanism of geotechnical structures, establishing the kinematically admissible velocity fields, there is no need to construct a constitutive model of geotechnical structures. If solved directly, the stress will be eliminated during the calculation process, only the structure’s upper-bound solution and corresponding velocity field can be obtained. Nevertheless, by transforming the original linear programming into a set of Kuhn–Tucker conditions, while obtaining the objective function value, the velocity and stress fields of the geotechnical structure can also be obtained. Therefore, when establishing a nonlinear mathematical model of the karst cave roof’s ultimate bearing capacity, the external load, is set as the objective function, and its minimum value is calculated. The stability factor of the karst foundation is obtained by reducing the shear strength parameters of rock while solving the nonlinear programming model. According to Equation (21), the objective function can be written as:

(23)$\begin{array}{cc}Minimize& Q=\end{array}{\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}u+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}$

3.7. Nonlinear Optimization Mathematical Model

In the previous section, the basic formula of the UBM was introduced, including the discretization of elements, flow rule constraints, velocity boundary conditions, constraint conditions for equally internal and external power, and the setting of the objective function. Then we assemble the Equations (4), (5), (14)–(16) and (23). Finally, a mathematical model is obtained, which can be used to calculate the ultimate bearing capacity of the karst foundation:

(24)$\begin{array}{ll}Minimize& Q={\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}\mathit{u}+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}\\ Subject to& \left\{\begin{array}{l}\mathbf{B}\mathit{u}-{\displaystyle \sum _{e=1}^E{\lambda }^e\nabla f^e(\mathit{\sigma })=0};e=(1,\cdots ,{\mathrm{N}}_1)\\ {\mathbf{A}}_u^d{\mathit{u}}^d+{\mathbf{A}}_v^d{\mathit{v}}^d=\mathbf{0}\mathbf{;}d=(1,\cdots ,{\mathrm{N}}_2)\\ {\mathbf{A}}^b{\mathit{u}}^b={\mathbf{B}}^b;b=(1,\cdots ,{\mathrm{N}}_3)\\ {\lambda }^e{f}^e(\mathit{\sigma })=0;e=(1,\cdots ,{\mathrm{N}}_1)\\ f^e(\mathit{\sigma })\le 0;e=(1,\cdots ,{\mathrm{N}}_1)\\ {\lambda }^e\ge 0;e=(1,\cdots ,{\mathrm{N}}_1)\\ v^d\ge 0;d=(1,\cdots ,{\mathrm{N}}_2)\\ \end{array}\right.\end{array}$

4. Solution

The feasible region of the original linear programming problem is a convex polyhedron, therefore, the solution, satisfying the Kuhn–Tucker condition, must be the optimal solution of the original linear programming problem. Note the constraint equation of Equation (24)

(25)$\left\{\begin{array}{l}\mathbf{B}\mathit{u}-{\displaystyle \sum _{e=1}^E{\lambda }^e\nabla f^e(\mathit{\sigma })=0};e=(1,\cdots ,{\mathrm{N}}_1)\\ {\lambda }^e{f}^e(\mathit{\sigma })=0;\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} e=(1,\cdots ,{\mathrm{N}}_1)\\ f^e(\mathit{\sigma })\le 0;\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} e=(1,\cdots ,{\mathrm{N}}_1)\\ {\lambda }^e\ge 0;\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} e=(1,\cdots ,{\mathrm{N}}_1)\end{array}\right.$

is the Kuhn–Tucker optimality condition for the nonlinear programming problem (26)

(26)$\begin{array}{l}Maximize Q={\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}\mathit{u}+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}(on stress field)\\ Subject to f^e(\mathit{\sigma })\le 0;e=(1,\cdots ,{\mathrm{N}}_1)\end{array}$

Simultaneous Equations (24) and (25) obtain a nonlinear programming problem of the double objective function, which can simultaneously solve the stress field and velocity field of the karst foundation in the critical state.

(27)$\begin{array}{l}Maximize Q={\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}\mathit{u}+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}(on stress field)\\ Minimize Q={\mathit{\sigma }}^{\mathrm{T}}\mathbf{B}\mathit{u}+{\mathbf{c}}_u^{\mathrm{T}}\mathit{u}+{\mathbf{c}}_v^{\mathrm{T}}\mathit{v}(on velocity field)\\ Subject to\left\{\begin{array}{l}{\mathbf{A}}_u^d{\mathit{u}}^d+{\mathbf{A}}_v^d{\mathit{v}}^d=\mathbf{0}\mathbf{;}d=(1,\cdots ,{\mathrm{N}}_2)\\ {\mathbf{A}}^b{\mathit{u}}^b={\mathbf{B}}^b;\hspace{1em}\hspace{1em} b=(1,\cdots ,{\mathrm{N}}_3)\\ f^e(\mathit{\sigma })\le 0;\hspace{1em}\hspace{1em} e=(1,\cdots ,{\mathrm{N}}_1)\\ f^e(\mathit{v})\le 0;\hspace{1em}\hspace{1em} e=(1,\cdots ,{\mathrm{N}}_1)\end{array}\right.\end{array}$

Writing the constraint ${\mathit{v}}^d\ge \mathbf{0}$ as $f({\mathit{v}}^d)\le 0$, the Kuhn–Tucker optimality condition of Equation (27) can be expressed as

(28)$\left\{\begin{array}{c}\mathbf{B}{\mathbf{u}}^{*}-\sum _{e=1}^E{\lambda }_{\sigma }^E\nabla f^{\sigma }\left({\mathit{\sigma }}^{*}\right)=0;e=(1,\cdots ,{\mathrm{N}}_1)\hfill \\ {\mathbf{c}}_v+{\left({\mathbf{A}}_v^d\right)}^{\mathrm{T}}\mathsf{\mu }+\sum _{e=1}^E{\lambda }_v^e{f}^e\left({\mathit{v}}^{*}\right)=0;e=(1,\cdots ,{\mathrm{N}}_1)\hfill \\ {\mathbf{B}}^{\mathrm{T}}{\mathit{\sigma }}^{*}+{\mathbf{c}}_u+{\mathbf{A}}_u^{\mathrm{T}}\mathsf{\mu }=0\hfill \\ {\mathbf{A}}_u^d{\mathit{u}}^d+{\mathbf{A}}_v^d{\mathit{v}}^d=\mathbf{0}\mathbf{;}d=(1,\cdots ,{\mathrm{N}}_2)\hfill \\ {\mathbf{A}}^b{\mathit{u}}^b={\mathbf{B}}^b;b=(1,\cdots ,{\mathrm{N}}_3)\hfill \\ {{\lambda }_{\sigma }}^e\nabla f^e\left({\mathit{\sigma }}^{*}\right)=0;e=(1,\cdots ,{\mathrm{N}}_1)\hfill \\ {{\lambda }_v}^{e}f\left({\mathit{v}}^{*}\right)=0;d=(1,\cdots ,{\mathrm{N}}_2)\hfill \\ f^e\left({\mathit{\sigma }}^{*}\right)\le 0;e=(1,\cdots ,{\mathrm{N}}_1)\hfill \\ f^e\left({\mathit{v}}^{*}\right)\le 0;e=(1,\cdots ,{\mathrm{N}}_1)\hfill \\ {\lambda }_{\sigma }^e\ge 0,{\lambda }_v^e\ge 0\hfill \end{array}\right.$

where ${\lambda }_{\sigma }^e$, ${\lambda }_v^e$ are the non-negative plastic multipliers on the velocity and stress fields of each element, respectively, $({\mathbf{u}}^{*},{\mathbf{v}}^{*},{\mathsf{\sigma }}^{*})$ is the optimal solution of the original linear programming problem, and $\mathsf{\mu }$ is a coefficient vector for the discontinuity variables.

Following the discussion above, the initial nonlinear programming problem is converted into a series of Kuhn–Tucker optimization conditions, which are subsequently resolved via the two-stage quasi-Newton algorithm. The fundamental premise of the quasi-Newton and Newton methods is to linearize the nonlinear equation by employing the Taylor expansion at the iteration point. During each iteration process, the iteration direction is aligned with the decrease in function value at the current point. The iteration continues until the minimum value of the objective function, which fulfills the accuracy requirement, is obtained. However, it should be noted that the Newton method necessitates the calculation of the Hessian matrix or inverse matrix during each iteration, which is time-consuming and may also be challenging or infeasible to compute. In contrast, the quasi-Newton method employs a positive definite symmetric matrix to approximate the Hessian matrix or its inverse matrix, facilitating a more manageable resolution of the problem.

The UBM program, designed for the analysis of the roof’s bearing capacity, was ultimately compiled in Python. The detailed analysis procedure is illustrated in Figure 2. The computational hardware utilized in calculations comprises AMD Ryzen Threadripper 2990 WX (AMD, Santa Clara, CA, USA), with a primary frequency of 3.0 GHz, 32 cores, and 128 GB memory.

5. Validation and Application

5.1. Numerical Example

A typical karst foundation example is selected for analysis to verify the effectiveness of the method proposed in this paper. The schematic diagram of the foundation model is shown in Figure 3. The length L of the model is 30 m, the width B is 30 m, and the height H is (T + h + 3) m. When other conditions are unchanged, the height of the karst cave has little influence on the karst foundation’s ultimate bearing capacity [25,26], the height h of the karst cave is taken to be equal to its width W, the type of the rock is limestone with a standard tensile strength value of 0.6 MPa [27]. Above the rock layer, with a 10 m clay layer, the self-weight of the clay is directly applied to the rock layer in the form of a uniform load. The pile foundation load is 400 MPa. The calculation of the mechanical parameters of rock and soil mass is shown in Table 1.

Unlike the previous research by Olmeda, J.L.B. [2], this paper considers the influences of different pile diameters D, dimensionless parameters T/D (ratio of the roof’s thickness T to pile diameter D), and W/D (ratio of the karst cave’s width W to pile diameter D) on the karst foundation’s ultimate bearing capacity. The calculated pile diameters D are 0.8 m, 1.0 m, 1.5 m, and 2.0 m, respectively; The roof’s thickness T varies from 1 D to 7 D, which takes integer multiples of the pile diameter. The karst cave’s width W varies from 1 m to 8 m, which takes an integer. A total of 224 different parameters are used to analyze the karst foundation’s ultimate bearing capacity. Furthermore, the mesh shown in Figure 4, the adaptive mesh, is used for refinement, which is divided into three steps. The velocity boundary conditions are:

• (1). The velocities around the model and at the bottom are 0;

• (2). The velocity in the z direction of the upper boundary is −1;

• (3). The velocities in the x and y directions are both 0.

5.2. Horizontal Stress Analysis

The capacity of a rock mass to withstand tensile deformation is commonly referred to as its tensile strength. As a result, investigating the tensile behavior of rock masses is of significant importance in the context of surrounding rock in underground engineering caverns [21]. The horizontal stress of the roof, as determined by the analysis presented above, is derived from the stress field within the calculation results. Strength analysis of the roof is then conducted with varying pile diameters of 0.8 m and 1.5 m, as illustrated in Figure 5 and Figure 6, respectively, depicting the influence curves of different T/D and W/D on the roof’s horizontal stress. It is evident that the junction of the pile and rock mass is subjected to extreme tensile stress, while the surrounding area experiences compressive stress. This local state of tension and compression is unfavorable for the roof’s stability; however, the reduction of horizontal stress can be attained through the alteration of the T/D and W/D parameters. As illustrated in Figure 5, a decrease in the roof’s tensile and compressive stress can be achieved by increasing T/D when T/D is less than 3. For instance, in Figure 5, an increase in T from 1 D to 2 D results in over 70% reduction in the extreme tensile stress values of the roof, decreasing by 2 MPa and 0.6 MPa, respectively. Similarly, the extreme compressive stress values decrease by 1 MPa and 0.2 MPa, respectively, both decreasing by more than 60%. However, the influence curves of different T/D values overlap, indicating a feeble impact on the roof’s horizontal stress when T/D is greater than 3. Nonetheless, the extreme values of the roof’s horizontal tensile stress are less than the tensile strength of the rock mass by 0.6 MPa, implying that the horizontal stress does not affect the overall stability of the roof. Figure 6 demonstrates that the reduction of the karst cave’s horizontal stress can also be achieved by decreasing W/D, particularly for W spans greater than 6 m. The alteration of W/D has the most substantial influence on the roof’s horizontal stress value. Of significance is the fact that even when T/D is equal to 3, and the karst cave’s width, W is equal to 8 m, the extreme value of the roof’s horizontal stress still exceeds the tensile strength of the rock mass. Henceforth, meticulous evaluation of the W/D and T/D parameters is imperative in ensuring the durability of the karst cave ceiling.

We select the maximum horizontal stress value of the roof, and subsequently produce stress contour cloud maps under various T/D and W/D parameters, as depicted in Figure 7. In the stress contour map, the upper left region which lies under the −0.6 MPa stress contour line, denotes the extreme tensile stress values of the roof. It should be noted that said values are all lesser than the roof’s corresponding tensile strength. Thus, the selection of any T/D and W/D combination falling within this specific region would ensure that the roof’s stability remains unaffected by horizontal stress. Notably, as the pile diameter D increases, the aforementioned −0.6 MPa stress contour line gradually shifts toward the lower right of the diagram. This constitutes that given equal W/D values, a larger pile diameter D results in a larger value of T/D required to prevent the extreme horizontal tensile stress from surpassing the tensile strength value of the roof. For instance, if W/D values are both 6, the minimum T/D value necessary to prevent the roof’s extreme tensile stress from surpassing its corresponding tensile strength value is 2.5 when the pile diameter is 0.8 m and 3.5 when the pile diameter is 1.5 m, as shown in Figure 7a,b, respectively.

5.3. Stability Analysis

Based on the previous discourse, the stability factor k is derived through the solution of a nonlinear mathematical programming model that determines the karst cave roof’s ultimate bearing capacity under various scenarios. Table 2, Table 3, Table 4 and Table 5 illustrate that a change in T/D significantly affects the stability factor k when T/D is lesser than 3, although all other conditions remain constant. As T/D increases, the stability factor k shows a pronounced upward trend. When T/D reaches 3, any alterations in the stability factor k due to T/D become negligible, indicating that T/D barely impacts the roof’s ultimate bearing capacity. Furthermore, an increase in W/D weakens the karst cave roof’s ultimate bearing capacity, as evidenced by the gradual decline in the stability factor k with an increase in W/D. It is noteworthy that when the T/D is 1, and the karst cave’s span W exceeds 6 m, the roof’s stability factor k for each of the four varying pile diameters is less than 1, suggesting a propensity for damage. Consequently, during actual engineering, special attention must be paid to karst caves with spans exceeding 6 m, consistent with the conclusion reached through numerical simulation in reference [28].

In order to better illustrate the impact of two dimensionless parameters (T/D and W/D) on the stability factor k, a graphical representation has been created based on the calculation results provided in Table 4, as depicted in Figure 8. The findings derived from Figure 8 indicate that the stability factor k of the roof can be significantly influenced by low (W/D < 2) or high (W/D > 5) values of W/D. When W/D is less than 2 and T/D is less than 3, the roof’s stability factor k is relatively high due to the load-sharing attribute of the rock mass surrounding the karst cave with the karst cave roof, which in turn enhances the overall roof stability. On the other hand, when W/D exceeds 5 and T/D exceeds 3, the stability factor k of the roof is also higher due to the limited stress diffusion range within the karst cave roof when the span is too large. However, the high T/D value reduces the potential for pile foundation load to cause damage to the roof, resulting in a skewed calculation outcome that overemphasizes safety in the evaluation of karst foundation stability.

Figure 9 displays the contour cloud map that indicates the stability factor k of the roof under varying T/D and W/D. It is imperative that the stability of the karst cave roof is designed with an appropriate safety reserve. Figure 9 can be referred to in order to obtain the T/D and W/D values that comply with the design requirements. For example, if the design necessitates a safety factor k of 2.0 for the karst cave roof and the width of the cave is 5 m, the values of T/D and W/D can be selected from the upper left area of the “safety factor 2.0 contour line” in Figure 9.

5.4. Discussion

Through the analysis of the horizontal stress present in the roof, it has been determined that the roof experiences both tensile and compressive stress when supporting the load of the pile foundation. In order to minimize the destabilizing effects of this stress on the roof, increasing the T/D has been found to be effective. However, it is important to consider the W/D to ensure the roof’s stability.

The stability factor of the roof displays a significant increase during the initial stages (T/D < 3), ultimately approaching a fixed value during the later stages (T/D > 3). This pattern can be attributed to the different failure modes exhibited by the roof when T/D is less or greater than 3 [29]. Wu, B. et al. [24] conducted a study utilizing the strength reduction method to examine the safety factor of underground caverns. Their research revealed that as the rock mass deformation in the vicinity of the cavern reaches a certain threshold, the reduction factor gradually increases until it reaches a steady state. These findings bear similarity to those presented in this article.

It is important to highlight that in this paper, the stability factor of the roof has been considerably increased in both cases (T/D > 3 and W/D > 5; T/D < 2 and W/D < 1) through the calculation of 224 distinct working scenarios. The stress diffusion range of the roof precludes load-induced cave damage, thus increasing the cave’s stability.

The paper utilizes a shear strength reduction parameter (stability factor) to convey the bearing capacity of karst foundations, a method that is not widely employed in the stability analysis of karst foundations. To ensure the accuracy of this method, the method is compared to the three-dimensional FEM and the two-dimensional limit analysis [22]. The safety factor calculation outcomes of the three methods are presented in Table 6 (with D = 0.8 m, W = 3.2 m, and the same physical and mechanical parameters as in Table 1). The calculation results of this method are in alignment with the other two methods for T/D values of less than 3, signifying the accuracy of this method. However, for T/D values greater than 3, the calculation errors of the three methods begin to increase, potentially resulting from four reasons: (1) the similarity calculation of the conservative model and display function; (2) the linearization of nonlinear constraints in nonlinear mathematical programming models using the UBM; (3) the division of the mesh; (4) the influence of three-dimensional spatial effect. The two-dimensional and three-dimensional limit analyses have the largest calculation errors that continue to remain within 10% (when T/D is equal to 4), validating the effectiveness of the proposed method in this paper for the stability analysis of three-dimensional karst foundations.

6. Conclusions

This paper examines the stability of a three-dimensional karst foundation by employing the UBM of finite limit analysis. The primary conclusions are outlined below:

• (1). The horizontal stress experienced by the roof due to the pile foundation load can be classified into tensile and compressive stresses. In order to mitigate the detrimental effects of these stresses on the bearing capacity of the roof, increasing the T/D proves to be a more effective approach than decreasing the W/D. This is especially true when the T/D is less than 3, as a higher T/D is associated with a more remarkable reduction in the horizontal stress value of the roof.

• (2). When the W/D is equal, an increase in the pile diameter will necessitate a corresponding higher T/D to prevent the surpassing of horizontal stress over the tensile strength of the rock mass.

• (3). The stability factor of the roof demonstrates an initial linear increase with the elevation of the T/D value (T/D < 3), followed by a gradual stabilization (T/D > 3). This suggests that when T/D exceeds 3, the impact of roof thickness on bearing capacity becomes negligible. Furthermore, the stability factor of the roof only experiences further enhancement when W/D exceeds 5.

• (4). The stability factor of the three-dimensional karst foundation is more cautious compared to the outcomes obtained under two-dimensional conditions. Upon reaching T/D greater than 3, the calculation errors under two-dimensional and three-dimensional circumstances start to increase. Therefore, considering the three-dimensional spatial effects, examining the bearing capacity of the karst foundation aligns better with engineering practice.

• (5). This paper solely looks into the stability of regular rectangular karst caves’ end-bearing pile karst foundation. However, real-life engineering projects often encounter more complex geological conditions such as karst caves with different shapes and multi-layer karst caves. Furthermore, the bearing capacity of rock socketed pile karst foundations has rarely been considered, and hence, these issues merit further investigation. Adding experimental results to scrutinize the research content will enhance its reliability and substance.

Author Contributions

Conceptualization, Z.L., K.L., W.L., H.W., P.P. and H.X.; methodology, K.L.; software, Z.L. and K.L.; validation, W.L., H.W. and K.L.; formal analysis, Z.L.; investigation, Z.L.; resources, W.L.; data curation, W.L.; writing—original draft preparation, K.L.; writing—review and editing, Z.L. and P.P.; visualization, H.X.; supervision, W.L. and Z.L.; project administration, H.X. and Z.L.; funding acquisition, Z.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

• The following symbols are used in this paper:

$c$	=	cohesion before strength reduction
$c^{\prime }$	=	cohesion after strength reduction
$\phi$	=	internal friction angle before strength reduction
${\phi }^{\prime }$	=	internal friction angle after strength reduction
k	=	stability factor
${\mathit{u}}^i$	=	velocity of the node i
$a_i,b_i,c_i,d_i$	=	coefficients related to the geometric position of the node i.
${\theta }_{\sigma }$	=	stress Lode angle
$I_1$	=	the first invariant of the stress deviator
$J_2$	=	the second invariant of the stress deviator
${\mathbf{B}}^e$	=	compatibility matrix of the element
$\mathbf{B}$	=	global compatibility
${\mathit{u}}^e$	=	continuous velocity vector of the element
$\mathit{u}$	=	continuous velocity vector
$\mathit{v}$	=	discontinuous velocity vector
${\mathit{\sigma }}^e$	=	stress vector of the element
$\mathit{\sigma }$	=	global stress vector
${\lambda }^e$	=	non-negative plastic multiplier of the element
${\lambda }_{\sigma }^e$	=	non-negative plastic multipliers on the stress filed of the element
${\lambda }_v^e$	=	non-negative plastic multipliers on the velocity filed of the element
$\Delta {{\mathit{u}}^{\prime }}^{lm}$	=	velocity jump vector of each pair of nodes $(l,m)$ on the discontinuity
$\Delta {{\mathit{u}}^{\prime }}_t^{lm}$	=	tangential velocity jump vector of each pair of nodes $(l,m)$ on the discontinuity
${\mathit{v}}^{+lm}$, ${\mathit{v}}^{-lm}$	=	non-negative velocity vectors of tangential velocity jump
$\Delta u_t$	=	tangential velocity jump
$\Delta u_n$	=	normal velocity jump
${{\mathbf{A}}_u}^d$	=	matrix of equality constraint coefficients for the continuous velocity
${{\mathbf{A}}_v}^d$	=	matrix of equality constraint coefficients for the discontinuous velocity
${\mathit{u}}^d$	=	continuous velocity vector of common faces
${\mathit{v}}^d$	=	discontinuous velocity vector of common faces
$\mathsf{\beta }$	=	direction cosine matrix of velocity jump between local coordinate system and global coordinate system
${\mathbf{B}}^b$	=	prescribed velocity distribution in the global coordinate system
${\mathbf{A}}^b$	=	coordinate transformation matrix between the global and local coordinate system
${\mathit{u}}^b$	=	components of known velocity vectors in a local coordinate system
$W_c^e$	=	dissipated power within the element
$V^e$	=	volume of the element
${\mathit{\sigma }}^e$	=	stress vector of the element
$W_d^s$	=	dissipated power on the velocity discontinuity
$S_d$	=	area of the velocity discontinuity
${\mathbf{g}}^e$	=	column vector of the load act on the node by the self-weight
$\mathrm{A}q$	=	area of load action
$Q$	=	external load
$\mathsf{\mu }$	=	coefficient vector for the discontinuity variables

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Upper-bound tetrahedral element: (a) Velocity and Stress of tetrahedral element, (b) Velocity discontinuity.

Figure 2. Program block diagram of nonlinear mathematical programming model for UBM of karst foundation’s ultimate bearing capacity.

Figure 3. Size of karst foundation calculation model: (a) Main View, (b) Vertical view.

Figure 4. Karst foundation model mesh: (a) Overall mesh, (b) Profile mesh.

Figure 5. Influence of different T/D on the roof’s horizontal stress: (a) D = 0.8 m, (b) D = 1.5 m.

Figure 6. Influence of different W/D on the roof’s horizontal stress: (a) D = 0.8 m, (b) D = 1.5 m.

Figure 7. Cloud chart of roof’s horizontal stress contour under different T/D and W/D: (a) D = 0.8 m, (b) D = 1.5 m.

Figure 8. Influence curve between k and two dimensionless parameters; (a) T/D and k; (b) W/D and k.

Figure 9. Contour cloud map of the roof’s stability factor k under different T/D and W/D: (a) D = 0.8 m, (b) D = 1.0 m, (c) D = 1.5 m, (d) D = 2.0 m.

Figure 9. Contour cloud map of the roof’s stability factor k under different T/D and W/D: (a) D = 0.8 m, (b) D = 1.0 m, (c) D = 1.5 m, (d) D = 2.0 m.

Table 1

Calculation mechanical parameters of rock and soil mass.

Material Type	$\mathit{c}$/(kPa)	$\mathit{\phi }$/(°)	$\mathit{\gamma }$/(kN/m3)
clay	10	20	20
Concrete pile	50,000	42	27
limestone	700	36	25

Table 2

Stability factor k of karst cave roof when pile diameter D = 0.8 m.

	k
	W = 1 m	W = 2 m	W = 3 m	W = 4 m	W = 5 m	W = 6 m	W = 7 m	W = 8 m
T = 0.8 m (1D)	1.94	1.23	1.00	1.00	1.00	0.97	0.93	0.86
T = 1.6 m (2D)	2.73	2.26	1.95	1.95	1.95	1.95	1.97	1.94
T = 2.4 m (3D)	2.76	2.75	2.73	2.70	2.73	2.78	2.73	2.75
T = 3.2 m (4D)	2.73	2.75	2.75	2.74	2.78	2.78	2.80	2.86
T = 4.0 m (5D)	2.73	2.74	2.75	2.75	2.73	2.80	2.83	2.80
T = 4.8 m (6D)	2.74	2.75	2.74	2.74	2.77	2.77	2.79	2.98
T = 5.6 m (7D)	2.72	2.76	2.75	2.73	2.75	2.80	2.79	2.85

Table 3

Stability factor k of karst cave roof when pile diameter D = 1.0 m.

	k
	W = 1 m	W = 2 m	W = 3 m	W = 4 m	W = 5 m	W = 6 m	W = 7 m	W = 8 m
T = 1 m (1D)	1.94	1.23	1.00	1.00	1.00	0.99	0.92	0.87
T = 2 m (2D)	2.76	2.23	1.94	1.95	1.95	1.97	1.97	1.92
T = 3 m (3D)	2.76	2.75	2.73	2.71	2.71	2.73	2.70	2.77
T = 4 m (4D)	2.63	2.74	2.75	2.79	2.80	2.81	2.98	2.84
T = 5 m (5D)	2.75	2.75	2.75	2.78	2.79	2.83	2.92	3.14
T = 6 m (6D)	2.75	2.73	2.77	2.76	2.79	2.87	3.02	2.95
T = 7 m (7D)	2.75	2.76	2.77	2.78	2.83	2.85	2.94	3.36

Table 4

Stability factor k of karst cave roof when pile diameter D = 1.5 m.

	k
	W = 1 m	W = 2 m	W = 3 m	W = 4 m	W = 5 m	W = 6 m	W = 7 m	W = 8 m
T = 1.5 m (1D)	1.94	1.23	1.01	1.00	0.99	0.96	0.92	0.85
T = 3.0 m (2D)	2.73	2.22	1.91	1.90	1.95	1.94	1.94	1.86
T = 4.5 m (3D)	2.75	2.75	2.74	2.66	2.66	2.66	2.71	2.69
T = 6.0 m (4D)	2.75	2.76	2.75	2.79	2.87	2.89	2.87	3.04
T = 7.5 m (5D)	2.75	2.78	2.75	2.76	2.79	2.84	2.90	3.08
T = 9.0 m (6D)	2.67	2.75	2.74	2.76	2.81	2.85	2.82	3.07
T = 10.5 m (7D)	2.72	2.76	2.75	2.73	2.75	2.80	2.79	2.85

Table 5

Stability factor k of karst cave roof when pile diameter D = 2.0 m.

	k
	W = 1 m	W = 2 m	W = 3 m	W = 4 m	W = 5 m	W = 6 m	W = 7 m	W = 8 m
T = 2.0 m (1D)	1.94	1.22	0.99	0.98	0.98	0.94	0.90	0.83
T = 4.0 m (2D)	2.74	2.21	1.90	1.90	1.89	1.91	1.86	1.75
T = 6.0 m (3D)	2.74	2.76	2.67	2.62	2.60	2.64	2.65	2.56
T = 8.0 m (4D)	2.74	2.76	2.74	2.77	2.80	2.85	2.86	2.85
T = 10.0 m (5D)	2.77	2.76	2.74	2.77	2.91	2.87	2.79	3.06
T = 12.0 m (6D)	2.74	2.74	2.74	2.74	2.77	2.84	2.98	2.91
T = 14.0 m (7D)	2.73	2.74	2.76	2.74	2.76	2.87	2.91	3.51

Table 6

Stability factor k of karst cave roof using different methods.

T/D (D = 0.8 m)	FEM (Three-Dimensional)	Reference [22] (2021, Two-Dimensional)	This Paper (Three-Dimensional)
1	1.02	1.03	1.00
2	1.95	1.96	1.95
3	2.75	2.76	2.73
4	2.79	2.84	2.75
5	2.83	2.94	2.75

References

1. Wang, P.S.; Ding, H.Y.; Zhang, P.Y. Influence of Karst Caves at Pile Side on the Bearing Capacity of Super-Long Pile Foundation. Math. Probl. Eng.; 2020; 2020, 4895735. [DOI: https://dx.doi.org/10.1155/2020/4895735]

2. Olmeda, J.L.B.; Robles, J.M.; Perez, E.S.; Maranon, C.O. Influence of Natural Cavities on the Design of Shallow Foundations. Appl. Sci.; 2020; 10, 1119. [DOI: https://dx.doi.org/10.3390/app10031119]

3. Zhang, Q.Q.; Qiao, S.S.; Xing, Y.C.; Zhang, K.; Cui, W.; Wang, Z.Y. Study on Bearing Characteristics of Single Pile Crossing Unfilled Karst Cave. J. Hunan Univ. Nat. Sci.; 2022; 49, pp. 186-196. [DOI: https://dx.doi.org/10.16339/j.cnki.hdxbzkb.2022087]

4. Xie, P.; Duan, H.C.; Wen, H.J.; Yang, C.; Ma, S.K.; Yue, Z.R. Study on a Quantitative Indicator for Surface Stability Evaluation of Limestone Strata with a Shallowly Buried Spherical Karst Cave. Mathematics; 2022; 10, 2149. [DOI: https://dx.doi.org/10.3390/math10122149]

5. Xie, P.; Duan, H.C.; Wen, H.J.; Yang, C.; Ma, S.K.; Yue, Z.R. The bearing capacity analysis of limestone strata roof containing a shallow-buried cylinder Karst cave. Mech. Adv. Mater. Struct.; 2022; 29, pp. 3421-3428. [DOI: https://dx.doi.org/10.1080/15376494.2021.1900466]

6. Lai, V.; Shiau, J.; Van, C.N.; Tran, H.D. Bearing capacity of conical footing on anisotropic and heterogeneous clays using FEA and ANN. Mar. Georesources Geotechnol.; 2022; pp. 1-18. [DOI: https://dx.doi.org/10.1080/1064119X.2022.2113485]

7. Liang, J.H.; Fan, Q.Y.; Qin, K. Influence of Karst Caves on the Pile’s Bearing Characteristics—A Numerical Study. Front. Earth Sci.; 2022; 9, 754330. [DOI: https://dx.doi.org/10.3389/feart.2021.754330]

8. Xu, H.H.; Liu, K.M.; Liu, W.L.; Zhu, C.; Zhu, S.T.; Dong, J.X. Genesis Analysis and Stability Evaluation of Karst Area in the Mudengdong Village Based on the Geological Investigation and Numerical Simulation Methods. Shock Vib.; 2022; 2022, 7494264. [DOI: https://dx.doi.org/10.1155/2022/7494264]

9. Li, L.; Zeng, Z.L.; Wang, Z.B.; Li, T.Y. Upper bound punching failure analysis of circular foundation overlying a cave in Karst area based on the Hoek-Brown failure criterion. J. Cent. South Univ. Sci. Technol.; 2020; 51, pp. 134-144. [DOI: https://dx.doi.org/10.11817/j.issn.1672-7207.2020.01.016]

10. Nie, Q.K.; Li, X.L.; Yuan, W.; Wang, A.L.; Wang, W. Calculation Method for the Critical Thickness of a Karst Cave Roof at the Bottom of a Socketed Pile. Adv. Mater. Sci. Eng.; 2021; 2021, 1669410. [DOI: https://dx.doi.org/10.1155/2021/1669410]

11. Lei, Y.; Li, P.J.; Liu, Z.Y.; Liu, J.Z. Method for calculation of buckling critical load of pile foundation crossing karst cave in karst area. Rock Soil Mech.; 2022; 43, pp. 3347-3356. [DOI: https://dx.doi.org/10.16285/j.rsm.2022.0377]

12. Zhao, M.; Chen, Y.; Xiao, Y.; Yang, C. Stability Analysis of Cave Roof under Rock-socketed Pile in Karst Area based on Three-hinged Arch Catastrophe Model. J. Disaster Prev. Mitig. Eng.; 2020; 40, 167. [DOI: https://dx.doi.org/10.13409/j.cnki.jdpme.2020.02.004]

13. Keawsawasvong, S.; Shiau, J. Stability of Spherical Cavity in Hoek-Brown Rock Mass. Rock Mech. Rock Eng.; 2022; 55, pp. 5285-5296. [DOI: https://dx.doi.org/10.1007/s00603-022-02899-4]

14. Liu, Z.Z.; Cao, P.; Lin, H.; Meng, J.J.; Wang, Y.X. Three-dimensional upper bound limit analysis of underground cavities using nonlinear Baker failure criterion. Trans. Nonferrous Met. Soc. China; 2020; 30, pp. 1916-1927. [DOI: https://dx.doi.org/10.1016/S1003-6326(20)65350-X]

15. Li, C.G.; Li, C.H.; Sun, C. Upper Bound Limit Analysis of Anisotropic Soils. CMC Comput. Mater. Contin.; 2020; 65, pp. 2607-2621. [DOI: https://dx.doi.org/10.32604/cmc.2020.04662]

16. Da Silva, M.V.; Deusdado, N.; Antao, A.N. Lower and upper bound limit analysis via the alternating direction method of multipliers. Comput. Geotech.; 2020; 124, 103571. [DOI: https://dx.doi.org/10.1016/j.compgeo.2020.103571]

17. Lyamin, A.V.; Sloan, S.W. Upper bound limit analysis using linear finite elements and non-linear programming. Int. J. Numer. Anal. Methods Geomech.; 2002; 26, pp. 181-216. [DOI: https://dx.doi.org/10.1002/nag.198]

18. Fang, Z.D.; Yang, W.M.; Wang, J. Study on the minimum safe thickness of water-proof rock mass in front of deep-buried tunnels. J. Cent. South Univ. Sci. Technol.; 2021; 52, pp. 2805-2816. [DOI: https://dx.doi.org/10.11817/j.issn.1672-7207.2021.08.026]

19. Liang, J.H.; Fan, Q.Y.; Liu, M.H.; Li, T.Y. Physical Model Test on the Influence of Karst Cave under and in front of Pile on the Stability of Embedded End of Antislide Pile. Math. Probl. Eng.; 2022; 2022, 6963982. [DOI: https://dx.doi.org/10.1155/2022/6963982]

20. Chen, H.Y.; Feng, Z.J.; Li, T.; Bai, S.F.; Zhang, C. Study on the vertical bearing performance of pile across cave and sensitivity of three parameters. Sci. Rep.; 2021; 11, 17342. [DOI: https://dx.doi.org/10.1038/s41598-021-96883-7]

21. Liu, Z.; Ma, H.; Mou, Y. Numerical simulation analysis and evaluation of stability of the karst foundation with developed joints and fissures. Carsologica Sin.; 2022; 41, pp. 100-110. [DOI: https://dx.doi.org/10.11932/karst20220105]

22. Li, Z.; Chen, Y.; Guo, Y.K.; Zhang, X.Y.; Du, S.G. Element Failure Probability of Soil Slope under Consideration of Random Groundwater Level. Int. J. Geomech.; 2021; 21, 04021108. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0002063]

23. Peng, P.; Li, Z.; Zhang, X.Y.; Liu, W.L.; Sui, S.G.; Xu, H.H. Slope Failure Risk Assessment Considering Both the Randomness of Groundwater Level and Soil Shear Strength Parameters. Sustainability; 2023; 15, 7464. [DOI: https://dx.doi.org/10.3390/su15097464]

24. Wu, B.; Lan, Y.; Yang, S.; Yang, J.; Pang, X. Study on Stability of Surrounding Rock Based on Strength Reduction Dynamic Analysis Method. Mod. Tunn. Technol.; 2020; 57, pp. 56-64. [DOI: https://dx.doi.org/10.13807/j.cnki.mtt.2020.03.009]

25. Dai, Z.H.; Fan, X.L.; Lu, C.J. Numerical analysis of stability of highway embankments and karst cave roofs in karst region. Rock Soil Mech.; 2014; 35, pp. 382-390. [DOI: https://dx.doi.org/10.16285/j.rsm.2014.s1.042]

26. Hu, Q.G.; Zhang, K.N.; Yang, J.S. Finite element analysis of ultimate bearing capacity of strip footing above Karst cave. J. Cent. South Univ. Sci. Technol.; 2005; 36, pp. 694-697.

27. Chang, S.P.; Zhang, S.M. Engineering Geology Manual; China Construction Industry Press: Beijing, China, 2006; 170.

28. Feng, Z.J.; Guo, S.Z.; He, J.B.; Huang, Z.Y.; Dong, Y.X.; Jiang, D.R. Study on the influence of karst size on the vertical bearing capacity of pile foundation. IOP Conf. Ser. Mater. Sci. Eng.; 2020; 780, 022005. [DOI: https://dx.doi.org/10.1088/1757-899X/780/2/022005]

29. Wang, J.; Zhang, D.; Sun, Z.; Fang, H.; Liu, C. Failure mechanism and scope prediction model of horizontal interbedded surrounding rock tunnel. Chin. J. Theor. Appl. Mech.; 2022; 54, pp. 2835-2849. [DOI: https://dx.doi.org/10.6052/0459-1879-22-133]