Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 825483 10.3389/feart.2021.825483 Earth Science Original Research Calculation Method for Optimizing the Bearing Capacity of Two-Layer Clay Foundation Based on Upper-Bound Limit Analysis Lin et al. Bearing Capacity Upper-Bound Limit Analysis Lin Xing-Chao 1 Zhang Qiang 1 * Wang Xianzhi 2 Feng Yelin 2 Zhu Guina 3 1 State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, Beijing, China 2 PowerChina Kunming Engineering Corporation Limited, Kunming, China 3 Beijing Glory PKPM Technology Co., Ltd., Beijing, China

Edited by: Dongxing Wang, Wuhan University, China

Reviewed by: Han Bao, Chang’an University, China

Yuyuan Chen, Kyushu University, Japan

 *Correspondence: Qiang Zhang, zhangq@iwhr.com

This article was submitted to Geohazards and Georisks, a section of the journal Frontiers in Earth Science

 16 02 2022 2021 9 825483 30 11 2021 27 12 2021 Copyright © 2022 Lin, Zhang, Wang, Feng and Zhu. 2022 Lin, Zhang, Wang, Feng and Zhu

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

In this study, a foundation bearing capacity limit analysis upper-bound solution (FLU) was proposed under the theoretical framework of limit analysis. To be more specific, an optimized mathematical model based on a rigid block discrete system was first established, with the minimum bearing capacity as the objective function, block velocity as the main variable, and the satisfaction of velocity compatibility, associated flow rule, and functional equilibrium equations of adjacent blocks as main constraints. Then, using an optimization approach, the upper-bound value of the foundation’s bearing capacity was obtained. On this basis, the principle of establishing the value interval of the most dangerous slip depth of the two-layer clay foundation was developed by investigating the effects of varying depths on the bearing capacity of the two-layer clay foundation. Meanwhile, an approach for calculating the bearing capacity of the two-layer clay foundation was proposed to achieve the goal of reaching the foundation’s minimum bearing capacity. Furthermore, using a mathematical example, the proposed approach was proven to be rational.

 bearing capacity of foundation limit analysis upper bound solution optimized mathematical model two-layer clay foundation 1 Introduction

In the field of geotechnical engineering, calculating the ultimate bearing capacity of foundations has always been a great concern (Terzaghi and Peck 1967; Davis and Booker 1974; Chen 1975; Griffiths 1982; Michalowski and Lei 1996), with significant implications for underground space and engineering safety. Upper-bound limit analysis based on plastic mechanics is a powerful approach for calculating the ultimate bearing capacity of foundations (Lyamin and Sloan 2002; Huang and Qin 2009; Osman 2019; Shamloo and Imani 2020). It is primarily composed of the limit analysis finite element and the limit analysis approach based on discretization of a rigid block in terms of various discretization approaches used in the calculation model. This method does not need to introduce too many assumptions, has a strict theoretical basis, and has advantages in the calculation of foundation bearing capacity.

Recently, many researchers have proposed the limit analysis finite element upper-bound solution approach of foundation bearing capacity by combining the plastic mechanics of limit analysis theory and finite element (Alkhafaji et al., 2020; Pham et al., 2020; Shamloo and Imani 2020). The approach based on the upper-bound theorem of plastic mechanics can compensate for the traditional limit equilibrium approach’s imprecise theoretical foundation. However, it can lead to heavy workload and low calculation efficiency since the finite element discretization must be conducted on the research object, and the velocity of each element node is considered an unknown quantity. Furthermore, a large number of unknown quantities and a high degree of freedom make the calculation process difficult. Overall, its widespread adoption and application in complex practical engineering will be challenging (Wang et al., 2019; Wang et al., 2020).

In the foundation bearing capacity limit analysis approach based on the discretization of rigid blocks, on the premise that the velocity of adjacent block meets the velocity compatibility and the interface velocity satisfies the relation equation of the associated flow principle, the velocity field of the sliding body is first obtained by recursing the velocity field of the slip mass in recurrence, and then the ultimate load is solved in accordance with the functional equilibrium equation that the internal energy dissipation is equal to the external work. The approach also implies the assumption of tangential relative velocity direction, and it can deal with the problem of introducing a significant number of assumptions in the traditional limit equilibrium approach.

