Undrained Stability of Unsupported Rectangular Excavations: Anisotropy and Non-Homogeneity in 3D

Abstract

The stability of unsupported rectangular excavations in undrained clay is examined under the influence of anisotropy and heterogeneity using the three-dimensional finite element upper and lower bound limit analysis with the Anisotropic Undrained Shear (AUS) failure criterion. Three anisotropic undrained shear strengths are considered in the study, namely triaxial compression, triaxial extension, and direct simple shear. Special considerations are given to the study of the linearly-increased anisotropic shear strengths with depth. The numerical solutions are presented by an undrained stability number that is a function of four dimensionless parameters, i.e., the excavated depth ratio, the aspect ratio of the excavated site, the shear strength gradient ratio, and the anisotropic strength ratio. To the authors’ best knowledge, this is the first of its kind to present the stability solutions of 3D excavation considering soil anisotropy and heterogeneity. As such, this paper introduces a novel approach for predicting the stability of unsupported rectangular excavation in undrained clays in 3D space, accounting for soil anisotropy and non-homogeneity. Notably, it develops a basis to formulate a mathematical equation and design charts for estimating the stability factor of such type of excavation, which should be of great interest to engineering practitioners.

1. Introduction

Unsupported excavation does not require retaining wall systems, and it is considered one of the affordable construction methods that are widely employed in many shallow underground construction projects. Shallow underground structures such as pipelines, shallow tunnels, and underpasses can be constructed by utilizing this excavation technique. Other examples may include the construction of piers, footings, retaining structures, raft foundations, mat foundations, and water tanks. An unsupported excavation during construction, if not properly assessed, can lead to an eventual collapse of the excavation wall that could result in an injury or fatality. These unfortunate events can cost money and cause death. It is, therefore, imperative to assess the stability of such unsupported excavations to reduce the risk of soil failure, thereby improving site safety and preventing death. This study aims to contribute to reducing that risk by providing a novel approach that predicts the undrained stability of unsupported rectangular excavations in anisotropic and non-homogeneous clays.

In general, the excavation can have either cylindrical, conical, or rectangular shapes. Griffiths and Koutsabeloulis [1] used a displacement-based elastoplastic finite element analysis to study the stability of cylindrical excavations under axisymmetric conditions. The same problem was also examined by Britto and Kusakabe [2,3] using the plastic-bound theorems. The recent development of finite element limit analysis (FELA) is a powerful numerical method based on the lower bound (LB) theorem, the upper bound (UB) theorem, and the finite element technique, as demonstrated in [4,5,6,7,8,9,10]. The axisymmetric FELA was employed by [11,12,13,14,15,16] to obtain stability solutions for vertical circular excavations. Recently, the stability of unsupported conical excavations was investigated by [17,18,19,20,21].

Among the various shapes of excavation, rectangular and cylindrical shapes are the most common in practice. Although cylindrical shapes may lead to smaller amounts of excavated material and they are more stable due to the arching effect [22], rectangular excavation is more widely used because it is less complex to build due to its shape and it follows a similar shape to most common subsoil structures being built within the excavation (e.g., footings, pile caps, piers, mat foundations). For the problem of unsupported rectangular excavations, stability solutions were reported by Ukritchon et al. [22] using 3D FELA. Their solutions are based on the Tresca failure criterion, which is limited to isotropic clays. However, it is common knowledge that soil, particularly clays, normally exhibits anisotropy and heterogeneity due to depositional geologic processes. It is generally recognized that soil anisotropy can have a substantial influence on clay stability, e.g., [23,24,25,26]. Ladd [23,24] reported that partial strength anisotropy in natural clays is generated through the processes of deposition and sedimentation with favored particle orientation. It was also demonstrated in the same paper that the anisotropic shear strengths of clays are very much dependent on the different shearing modes as well as the depositional axis. Thus, including anisotropy and non-homogeneity in the stability solution of unsupported excavation will provide a more reliable and realistic solution to excavation problems. Some studies have explored the problem of excavation in anisotropic clays, e.g., [27,28], but most are braced or supported.