In this study, an optimized mathematical model for the foundation bearing capacity was directly established in accordance with the upper-bound theorem instead of adopting the thought of introducing a significant number of assumptions in the traditional limit equilibrium approach. In this way, the foundation bearing capacity problem can be transformed into an optimization solution problem. Based on this, a foundation bearing capacity limit analysis upper-bound solution (FLU) was proposed to solve the foundation bearing capacity using a nonlinear mathematical programming solver. Meanwhile, the principle of determining the value interval of the most dangerous slip depth of the two-layer clay foundation was proposed by studying the effects of various depths on the bearing capacity of the two-layer clay foundation. Similarly, the calculation approach of iterative optimization using the dichotomy was proposed for the bearing capacity of the two-layer clay foundation in order to achieve the goal of reaching the minimum foundation bearing capacity. Furthermore, the proposed method was proven to be rational using a mathematical example.

 2 Foundation Bearing Capacity Limit Analysis Upper-Bound Solution (FLU) 2.1 Discretization of Rigid Block

As shown in Figure 1, this study proposes an improved rigid block discrete approach that can automatically generate a block division model based on geological conditions and failure mode of the actual foundation. Through the interface segment or block set that constitutes the block division system, the block division system model can be defined as the following: { V = L 1   L 2 … L j … L h l − 1 L h l V = B 0 B 1   B 2 … B k … B h b − 1 B h b where j is the index number of the interface segments in the block division system; h l is the total number of interface segments in the block division system; h l ′ is the total number of all interface segments, excluding the ground surface; k is the index number of the blocks in the block division system; h b is the total number of blocks contained in the slip mass; B 1 ∼ B h b are the failure mechanisms; and B 0 indicates the soil without plastic failure.

Improved rigid block division model.

The normal vector in the interface segment is defined as facilitating the establishment of an optimized mathematical model of the limit analysis upper-bound solution in the block division system. In addition, the inner normal vector’s calculation approach is briefly introduced below.

Any block B k in the block division system can be formed by a series of nodes or interface segments arranged in sequence, as demonstrated below in Eq. 2 { B K = a k 1 , a k 2 … a k i … a k h 1 − 1 , a k h 1 B K = a k 1 a k 2 … a k i a k i + 1 … a k h 1 − 1 a k h 1 where a k i are the vertices that compose the polygon B K and arranged clockwise, and the last point of the polygon B K coincides with the first point, or a k h 1 = a k 1 , ensuring the polygon’s closure; i is the index number of the vertex; h 1 is the total number of vertices; and a k i a k i + 1 is the interface segment.

The internal normal vector of the a k i a k i + 1 interface line segment is obtained using Eq. 3i. m k i = ( a k i + 1 − a k i ) × ( a k i − a k i − 1 ) × ( a k i + 1 − a k i ) / | m k i | where m k i is the a k i a k i + 1 internal normal vector of the interface segment in the block B K . The normal vector calculated using the above equation in the boundary segment always points to the inner part of the polygon, irrespective of the vertex’s rotation direction, and is unique in geometric space. Furthermore, as shown in Figure 2, the internal normal vectors of the interface shared by two adjacent blocks are equal in magnitude and opposite in direction.

Definition of the normal vector in the interface segment.

 2.2 Variable System

Variables involved in the optimization model of the foundation bearing capacity upper-bound solution consist of the ultimate bearing capacity (P), the velocity V B k ( V B x k , V B y k ) of each block in the block division system, and the intermediate variable of the relative velocity V j ( V T j , V N j ) of the interface segment introduced for simplifying the optimization solution. To describe the directional characteristics of the interface force, a local coordinate system is established on all interfaces, with the direction of the segment as the positive axis and the direction vector as t j . As shown in Figure 3, the N axis is perpendicular to the boundary segment, with the direction vector of n j = l z × t j , where l z = ( 0,0,1 ) . Then, the relative velocity of the interface can be expressed using Eq. 3ii: { V T j = V T j ∗ t j V N j = V N j ∗ n j

Approaches for representing the force on the interface used by two adjacent blocks in common.

In the block, the direction of the local coordinate system in B k should be determined using the internal normal vector of the interface segment L j , as shown in Eq. 4: { V T k i = V T j ∗ t j ∗ ( n j · m k i ) V N k i = V N j ∗ n j ∗ ( n j · m k i )

When n j and m k i are in the same direction, then n j · m k i = 1 . When n j and m k i are in the opposite direction, then n j · m k i = − 1 . According to this definition, the local coordinate direction vectors established on the common interface l 12 are opposite to each other for two adjacent blocks (such as B 6 and B 8 in Figure 3). In that case, only a set of variables V 12 ( V T 12 , V N 12 ) on the common interface l 12 are required for describing the characteristics of acting force between blocks equal in magnitude and opposite in direction.

Obviously, the variables required by the variable system established by this method are less than the finite element of limit analysis. The main reasons are as follows: 1) the number of blocks after block discretization is less than the number of finite element elements; 2) the method in this paper is based on the block interface, the limit analysis finite element is based on nodes and needs to consider the situation of overlapping nodes, and the number of block interfaces is less than the number of nodes.

 2.3 Constraint equations

According to the upper-bound theorem of limit analysis, the constraint equations comprise associated flow rule, velocity compatibility, and functional equilibrium equations.

2.3.1 Associated Flow Rule

The interface’s relative velocity must also conform to the associated flow rule, or the angle between the relative kinematic velocity and the interface should be equal to the internal angle φ j of friction: V N j = tan ⁡ φ j | V T j |

To better describe the characteristics of | V T j | and reduce the difficulty of finding the optimization solution, two intermediate variables, V T a j and V T b j , are introduced and defined as the following: { | V T j | = V T a j + V T b j V T j = V T a j − V T b j V T a j > 0 V T b j > 0 V T a j V T b j = 0

Equation 5 can be rewritten as the following: V N j = tan ⁡ φ j ( V T a j + V T b j )

The upper-bound solution proposed by Chen et al. (2003), for example, assumed that the direction of V T j V T j on the bottom slip surface is determined. However, the rationality of the assumption of the V T j direction on the interface between blocks cannot be guaranteed in some complicated engineering problems. Furthermore, unreasonable direction assumptions can lead to inaccurate results. The direction of V T j in the proposed approach is determined directly through optimization rather than making assumptions, which is theoretically more rigorous.

 2.3.2 Velocity Compatibility Condition

Adjacent blocks should conform with the velocity compatibility condition, that is, a velocity triangle is made up by the relative velocity of the interface velocity and the velocity of two adjacent blocks and can be expressed mathematically as the following: { V N j = ( V j 1 − V j 2 ) · n j V T j = V T a j − V T b j = ( V j 1 − V j 2 ) · t j where V j 1 and V j 2 are the blocks bounded by the common interface segment l j . In V j 1 , the internal normal vector of the internal interface segment l j is m j 1 = n j . In V j 2 , the internal normal vector of the internal interface segment l j is m j 1 = − n j . Regarding the boundary segment, V j 2 = V B 0 = ( 0,0,0 ) on the slip surface, V j = V j 1 , as shown in Figure 4.

Velocity compatibility diagram.

 2.3.3 Functional Equilibrium Equation

According to the associated flow law, the internal energy dissipation on the interface polygon L j can be expressed as: E 1 = ∑ j = 1 h l ′ | V T j | c j l j = ∑ j = 1 h l ′ ( V T a j + V T b j ) c j l j where c j and l j are the cohesive force and length of the interface segment L j , respectively.

The gravity acting of the block is E 2 = ∑ k = 1 h b V k · G d k where G is the gravity direction vector, taking (0, −9.81, 0) in general, and d k is the density of the block, B k .

The external work is E 3 = ∑ i = 1 h q Q i · V i where Q i is the external load that includes the unknown variables such as foundation bearing capability ( P ) and overlying load q ; i is the index number of external load; hq is the total number of external loads; and V i is the velocity of the block where Q i is located.

On this basis, a functional equilibrium equation for the whole slip mass can be obtained: ∑ j = 1 h l ' ( V T a j + V T b j ) c j l j = ∑ k = 1 h b V k · G d k + ∑ i = 1 h q Q i · V i

2.3.4 Other Constraints

To ensure that the blocks are not embedding, the relative normal velocity must also satisfy the constraint equation of Eq. 13: V N j ≥ 0

2.4 Optimization Model

According to the upper-bound theorem, the foundation bearing capacity calculated by any failure mode satisfying the constraints is greater than the real foundation bearing capacity. The smaller the calculated foundation bearing capacity is, the closer it is to the real solution. Therefore, the objective function of the optimization model is to minimize the bearing capacity Eq. 14. The objective function Eqs 5, 6, 9, 14 as standard optimization models of constraint conditions are shown in Figure 5. F L U = M i n ( P )

Optimization model.

 2.5 Verification of Equation for Foundation Bearing Capacity

For foundations with agravic medium, Prandtl–Reissner assumes that the foundation bed is completely smooth with the failure mode (as shown in Figure 6A) and the calculation equation for ultimate bearing capacity given, using the characteristic approach: { P = q N q + c N c N q = tan 2 ( π 4 + φ 2 ) e π ⁡ tan ⁡ φ N c = ( N q − 1 ) cot ⁡ φ where N q is the Earth pressure coefficient for the overlying load q ; N c is the Earth pressure coefficient; and c is the soil cohesion.

Prandtl–Reissner failure mode and rigid block discrete. (A) Prandtl–Reissner failure mode. (B) Rigid block discrete models of φ = 0 ° and φ = 30 °

As shown in Figure 6B, the rigid block discrete model is generated automatically according to the geometric relationship in the failure mode (when the friction angles in the soil mass are 0° and 30°). Then, an optimized mathematical model is established to obtain the ultimate bearing capacity of the foundation with an agravic medium under various parameters (Table 1). For the simple foundation with agravic medium, the calculation results obtained by the FLU proposed in this study are consistent with the theoretical solution, proving the rationality of the proposed approach.

Comparisons of ultimate bearing capacity and characteristic parameters under different internal friction angles.

φ ( ° )
0 5 10 15 20 25 30 35
P 1 c b ( q = 0 ) ① 5.14 6.49 8.34 10.98 14.83 20.73 30.15 46.14
② 5.14 6.49 8.34 10.98 14.83 20.72 30.14 46.12
P 2 c b ( q = 40 ) ① 5.54 7.12 9.33 12.55 17.40 25.00 37.51 59.47
② 5.54 7.12 9.33 12.55 17.39 24.99 37.50 59.44

①: FLU (this study). The cup practice of calculating the foundation bearing capacity by this method is less than 1 s.

②: Prandtl–Reissner.

 3 Research on the Value Range of Slip Depth of Two-Layer Clay Foundation

The bearing capacity of the two-layer clay foundation is difficult to calculate. Many researchers have conducted analytical studies on it. To better present the advantages of the rigid block discrete in this study, slip modes similar to those used by Prandtl–Reissner and Terzaghi are for the discretion of rigid body, as shown in Figure 7. The slip depth coefficient H / b in Figure 7 is the layer depth; h / b is the slip depth. An accurate block discrete model can be obtained according to the characteristics of the slip mode by determining h / b and automatically generated using the proposed approach; c 1 / c 2 is the ratio of cohesive forces of the upper and lower soil; and h b is the number of blocks.

Calculation model of bearing capacity of the two-layer clay foundation.

To determine the value range of slip depth for two-layer clay foundations, the effect of slip depth h b on foundation bearing capacity is investigated under different combinations of c 1 c 2 = 1 , c 1 c 2 < 1 ( 0.2 , 0.25 , 0.33 , 0.5 ) , c 1 c 2 > 1 ( 2,3,4,5 ) , and H b = 0.2 、 0.5 、 1.0 、 1.5 , as shown in Figure 8. According to the calculation results:

1) When c 1 c 2 = 1 , the problem of calculating the bearing capacity of two-layer clay foundations is reduced to the problem of calculating the bearing capacity of uniform and non-cohesive foundations. The coefficient P c 1 b of the foundation bearing capacity first decreases and then increases with h b , reaching a minimum value of 5.14 at h b = 0.71 , which is the Prandtl theoretical solution.

2) When c 1 c 2 < 1 and H b < 0.71 ( H b = 0.2 、 0.5 ) , local failure occurs in the upper soil due to its low strength.

3) When c 1 c 2 < 1 and H b > 0.71 ( H b = 1.0 、 1.5 ) , the variation law before reaching the layered interface H b is the same as that of the uniform foundation, and P c 1 b reaches 5.14, the minimum value, at h b = 0.71 . No low value of P c 1 b can be obtained within the range of the lower soil because the strength of the lower soil is greater than that of the upper soil as h b increases. In that case, it can be considered a uniform foundation.

4) When c 1 c 2 > 1 and H b < 0.71 ( H b = 0.2 、 0.5 ) , the bearing capacity coefficient of the foundation first decreases and then increases with h b . Furthermore, the minimum value will be lagging due to the lower strength of lower soil, h b > 0.71 .

5) When c 1 c 2 > 1 and H b < 0.71 ( H b = 1.0 、 1.5 ) , a downward and upward trend can be seen on both sides of H b , with two distinct troughs. In the actual search, P c 1 b values corresponding to two troughs should be compared, with the smaller one being used as the ultimate load.

h b − P c 1 b . Curve. (A) c 1 c 2 = 1 , (B) c 1 c 2 > 1 , (C) c 1 c 2 > 1 .

According to the value range of h b provided by the above results (Table 2), reducing the search range while calculating the ultimate load can enhance the calculation efficiency.

Range of slip depth.

c 1 c 2 = 1 c 1 c 2 < 1 c 1 c 2 > 1
H b < 0.71 H b > 0.71 H b < 0.71 H b > 0.71
h b Slip depth h b = 0.71 h b = H b h b = 0.71 h b > 0.71 h b = 0.71 or H b < h b
4 Optimized Calculation Method for the Bearing Capacity of Two-Layer Clay Foundation

Based on Table 2 and FLU, an optimized calculation approach for the bearing capacity of the two-layer clay foundation is proposed. The proposed approach first calculates the range of slip depth according to Table 2 and then obtains the final slip depth and foundation bearing capacity using the dichotomy in iteration. Figure 9 shows the specific process.

Optimized calculation approach for the bearing capacity of the two-layer clay foundation. (A) Flowchart, (B) Iternative solution using the dichotomy

As shown in Figure 9, L (min, max) is an optimization model of an iterative solution using the dichotomy, where min is the minimum of the initial iteration interval, and max is the maximum of the initial iteration interval. FLU ( h b ) is the ultimate bearing capacity corresponding to the slip depth h b (as obtained by the upper-bound solution of the limit analysis), FLU ( Δ h b ) = FLU ( h b ) − FLU ( h b + Δ ) , and Δ = 0.001 is the calculation analysis increment.