Indeed, there are three undrained shear strengths that can be obtained in a laboratory: (1) triaxial compression (TC), (2) triaxial extension (TE), and (3) direct simple shear (DSS). The three undrained shear strengths have contributed to the development of mathematical forms of failure criteria for anisotropic soils, e.g., [9,10,25,29,30,31,32,33]. Recently, Krabbenhoft and Lyamin [34] developed a unique failure criterion for anisotropic clays, known as AUS (Anisotropic Undrained Shear), by adopting the Generalized Tresca (GT) criterion for undrained total stress analysis. Even though both Davis and Christian’s (DC) failure criterion [29] and the AUS failure criteria consider an empirical correlation of the undrained strength (s_u) of clay in triaxial compression (TC), direct simple shear (DSS), and triaxial extension (TE), the explicit form of the DC failure criterion cannot be applied to 3D problems since it was developed under plane strain condition. Unlike the DC model, the AUS model was developed under 3D coordinates, which can be used to simulate 3D problems in the Cartesian coordinates. As a result, the AUS model is preferred in this paper to investigate the stability of 3D unsupported excavations.

The AUS model has recently been included in the 3D FELA software, OptumG3 [35], and it has been successfully applied to the stability problems of plate anchors [36] and caissons [37]. Apart from the recent AUS studies, FELA has previously been adopted to report numerical results for various 3D geotechnical problems, such as determining the capacity of a rigid pile with a pile cap in Zhou et al. [38], the trapdoor stability problem in Shiau et al. [39], the bearing capacity of footings on slopes in Yang et al. [40], and the tunnel stability problem in Shiau and Al-Asadi [41,42,43].

A thorough search of the relevant literature shows that the undrained stability numbers for unsupported rectangular excavations considering both anisotropy and heterogeneity have never been reported in the literature. The most recent paper by Yodsomjai et al. [44], which has close similarities to this current study, tackled the undrained stability of unsupported conical slopes in anisotropic clays, which was similarly analyzed using the AUS failure criterion. However, due to its axisymmetric condition, it becomes a 2D plane strain problem rather than 3D. Other than the study undertaken by Ukritchon et al. [22] on the 3D undrained stability of unsupported rectangular excavations in non-homogeneous clays, which is also similar to the current study but without considering soil anisotropy, most of the other studies in the literature that dealt with anisotropy and heterogeneity in clays were related to other stability problems such as trapdoors [30], pile bearing capacity [32], unlined square tunnels [33], anchors [36], and suction caissons [37]. Therefore, the aim of this paper is to study this underexplored subject on the 3D undrained stability of unsupported rectangular excavations in clays with linearly increasing anisotropic shear strength. The stability solutions were formulated by a dimensionless stability number that is a function of four dimensionless parameters: the excavated depth ratio, the aspect ratio of the excavated site, the shear strength gradient ratio, and the anisotropic strength ratio. The selected failure mechanisms of this problem were examined to demonstrate the effects of all four dimensionless parameters. With the development of accurate design equations, the study would assist practicing engineers in determining the soil stability of unsupported rectangular excavations in clays with anisotropy and heterogeneity.

2. Statement of the Problem and Modelling Technique

Figure 1 shows the problem of defining a 3D unsupported rectangular excavation. Due to the problem of symmetry, only a quarter of the model domain was used in the analysis. See Figure 1a for the model. The excavation depth is denoted by H, B is the excavation width, and L is the length.

The AUS failure criterion with the associated flow rule was used to study the 3D soil stability of the unsupported rectangular excavation. The three anisotropic undrained shear strengths obtained from triaxial compression (s_uTC), triaxial extension (s_uTE), and direct simple shear (s_uDSS) were the required strengths for this failure criteria. According to Krabbenhoft et al. [45], two anisotropic strength ratios can be defined using the three undrained shear strengths: (1) r_e = s_uTE/s_uTC and (2) r_s = s_uDSS/s_uTC. The relationship between r_e and r_s is the harmonic mean, which can be written as follows:

$$r_s=\frac{2r_e}{1+r_e}$$

(1)

As shown in Equation (1), the parametric analysis only used one anisotropic strength ratio, which is r_e. Note that r_s is a function of r_e, and the range of r_e should be between 0.5 and 1. A change in the r_e value may vary the AUS failure criterion’s failure surface, as shown in Figure 2 [34,45]. The form of the yield function of the AUS model with the harmonic mean of three undrained shear strengths can be expressed by Equation (2):

$$F_u={\sigma }_1-{\sigma }_3+\left(r_e-1\right)\left({\sigma }_2-{\sigma }_3\right)-2s_{uTC}=0$$

(2)

where σ₁ ≥ σ₂ ≥ σ₃ are the principal stresses (positive in compression), and F_u is the yield function. It should be noted that the AUS failure criterion becomes the Tresca failure criterion when r_e = 1, meaning the isotropic state, i.e., s_uTC = s_uTE = s_uDSS. Note that, for the AUS failure criterion, three undrained shear strengths were considered to be an empirical correlation of the undrained strength (s_u) of clay in triaxial compression (TC) for s_uTC, direct simple shear (DSS) for s_uDSS, and triaxial extension (TE) for s_uTE.

The increasing shear strength with depth, i.e., heterogeneous soil behaviors, is another important factor when determining soil stability. This variation in shear strength has been considered by many researchers for the problems of the face stability of tunnels [10,31,33,46], supported excavations [47], piles [48,49], floodwalls [50], and active trapdoors [30,51]. This study considered three anisotropic undrained shear strengths that linearly increase with depth. Mathematically, they are expressed in Equations (3)–(5).

$$s_{uTC}\left(z\right)=s_{uTC0}+\rho z$$

(3)

$$s_{uTE}\left(z\right)=s_{uTE0}+r_e\rho z$$

(4)

$$s_{uDSS}\left(z\right)=s_{uDSS0}+r_s\rho z$$

(5)

where s_uTC0, s_uTE0, and s_uDSS0 are the anisotropic undrained shear strengths at the ground level, z is the depth from the ground surface, and ρ is the linear strength gradient. See Figure 1b for the linear distributions of the three anisotropic undrained shear strengths.

Using the dimensional analysis [52], a function combining four dimensionless parameters that are variables of a stability number function can be expressed by Equation (6).

$$N=\frac{\gamma H}{s_{uTC0}}=f(\frac{B}{L},\frac{H}{B},r_e,m=\frac{\rho B}{s_{uTC0}})$$

(6)

where N is the stability number, B/L is the aspect ratio of the excavated site, H/B is the excavated depth ratio, r_e is the anisotropic strength ratio, and m is the strength gradient ratio.

In the lower bound analysis, a four-node tetrahedron element is used, where six unknown nodal stresses are used for each node of tetrahedral elements. The statically admissible stress discontinuities are allowed to produce the continuity of normal and shear stresses along the interfaces of all the elements. The conditions of stress equilibrium, stress boundary condition, and the AUS failure criterion are all constraints in a typical LB analysis, in which the objective function is to maximize the critical unit weight γ that yields an excavation collapse. In the upper bound theorem, a four-node tetrahedron element is also adopted for the upper bound analysis, where each node contains three unknown velocities that vary linearly within the tetrahedron element. The kinematically admissible velocity discontinuities are applied at the interfaces of all the elements. The material is set to obey the AUS failure criterion associated flow rule. The formulated objective function is to minimize the critical unit weight γ. The obtained critical γ from both LB and UB analyses were then used to compute the stability number in Equation (6). More details on the LB and UB FELA can be found in [5].

Figure 3 presents a typical 3D FELA mesh used for the analysis. The nodes around the sides of the model are fixed in the normal direction to the planes of the sides. The same boundary condition is applicable to the two symmetrical planes as well. At the bottom domain, the nodes are fixed in all directions. Both the ground surface and the excavation faces are free to move in all directions. The overall domain size is chosen to be sufficiently large such that the stability solutions are not affected by the boundary conditions, i.e., the effects of boundary size on the computed LB solutions are minimized by generating LB meshes with sufficient lateral and lower dimensions that produce a computed plastic yielding zone that does not intersect the boundary planes. Automatic adaptive mesh refinement is one of the advanced features of the 3D program. This technique is based on Ciria et al. [53], where the numbers of elements in sensitive zones (i.e., with very high plastic shear strains) are increased through successive iterations with adaptive mesh refinement. The required input for the adaptive scheme is the original and target number of elements, the number of adaptive iterations, and the control variable for error estimation (i.e., shear power in this paper). In this study, 5000 initial elements were employed, which was expanded to 10,000 elements after five iterations.