Different calculation approaches of optimizing the bearing capacity of the two-layer clay foundation are used to calculate ultimate load under various combinations of c 1 c 2 and H b according to the foundation bearing capacity of the two-layer clay foundation. Table 3 shows the ultimate load results and literature comparison. Figure 10 shows the typical failure modes obtained by the search. As can be observed from the calculation results, when c 1 c 2 < 1 and H b < 0.71 ( H b = 0.2 、 0.5 ) , local failure occurs in the upper soft soil, and the calculation result is between the Terzaghi solution (5.71) and Prandtl solution (5.14) of uniform and agravic clay, which is supported by the literature results. When c 1 c 2 < 1 and H b > 0.71 ( H b = 1.0 、 1.5 ) , the upper- and lower-bound solutions calculated using the proposed calculation approach are 5.14, which is within the range of the upper- and lower-bound solutions in the literature. The proposed approach’s rationality has been proven to a certain extent. When c 1 c 2 > and H b < 0.71 ( H b = 0.2 、 0.5 ) , the slip depth increases and the ultimate bearing capacity decreases rapidly; the corresponding ultimate bearing capacity coefficient P c 1 b corresponding to H b = 0.2 and c 1 c 2 = 5 is 1.80. When c 1 c 2 > and H b < 0.71 ( H b = 1.0 、 1.5 ) , there are two possible instability modes. One is that the buried depth of the soil is deep in the lower soft soil (such as H b = 1.5 ) and not much different from the upper soil (such as c 1 c 2 = 2 、 3 、 4 ). The lower soil does not influence the ultimate bearing capacity. In that case, the uniform foundation is considered with the ultimate load factor P c 1 b = 5.14 and slip depth h b = 0.71 . The other is that the weak lower soil affects the ultimate bearing capacity with deep slip generated. When H b = 1.5 and c 1 c 2 = 5 , the slip depth h b reaches 2.93.