Note that, the range of four dimensionless parameters in all studies of the paper are: (1) H/B = 0.5, 1, 2, 3, 4; (2) B/L = 1, 2/3, 1/2, 1/4, 1/8; (3) r_e = 0.5, 0.6, 0.7, 0.8, 0.9, 1; (4) m = 0, 4, 12, 25, 100. The ranges of H/B and B/L used in this study are based on the previous work by Ukritchon et al. [22]. For the range of r_e, Krabbenhoft et al. [45] suggested that the value of this parameter should be between 0.5 and 1, which corresponds to the natural ratios of compressive and tensile undrained shear strengths. The range of m or ρB/s_uTC0 constitutes the combined effect of the excavation size B, the compressive shear strength at the ground surface s_uTC0, and the linear strength gradient ρ. In practice, s_uTC0 and ρ depend on the geological nature of the sites where the excavated width B can range from 1 to 20 m in practice. Theoretically, the ρB/s_uTC0 parameter ranges from 0 (homogeneous case) to a large value (non-homogeneous case). The homogeneous cases correspond to a case with ρ = 0 and/or a very large value of s_uTC0. For the non-homogeneous cases, they represent the cases with a relatively low s_uTC0 and/or a relatively large value of ρ.

3. Comparison for Model Validation

In the first step of the investigation, the stability numbers, N, determined by the rigorous FELA solutions, were compared with the published results in Ukritchon et al. [22]. The comparison shown in Figure 4a is for the effect of H/B on the stability number N, as well as its effect on various B/L with isotropic (r_e = 1) and homogeneous (m = 0) clays. Note that r_e is the anisotropic strength ratio, and m is the strength gradient ratio. Moreover, note that the present solution is the average (Avg) results calculated from the UB and LB FELA solutions. In general, the stability number increases with the increasing depth ratio H/B. The increase can be either nonlinearly or linearly, depending on the value of B/L. When B/L is smaller (B/L = 1/4, 1/8), fewer 3D constraints are observed, and a linear relationship between N and H/B is presented.

Whilst in Figure 4b, the comparison is made for (r_e = 1 and m = 4). It is interesting to note that, for the large strength gradient ratio such as m ≥ 4, N increases linearly with an increase in H/B for all values of B/L. Overall, the numerical results have shown an excellent agreement between the two solutions. The neglectable numerical differences between the two results can be attributed to the use of the perfectly plastic Tresca failure criterion in Ukritchon et al. [22] as opposed to the AUS failure criterion, with r_e = 1 used in the present study. To the best knowledge of the authors, there are currently no other values of r_e to be compared since this is the first work to consider the stability of unsupported rectangular excavations in anisotropic and non-homogeneous soils.

4. Results and Discussion

The effects of H/B on the stability number N are presented in Figure 5 for various values of r_e (the anisotropic strength ratio). Those shown in Figure 5a–f are for B/L = (0.25, 1.0) and m = (0, 12, 100). The numerical results have shown that the stability number N increases linearly with an increase in the excavation depth ratio H/B, except for the case of (B/L = 1.0 and m = 0). See Figure 5b for this special case of a square (B/L = 1.0) excavation in homogeneous (m = 0) clay, where N increases nonlinearly with the increasing H/B. One of the possible reasons could be attributed to the greater corner effects (geometrical arching). Note that the rate of increase in N (i.e., the gradient) increases as the strength gradient ratio m increases. Furthermore, note that a decrease in the anisotropic ratio r_e results in a decrease in the stability number. The selected failure mechanisms (shear dissipation) are presented in Figure 6 for the different values of H/B = (0.5, 1, 2, 3, 4). The comparison is based on the case of (r_e = 0.7, m = 4 and B/L = 1), and the results of the shear dissipation contour plots have shown a toe-failure mode for the shallow cases of H/B = (0.5 and 1). On the other note, for H/B > 1, a face-failure mode is obtained owing to the effect of the strength gradient ratio m.

Figure 7 shows the effects of B/L (the aspect ratio of the excavated site) on the stability number N for the various values of r_e (the anisotropic strength ratio). All of the values of m (m = 0, 4, 12, 25, 100) are considered for the chosen depth ratio H/B = 3, and they are presented in Figure 7a–e respectively. The numerical results have shown that N increases nonlinearly with the increasing B/L for all values of r_e. The gradient of the nonlinear curves becomes smaller as the strength gradient ratio m increases (see Figure 7a–d)—a linear relationship is observed for the case with m = 100. It is also noted that the stability number N decreases as the anisotropic strength ratio r_e decreases (transforming from isotropic to anisotropic soils). The comparison of five failure mechanisms for the various B/L = (1/8, 1/4, 1/2, 2/3, 1) is shown in Figure 8. The chosen plots are for H/B = 1 (r_e = 0.7, and m = 4). The shear dissipation contour plot of B/L = (1/8, 1/4) has shown a mechanism that resembles a 2D plane strain condition (see Figure 8a,b). As the value of B/L increases (so as the stability number N), a stronger system is presented, owing to full 3D corner effects (see Figure 8e for B/L = 1). Interestingly, a two-way failure mechanism is found in Figure 8e for B/L = 1. It should also be noted that the failure patterns are for the toe-failure mode in this shallow case of H/B = 1.

Figure 9 shows the relationship between the stability number N and the strength gradient ratio m for the various values of r_e (the anisotropic strength ratio). The presentations are for B/L = (1/8, 1) and H/B = (0.5, 1.0, 4.0). In general, an increase in m results in an increase in N. A linear relationship between N and m is observed in all investigated cases. Same as the previous discussions, the smaller the r_e, the smaller the stability number N. The chosen case for the failure mechanism comparison is presented in Figure 10 for (r_e = 0.7, H/B = 1, B/L = 1) with different values of m = (0, 4, 12, 25, 100). It should be noted that the size of the failure zone decreases as m increases. As a result, the failure mechanism changes from a toe-failure mode to a face-failure mode when m is larger than 4.

Figure 11 shows the relationships between the stability number N and the anisotropic strength ratio r_e for various values of m = (0, 4, 12, 25, 100). The plots are for the selected ratios of H/B = (0.5, 4) and B/L = (1/8, 1/2, 1). The numerical results have shown that the larger the m, the greater the stability number N. Overall, the stability number N varies linearly with the increase in the anisotropic ratio r_e. The rate of increase (gradient of the line) in N is dependent on the value of m. The larger the m, the greater the gradient of the line. Figure 12 shows a comparison of failure mechanisms among the various anisotropic ratios, r_e = (0.5–1). The comparison is for the excavation problem of (m = 4, H/B = 1, B/L = 1). The results have shown that the failure patterns are all in a toe-failure mode, and the variation of anisotropic ratio r_e does not seem to affect the failure size of the problem. The same conclusion can be drawn from Figure 13, where an additional study of m = 100 is presented. Indeed, as discussed previously, the face-failure mode is always the one observed for the large strength gradient ratio such as m = 100. It should be noted that all of the numerical results of this paper study are summarized in

Table 1. Stability numbers, N (r_e = 1.0 and 0.9).

m	re = 1						re = 0.9
	H/B	B/L					H/B	B/L
		1	2/3	1/2	1/4	1/8		1	2/3	1/2	1/4	1/8
0	0.5	4.559	4.372	4.234	3.959	3.860	0.5	4.299	4.128	4.027	3.764	3.460
	1	5.291	4.953	4.677	4.153	3.955	1	4.958	4.661	4.420	3.931	3.759
	2	6.420	5.969	5.553	4.637	4.125	2	5.968	5.560	5.201	4.384	3.917
	3	7.170	6.707	6.279	5.135	4.362	3	6.657	6.243	5.849	4.833	4.137
	4	7.770	7.280	6.862	5.606	4.632	4	7.216	6.774	6.388	5.262	4.378
4	0.5	9.259	9.008	8.758	8.362	8.255	0.5	8.749	8.495	8.288	7.870	7.898
	1	15.758	14.648	13.891	12.832	12.437	1	14.778	13.806	13.189	12.151	11.775
	2	30.562	28.105	26.095	22.663	21.168	2	28.283	26.251	24.639	21.455	19.996
	3	45.966	42.212	39.161	33.462	30.330	3	42.365	39.288	36.776	31.691	28.730
	4	61.106	56.160	52.214	44.704	39.942	4	56.528	52.402	48.986	42.238	37.914