Statistics of ultimate bearing capacity coefficient.

H/b — c 1 / c 2
0.2 0.25 0.33 0.5 1 2 3 4 5
0.2 This paper h/b 0.2 0.2 0.2 0.2 0.71 0.91 1.05 1.18 1.29
Upper bound 5.72 5.72 5.72 5.72 5.14 3.17 2.44 2.06 1.81
Merifield et al. (1999) Lower bound 5.44 5.44 5.44 5.44 4.86 3.08 2.15 1.75 1.48
Upper bound 5.89 5.89 5.89 5.89 5.32 3.24 2.44 2.00 1.73
0.5 This paper h/b 0.5 0.5 0.5 0.5 0.71 1.13 1.40 1.62 1.80
Upper bound 5.72 5.72 5.72 5.72 5.14 3.89 3.28 2.94 2.70
Merifield et al. (1999) Lower bound 4.86 4.86 4.86 4.86 4.86 3.52 2.84 2.44 2.16
Upper bound 5.31 5.31 5.31 5.31 5.31 3.89 3.16 2.74 2.44
1.0 This paper h/b 0.71 0.71 0.71 0.71 0.71 1.39 1.80 2.13 2.41
Upper bound 5.14 5.14 5.14 5.14 5.14 4.92 4.70 4.20 3.98
Merifield et al. (1999) Lower bound 4.94 4.94 4.94 4.94 4.94 4.44 3.89 3.46 3.10
Upper bound 5.32 5.30 5.30 5.30 5.30 4.82 4.24 3.83 3.54
1.5 This paper h/b 0.71 0.71 0.71 0.71 0.71 0.71 0.71 0.71 2.93
Upper bound 5.14 5.14 5.14 5.14 5.14 5.14 5.14 5.14 5.12
Merifield, Sloan et al. (1999) Lower bound 4.94 4.94 4.94 4.94 4.94 4.87 4.69 4.24 3.89
Upper bound 5.30 5.30 5.30 5.30 5.32 5.31 5.15 4.84 4.56