12	0.5	18.213	17.662	17.175	16.735	16.561	0.5	17.156	16.730	16.444	15.785	15.516
	1	34.964	32.841	31.468	29.458	28.713	1	32.745	31.043	29.716	27.916	27.309
	2	69.931	65.641	62.571	57.045	53.981	2	65.470	61.913	59.186	53.975	51.219
	3	104.996	98.597	93.842	85.662	80.372	3	98.151	92.897	88.782	81.116	76.146
	4	140.088	131.336	125.132	114.290	107.172	4	130.918	124.024	118.326	108.284	101.830
25	0.5	32.254	31.504	30.613	30.072	29.207	0.5	30.536	29.713	29.120	28.351	28.198
	1	64.231	61.253	59.346	56.412	55.183	1	60.475	57.966	56.147	53.315	51.845
	2	128.460	122.685	118.597	111.345	107.223	2	120.981	115.931	112.269	104.690	101.445
	3	192.627	183.920	178.077	167.298	160.784	3	181.659	174.002	168.429	158.561	152.492
	4	256.720	245.260	237.102	223.242	214.300	4	241.932	231.924	224.580	211.538	203.318
100	0.5	111.392	109.320	107.744	106.143	105.742	0.5	105.461	103.419	101.948	100.684	100.869
	1	223.302	218.275	215.000	209.382	206.320	1	211.114	206.404	203.672	198.231	196.114
	2	446.620	436.122	429.972	415.484	411.352	2	423.437	413.457	407.257	396.090	389.648
	3	671.124	654.846	645.300	627.836	617.555	3	632.489	620.067	611.184	594.366	585.180
	4	894.322	873.216	859.562	837.478	823.428	4	845.674	827.874	813.880	793.196	780.282

,

Table 2. Stability numbers, N (r_e = 0.8 and 0.7).

m	re = 0.8						re = 0.7
	H/B	B/L					H/B	B/L
		1	2/3	1/2	1/4	1/8		1	2/3	1/2	1/4	1/8
0	0.5	3.998	3.880	3.807	3.531	3.255	0.5	3.696	3.571	3.503	3.273	3.043
	1	4.594	4.327	4.129	3.697	3.525	1	4.206	3.983	3.807	3.415	3.257
	2	5.502	5.126	4.825	4.103	3.668	2	5.006	4.694	4.418	3.790	3.398
	3	6.120	5.741	5.400	4.518	3.885	3	5.576	5.232	4.935	4.155	3.590
	4	6.608	7.654	5.868	4.896	4.100	4	6.010	5.658	5.350	4.494	3.788
4	0.5	8.192	8.020	7.753	7.399	7.252	0.5	7.574	7.352	7.212	6.878	6.741
	1	13.733	12.929	12.311	11.406	11.087	1	12.572	11.952	11.463	10.533	10.203
	2	26.032	24.326	22.942	20.133	18.797	2	23.642	22.283	21.126	18.623	17.413
	3	38.835	36.234	34.224	29.712	26.985	3	35.222	33.110	31.422	27.500	25.028
	4	51.836	48.330	45.558	39.574	35.560	4	46.896	44.168	41.810	36.482	32.908
12	0.5	16.080	15.921	15.174	14.900	14.827	0.5	14.856	14.759	14.108	13.673	13.539
	1	30.498	29.064	27.867	26.213	25.641	1	28.015	26.842	25.829	24.231	23.713
	2	60.747	57.799	55.401	50.671	47.946	2	55.688	53.260	51.218	46.939	44.337
	3	91.143	86.838	83.121	76.007	71.582	3	83.531	79.896	76.892	70.436	66.102
	4	121.484	115.712	110.738	101.452	95.462	4	111.290	106.608	102.410	93.984	88.160
25	0.5	28.631	27.955	27.377	26.578	25.636	0.5	26.442	25.683	25.257	24.657	24.445
	1	56.445	54.317	52.686	50.016	49.108	1	51.882	50.107	48.741	46.161	45.508
	2	112.831	108.652	105.268	99.169	95.067	2	103.890	100.505	97.455	91.487	87.891
	3	169.331	162.911	157.857	148.704	142.529	3	155.756	150.548	146.213	137.873	132.387
	4	225.288	231.924	210.434	198.498	190.700	4	207.658	200.882	194.796	183.884	176.652
100	0.5	98.804	97.060	95.754	94.274	92.989	0.5	91.379	89.479	88.734	87.369	82.820
	1	197.710	193.955	190.910	183.274	183.608	1	182.879	179.463	176.887	171.749	169.835
	2	395.618	387.909	382.134	372.104	364.822	2	365.556	359.051	353.148	340.475	337.916
	3	594.897	581.319	573.506	558.114	548.117	3	549.108	538.401	530.817	517.062	502.274
	4	794.012	775.738	763.674	744.168	731.800	4	733.122	719.498	706.896	689.530	677.690

and

Table 3. Stability numbers, N (r_e = 0.6 and 0.5).

m	re = 0.6						re = 0.5
	H/B	B/L					H/B	B/L
		1	2/3	1/2	1/4	1/8		1	2/3	1/2	1/4	1/8
0	0.5	3.317	3.256	3.181	2.955	2.771	0.5	2.944	2.882	2.820	2.628	2.319
	1	3.774	3.586	3.450	3.117	2.969	1	3.314	3.155	3.029	2.758	2.631
	2	4.479	4.207	3.977	3.444	3.097	2	3.901	3.674	3.495	3.050	2.743
	3	4.979	4.679	4.421	3.753	3.254	3	4.328	4.085	3.857	3.309	2.898
	4	5.366	5.048	4.784	4.054	3.440	4	4.666	4.390	4.188	3.564	3.052
4	0.5	6.823	6.675	6.550	6.294	6.027	0.5	6.038	5.960	5.753	5.559	5.491
	1	11.358	10.816	10.409	9.624	9.369	1	9.967	9.521	9.226	8.561	8.277
	2	21.138	19.958	19.048	16.942	15.852	2	18.563	17.571	16.751	15.037	14.072
	3	31.406	29.540	28.229	24.951	22.737	3	27.335	25.931	24.728	22.184	20.210
	4	41.728	39.450	37.670	33.366	29.942	4	36.506	34.454	32.832	29.342	26.568
12	0.5	13.475	13.208	12.811	12.563	11.985	0.5	11.915	11.702	11.428	11.133	10.896
	1	25.262	24.302	23.432	22.134	21.623	1	22.330	21.468	20.794	19.667	19.228
	2	50.092	48.210	46.566	42.644	40.585	2	44.187	42.254	41.138	37.943	35.950
	3	75.329	72.156	69.764	64.122	60.272	3	66.371	63.263	61.695	56.939	53.634
	4	100.316	96.282	92.980	85.524	80.114	4	88.410	84.236	82.204	75.902	71.372
25	0.5	24.129	23.467	23.007	22.585	22.364	0.5	21.192	20.744	20.493	19.967	19.757
	1	47.000	45.646	44.432	42.214	41.395	1	41.437	40.252	39.182	37.430	36.636
	2	93.925	91.067	88.624	83.521	80.038	2	82.907	80.487	78.697	74.134	70.683
	3	141.065	136.667	132.840	125.496	120.269	3	124.328	120.735	118.287	111.372	107.109
	4	188.000	182.288	177.256	167.254	160.868	4	165.788	160.782	157.210	148.552	141.700
100	0.5	83.089	81.711	79.615	79.360	75.411	0.5	73.433	72.459	71.454	70.765	69.833
	1	165.317	163.306	160.071	156.043	154.564	1	147.068	144.681	142.492	139.251	137.323
	2	332.289	326.784	322.154	313.455	305.386	2	294.210	289.922	285.932	278.151	272.439
	3	500.166	489.830	483.281	470.489	462.878	3	440.799	433.592	428.594	416.861	408.983
	4	664.460	653.052	644.744	627.270	617.136	4	587.558	579.212	571.564	556.762	545.928

.

5. Design Equations