Typical failure modes.

 5 Conclusion

1) An optimized mathematical model is first established based on a rigid block discrete system with the minimum bearing capacity as the objective function, block velocity as the main variable, and the satisfaction of velocity compatibility, associated flow rule, and functional equilibrium equations of an adjacent block as main constraints in the rigid block discrete system. Then, FLU is proposed after the upper-bound value of the bearing capacity of the foundation is obtained through an optimization solution.

2) The influence of slip depth ( h b ) on the foundation bearing capacity is investigated for the complex two-layer clay foundation before proposing a value range of slip depth under different parameter combinations ( c 1 c 2 ) and varied layer depths ( H b ) .

3) Based on the above results, an optimized calculation approach for the bearing capacity of the two-layer clay foundation is proposed by introducing an iterative solution using a dichotomy. Furthermore, the approach’s rationality is verified using a calculation example.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

 Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

 Funding

This research was financially supported by the National Key R&D Program of China (No. 2018YFC0407000), the National Natural Science Foundation of China (Nos. 51809289, U1965204) and the IWHR Research and Development Support Program (Nos. GE0199A082021, GE110145B0022021). Research Project of China Three Gorges Corporation (Contract No. JG/19055J).

 Conflict of Interest

XW and YF were employed by the company PowerChina Kunming Engineering Corporation Limited, and GZ was employed by the company Beijing Glory PKPM Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Alkhafaji H. Imani M. Fahimifar A. (2020). Ultimate Bearing Capacity of Rock Mass Foundations Subjected to Seepage Forces Using Modified Hoek–Brown Criterion. Rock Mech. Rock Eng. 53 (1), 251–268. 10.1007/s00603-019-01905-6 Chen J. Yin J.-H. Lee C. F. (2003). Upper Bound Limit Analysis of Slope Stability Using Rigid Finite Elements and Nonlinear Programming. Can. Geotech. J. 40(4), 742–752. Chen W. F. (1975). Limit Analysis and Soil Plasticity. New York: Elsevier Scientific Publishing Company. Davis E. H. Booker J. R. (1974). The Effect of Increasing Strength with Depth on the Bearing Capacity of Clays : 14F, 6R. GEOTECHNIQUE, V23, N4, DEC. 1973, P551–563. Int. J. Rock Mech. Mining Sci. Geomechanics Abstr. 11 (3), 62. 10.1016/0148-9062(74)91688-x Griffiths D. V. (1982). “Computation of Bearing Capacity on Layered Soils,” in Proc., 4th Int. Conf. on Num Meth. In Geomechanics, Vol. 1, 163–170. Huang M. Qin H.-L. (2009). Upper-bound Multi-Rigid-Block Solutions for Bearing Capacity of Two-Layered Soils. Comput. Geotechnics 36 (3), 525–529. 10.1016/j.compgeo.2008.10.001 Lyamin A. V. Sloan S. W. (2002). Upper Bound Limit Analysis Using Linear Finite Elements and Non-Linear Programming. Int. J. Numer. Anal. Methods Geomech. 26 (2), 181–216. Merifield R. S. Sloan S. W. Yu H. S. (1999). Rigorous Plasticity Solutions for the Bearing Capacity of Two-Layered Clays. Géotechnique 49 (4), 471–490. 10.1680/geot.1999.49.4.471 Michalowski R. L. Lei S. (1996). Bearing Capacity of Footings over Two-Layer Foundation Soils. J. Geotechnical Eng. 121 (5), 421–428. 10.1061/(asce)0733-9410(1995)121:5(421) Osman A. S. (2019). Upper Bound Solutions for the Shape Factors of Smooth Rectangular Footings on Frictional Materials. Comput. Geotechnics 115 (Nov), 103177. 10.1016/j.compgeo.2019.103177 Pham Q. N. Ohtsuka S. Isobe K. Fukumoto Y. (2020). Limit Load Space of Rigid Footing under Eccentrically Inclined Load. Tokyo: Soils and Foundations. Shamloo S. Imani M. (2020). Upper Bound Solution for the Bearing Capacity of Rock Masses Considering the Embedment Depth. Ocean Eng. 218 (6), 108169. 10.1016/j.oceaneng.2020.108169 Terzaghi K. Peck R. B. (1967). Soil Mechanics in Engineering Practice. second edition. New YorK: Wiley. Wang X. G. Lin X. C. Sun P. Li X. Ling Y. Y. (2020). 2D Slope Stability Analysis Based on Pan's Maximum Principle. Bull. Eng. Geology. Environ. 79 (8). 10.1007/s10064-020-01840-9 Wang X. Lin X. Li X. Sun P. Ling Y. (2019). 3D Slope Stability Analysis Method Based on Pan's Maximum Principle. Landslides 17 (1). 10.1007/s10346-019-01303-7