A mathematical equation is developed and presented in this section by using a trial-and-error method of curve fitting. Nonlinear regression with multiple variables to the Avg bound solutions is employed to develop design equations for estimating the stability factor of unsupported rectangular excavations in clays with anisotropy and heterogeneity, as shown in Equation (7).

$$N=a_1+a_2\left(\frac{\rho B}{s_{uTC0}}\right)\left(\frac{H}{B}\right)+\frac{B}{L}\left[\left(a_3+a_4\frac{H}{B}\right)+\sqrt{\frac{\rho B}{s_{uTC0}}}\left(a_5\frac{H}{B}-a_6\right)+\left(\frac{\rho B}{s_{uTC0}}\right)\left(a_7\frac{H}{B}+a_8\right)\right]$$

(7)

where a₁ to a₈ are constant coefficients. To determine the optimal value of the constant coefficients (a₁–a₈), the nonlinear least square regression [54] is utilized. The sum of the squares of the deviation in N between the computed Avg solutions shown in Table 1, Table 2 and Table 3 and the approximate solutions from Equation (7) is then minimized to obtain the optimal values of constant coefficients.

Note that, to achieve high accuracy, Equation (7) is a “step-wised” equation developed for the different values of r_e. Using the complete data in Table 1, Table 2 and Table 3, the optimal values of the coefficients a₁ to a₈ for the different values of r_e are computed and presented in

Table 4. Constant coefficients for the proposed design equation.

Constant Coefficients	re
	0.5	0.6	0.7	0.8	0.9	1
a1	3.16697	3.3617	3.73044	4.1027	3.81764	4.28435
a2	1.35068	1.52276	1.66728	1.80244	1.92334	2.02793
a3	−0.69739	−0.55146	−0.64589	−1.47779	−0.09287	−0.55003
a4	0.46763	0.6014	0.72694	0.99956	1.01361	1.18849
a5	1.96752	2.24958	2.55588	3.08414	3.30362	3.39274
a6	0.92331	0.99065	1.14742	1.57205	1.5951	1.54619
a7	−0.07001	−0.07579	−0.081989	−0.12147	−0.13622	−0.13032
a8	−0.05667	−0.05068	−0.06182	−0.10530	−0.10229	−0.09884
R2	99.99%	99.99%	99.99%	99.99%	99.99%	99.99%

. On the other hand, the comparisons of N between the computed Avg bound solutions and the approximate solutions from Equation (6) are shown in Figure 14a–f, respectively, for different values of r_e = (0.5 to 1.0). It is pleasing to see the highly accurate solutions of the equation development—the coefficient of determination (R²) = 99.99%.

6. Conclusions

Rigorous stability solutions of the unsupported rectangular excavation in anisotropic and heterogeneous clays have been successfully studied in the paper using 3D LB and UB FELA. The stability number (N) that is a function of the excavation aspect ratio, B/L, the excavated depth ratio, H/B, the strength gradient ratio, m = ρB/s_uTC0, and the anisotropic strength ratio, r_e, was presented throughout the paper. The following conclusions are drawn based on the study.

• The stability number, N, increases with an increase in all of the investigated parameters of B/L, H/B, m, and r_e. The increases can be either in a linear or nonlinear relationship. The linear relationship was obtained for all investigated cases except for cases with smaller values of m, where a nonlinear relationship exists between N and B/L.
• The failure patterns of unsupported rectangular excavation in anisotropic and heterogeneous clays are either in a toe-failure mode (for small values of H/B, i.e., H/B = 0.5, 1) or a face-failure mode (for large values of H/B > 1) due to the effect of the strength gradient ratio m. For large values of m > 4, the failure modes are predominately the face-failure mode. The variation in the anisotropic ratio, r_e, does not seem to affect the failure size of the unsupported rectangular excavation problem.
• A new equation for predicting the stability number, N, of the unsupported rectangular excavation in anisotropic and heterogeneous clays is proposed. With the coefficient of determination (R²) being 99.99%, the proposed equation is highly accurate and useful for practical uses.

The proposed study provides deeper contextualized insights into the understanding of 3D unsupported excavations in undrained clay under the influence of soil anisotropy and heterogeneity. Future work directions may include the seismic stability performance as well as the soil random field probabilistic analysis.