*Article* **A Parametric Numerical Study for Diagnosing the Failure of Large Diameter Bored Piles Using Supervised Machine Learning Approach**

**Mohamed E. Al-Atroush 1,\*, Ashraf M. Hefny <sup>2</sup> and Tamer M. Sorour <sup>3</sup>**


**Abstract:** The full-scale static pile loading test is without question the most reliable methodology for estimating the ultimate capacity of large diameter bored piles (LDBP). However, in most cases, the obtained load-settlement curves from LDBP loading tests tend to increase without reaching the failure point or an asymptote. Loading an LDBP until reaching apparent failure is seldom practical because of the significant amount of settlement usually required for the full shaft and base mobilizations. With that in mind, the supervised learning algorithm requires a huge labeled data set to train the machine properly, which makes it ideal for sensitivity analysis, forecasting, and predictions, among other unsupervised algorithms. However, providing such a huge dataset of LDBP loaded to failure tests might be very complicated. In this paper, a novel practice has been proposed to establish a labeled dataset needed to train supervised machine learning algorithms on accurately predicting the ultimate capacity of an LDBP. A comprehensive numerical parametric study was carried out to investigate the effect of both pile geometrical and soil geotechnical parameters on both the ultimate capacity and settlement of an LDBP. This study was based on field measurements of loaded to failure LDBP tests. Results of the 29 applied models were compared with the calibrated model results, and the variation in LDBP behavior due to change in any of the hyperparameters was discussed. Accordingly, three primary characteristics were identified to diagnose the failure of LDBPs. Those characteristics were utilized to establish a decision tree of a supervised machine learning algorithm that can be used to predict the ultimate capacity of an LDBP.

**Keywords:** large diameter bored pile; hyperparameters; supervised machine learning; finite element method; parametric study; load transfer; failure mechanism

#### **1. Introduction**

The design of foundations has to satisfy two principal necessities [1]. First, complete failure of the foundation must be prevented with a sufficient margin of safety. Usually, the safety factor is assumed in practice to obtain the maximum safe foundation load (loadbased design approach). Second, the relative and total settlements of the foundation must be kept within limits that the superstructure can tolerate. Therefore, the settlement of the foundation under the working load has to be estimated to ascertain its effect on the superstructure (settlement-based design approach). Traditionally, both load and settlement control approaches are usually applied separately in the design process. However, Ref. [2] stated that the allowable load on the pile foundation should be obtained through a combined approach considering both soil resistance and its settlement inseparably acting together and each influencing the value of the other.

**Citation:** Al-Atroush, M.E.; Hefny, A.M.; Sorour, T.M. A Parametric Numerical Study for Diagnosing the Failure of Large Diameter Bored Piles Using Supervised Machine Learning Approach. *Processes* **2021**, *9*, 1411. https://doi.org/10.3390/pr9081411

Academic Editor: Li Li

Received: 15 July 2021 Accepted: 13 August 2021 Published: 16 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Significant contributions have been proposed to develop a settlement-based design approach for large diameter bored piles (LDBP) alongside the conventional capacity-based design approaches. Despite that, several studies [3–6] reported that most of the available methods for forecasting either the bearing capacity or the load-settlement behavior of LDBP invariably are associated with various degrees of uncertainty resulting from several factors. The mechanisms of pile foundations and pile–soil interaction are ambiguous, complex, and not yet entirely understood [7–10]. Additionally, several factors influence the pile settlement in a semi-infinite mass and a finite layer, as indicated by [11], i.e., the pile length to diameter ratio, Poisson's ratio, and others. Based on results of numerical analyses studies, it was argued by [12] that the skin friction and end bearing load-transfer responses of large diameter bored piles may invariably depend on the pile shaft length to diameter ratio (l/d), the relative strength and stiffness of the existing soils along the pile shaft and in the vicinity of the pile base, also the applied load level, and the accomplished amount of pile settlement. Those conclusions pinpoint that pile proportional and geometrical parameters (length and diameter) have a significant effect not only on pile settlement but also on the ultimate capacity of the pile, i.e., the portion of the load carried by end bearing resistance may be somewhat dominant in the cases of either an embedded short pile in a predominately homogeneous granular soil or an embedded long pile in a stronger/stiffer layer overlain by a weaker/soft layer. Incidentally, Ref. [12] proposed an l/d limit of "10" to distinguish between short and long piles.

Consequently, by necessity, many of the available methods have been mainly based on simplifications and assumptions [13]. This has led to limited success in terms of providing consistent and accurate predictions (Refs. [8,14–17]). Therefore, the in-situ pile loading test has been commonly accepted as the method providing accurate bearing capacity and settlement predictions. For design purposes, the full load–settlement response of the pile has to be well predicted and simulated; the designer can thus decide the ultimate load and comply with the serviceability requirement. Based on that, the loading test has been widely recommended by several international design standards (e.g., Refs. [18,19] and others).

Fellenius [20] defined the ultimate pile capacity at failure state, based on the pile load settlement performance as the load at which a considerable increase in the pile settlement occurs under a sustained or slight increase in the applied load. However, in most cases, load-settlement curves obtained from full-scale load tests conducted on large diameter bored piles tend to increase without reaching the failure point or an asymptote. Despite the reliability, loading of this class of piles (large diameter bored piles) until reaching apparent failure is seldom practical. An enormous amount of pile settlement is usually compulsory to achieve full friction and bearing mobilization [21–24]. Therefore, several extrapolation techniques, such as [25,26], among others, were developed to interpret the pile failure load (ultimate pile capacity) using the pile loading test data. Most of those methods are associated with several degrees of uncertainty as they were developed using the results of different cases of pile loading tests with a variety of pile geometries and geological conditions, and hence occasionally, any two give the same failure load (Refs. [5,15,27]). Despite that, these methods are widely accepted in many international design codes (e.g., Refs. [18,19]) to predict the ultimate capacity of large diameter bored piles if the in-situ loading test measurements do not indicate an apparent failure point.

Several international codes (e.g., Refs. [18,19,28]) suggest settlement-based failure methods for estimation of the ultimate capacity of large diameter bored piles in caseit is impossible to perform pile loading tests at the design phase. In addition, various correlations have been proposed in several pile design codes to predict the bearing capacity of large diameter bored piles based on the results of in-situ soil tests such as the Standard Penetration Test (SPT) and the Cone Penetration Test (CPT). However, most of the available forecasting methods either for the bearing capacity or the load-settlement behavior of large diameter bored piles in clay soils are deterministic in the sense that average/representative values of the soil properties incorporated along the affected zone of soil by the pile load are used. However, these pile design methods neglected the potentially existing variation of the in-situ soil properties along the pile shaft (inherent vertical soil variability) (Refs. [5,29–32]).

With that in mind, the load transfer method is generally a simple analytical procedure that can be applied in many complex situations, such as variation in the sections along a pile shaft and an inhomogeneous layered soil system. Therefore, there is still a need to develop a theoretically sound method that predicts a reliable value for the ultimate capacity of large diameter bored piles (LDBP). However, with the new technologies available recently, such as neural networks, artificial intelligence, and machine learning, the influence of both soil resistance and its settlement can be more accurately considered inseparably, acting together and each influencing the value of the other. Nevertheless, the machine should first learn the effect of each parameter influencing the pile behavior. Therefore, each influencing factor should be investigated individually to be used in the machine learning process [33–37]. The parameters whose values are used to control the learning process are called hyperparameters.

Supervised and unsupervised machine learning are the two basic approaches of artificial intelligence (AI) and machine learning. The main distinction between the two approaches is the use of labeled datasets [34,37]. Supervised learning uses labeled input and output data, while an unsupervised learning algorithm does not. These datasets are designed to train or supervise algorithms into classifying data or predicting outcomes accurately. Using labeled inputs and outputs, the model can measure its accuracy and learn over time. On the other hand, unsupervised learning uses machine learning algorithms to analyze and cluster unlabeled data sets. These algorithms discover hidden patterns in data without the need for human intervention. Thus, supervised learning models are ideal for sensitivity analysis, forecasting, and predictions, among others. Despite that, they require a huge data set to train the machine properly and obtain accurate results [34]. With that in mind, it is hard to have a huge data set for loading tests for the pile cases, especially for LDBP that rarely reach failure.

Certainly, the full-scale static pile loading test is the most reliable methodology for estimating the ultimate capacity of large diameter bored piles (LDBP). That is why several international geotechnical codes and foundation design standards (e.g., Refs. [18,19,28]) recommend this method to study and investigate the load transfer and failure mechanisms of this class of piles. However, in most cases, the obtained load-settlement curves from such tests conducted on LDBP tend to increase without reaching the failure point or an asymptote. Loading LDBP until they reach apparent failure is seldom practical because of the significant amount of pile settlement that is usually required to fully mobilize the pile shaft and reach the ultimate base resistances [21–27]. Huge test loads and hence highcapacity reaction systems should be used to accomplish the required enormous settlements. Thus, the targeted failure load may not always be practical to achieve, as reported in many case studies [21–27].

In this paper, a novel practice is proposed to establish a labeled dataset needed to train a supervised learning algorithm on accurately predicting the ultimate capacity of the LDBP. A comprehensive parametric study is carried out to investigate the effect of both pile geometrical and soil geotechnical parameters on both the ultimate capacity and settlement of LDBP installed in clayey soils. The procedure followed in this parametric study aims to explore the characteristic effect of each factor affecting the behavior of LDBP by comparing results of the numerical models with the measurements of a well-documented loaded to failure pile load test and assess the variation in both pile settlement and ultimate capacity due to change in either pile geometrical or soil geotechnical factors included in this study. Results of this parametric study are utilized to develop the decision tree needed to train supervised machine learning algorithms.

#### **2. The Reference Case Study**

A full-scale and well-instrumented load test was carried out by [38] at the location of Alzey bridge (Germany) to investigate the behavior of a large diameter bored pile (LDBP).

The instrumentation utilized in the LDBP loading test is described in Figure 1a. The length and diameter of the investigated pile were 9.50 m and 1.30 m, respectively. This LDBP was installed in overconsolidated stiff clay soil and subjected to axial loading cycles until achieving failure, as shown in Figure 1b. Test setup and the soil characteristics at the site location are also given in Figure 1b.

**Figure 1.** Large diameter bored pile loading to failure test (modified from [38]). (**a**) Measuring devices and instrumentation. (**b**) Test arrangement and typical soil profile with mechanical properties. (**c**) Field measurements of the loading test.

The main measurements of the well-instrumented load test are shown in Figure 1c. As shown, the significant increase in the measured settlement indicated apparent failure at the end of this test when the last load increment was applied on the pile head. More details of the loading test are available in [38].

#### **3. The Reference Numerical Calibration Study**

A numerical study has been carried out to simulate the response of the LDBP of the Alzey bridge case history [39]. Figure 2a shows the numerical model established to simulate the drained condition of overconsolidated (OC) stiff clay soil. The micro-fissures associated with the OC stiff clay are the main reason behind using the drained condition. These micro-fissures usually provide avenues for local drainage; soil along fissures has softened (increased water content) and is softer than intact material (more comprehensive discussion for using the drained condition with OC stiff clay is provided in [39]). Three constitutive soil models have been utilized to simulate the drained condition of the overconsolidated (OC) stiff clay soil. It was found that for this case history, the modified Mohr-Coulomb (MMC) constitutive model was superior to Mohr-Coulomb (MC) and the soft soil (SS) model in the simulation of the soil behavior [40,41]. The secant stiffness non-linear convergence method has been utilized to provide numerical stability required for the software solver to obtain convergence at substantial strain results (at failure). In addition, the sensitivity analyses performed also highlighted the significant effects of the mesh size and geometric dimensions on the analysis results. Results of this calibration study showed that excellent agreement was obtained between finite element results and the in-situ measurements of both the pile load settlement and load transfer relationships, as presented in Figure 2b,c.

**Figure 2.** Numerical model established to simulate the response of the LDBP of the Alzey bridge case history (after [39]). (**a**) Details of boundary conditions. (**b**) Deformed shape of the finite element mesh (MMC) under the failure load. (**c**) Comparison between field measurements and the numerical results (after [39]).

Based on the numerical calibration analysis performed, the finite element method was capable of predicting not only the working capacity but also the ultimate capacity of the large diameter bored piles (LDBP). Moreover, the large induced pile settlement at the failure state could be determined. The procedures that followed to calibrate the numerical models are available in [39].

#### **4. Methodology of the Parametric Study**

Factors affecting the response of large diameter bored piles in clayey soils are classified in this study into two main categories (Figure 3). The procedure followed in this parametric study aims to explore the characteristic effect of each geometrical or geotechnical factor affecting the behavior of LDBP. Therefore, a particular sequence will be followed in this study; each factor is explored separately. Measurements of the reference case are used to assess the variation in both pile settlement and ultimate capacity due to change in the factor under investigation, and the influences of the other parameters are filtered out at this step. In this way, it is possible to examine the specific effect of each factor.

**Figure 3.** Factors affecting the behavior of large diameter bored pile in cohesive soil.

The first category (A) includes pile geometrical parameters (pile diameter D and pile length L). The effect of these parameters on the response of large diameter bored piles is assessed by implementing different values for pile diameter and length in the finite element analyses. At the same time, soil properties are kept with the same values adopted in the reference calibrated model [39] without any change to reveal the variance in pile behavior due to the change in either pile diameter or length. Because the pile's diameter is the pivotal parameter in this study, nine numerical models for piles with equal lengths and different diameters are established to reveal the change in ultimate bearing resistance under the pile base during the failure. In contrast, only three finite element pile models with equal diameters and various lengths ranging from double to triple the original pile case's length are established to investigate skin friction changes.

On the other hand, the second category (B) involves the cohesive soil geotechnical parameters: effective cohesion c', effective friction angle Ø', lateral earth pressure coefficient K0, soil Young's modulus E, and dilatancy angle Ψ. The effect of each of these soil parameters on the response of LDBP in clayey soil will be individually evaluated by considering different values for each factor and assessing the change in the obtained ultimate pile capacity. Worth noting in this category, pile geometry parameters are kept with the same values adopted in the reference study [39]. Seventeen (17) new numerical analytical trials are performed for this purpose and compared with the calibrated model results.

#### *4.1. Adopted Parameters in the Numerical Models*

Modified Mohr-Coulomb (MMC) constitutive model has been used to define the drained condition of the homogenous overconsolidated stiff clay soil in all of the twentynine established numerical models for both categories A and B. The adopted parameters in each numerical model included in this study are summarized in Table 1. It is fundamental to note that the overconsolidation ratio (OCR) was set at a constant value of 2.0 for all cases [38], and its effect was not considered a variable in the parametric analyses. Additionally, the adopted values for shear strength parameters were selected to cover the different OC stiff clay soil classes in various drained conditions.


**Table 1.** The adopted geometrical and geotechnical parameters in the parametric study.

D: pile diameter. L: pile length. Eelastic: Young's modulus of Pile material. μ: pile Poisson's ratio. Ø': soil effective friction angle. c': soil effective cohesion. *<sup>v</sup>*: soil Poisson's ratio. K0: lateral earth pressure coefficient.t <sup>ψ</sup>: dilatancy angle. *<sup>E</sup>ref* <sup>50</sup> : secant stiffness in the standard drained triaxial test. *Eref oed*: tangential stiffness in oedometer test loading. *<sup>E</sup>ref ur* : elastic modulus at unloading. R: interface strength reduction factor.

#### *4.2. Numerical Modeling and Sensitivity*

The 29 numerical models of this study were established using MIDAS GTS NX finite element package. For category A, it was essential to perform sensitivity analysis to unify the effect of mesh size on the examined finite element models. This sensitivity study has been performed with the same procedure presented in the numerical reference study [39]. Results of these analytical trials showed that a fine mesh should be utilized to represent the pile element with a size that enables having at least two or three mesh elements in the vicinity of the pile base. This is important for bearing resistance results to eliminate the effect of mesh dependence at the level of the pile base [40]. Further, the boundary locations should be far away from the pile element for enough distance, not less than 15 times pile diameter (15D) in width, and at least 4 times the pile length (4L) in depth. These dimensions should be considered to ensure that the locations of boundaries do not affect the analysis zone.

Based on the results of sensitivity analyses performed, 12 two-dimensional axisymmetric models are established to investigate the effect of the pile geometrical parameters. Nine of them are prepared to simulate 10 separate single piles with equal lengths of 9.50 m and different diameters of; 0.40 m, 0.50 m, 0.60 m, 0.70 m, 0.80 m, 1.0 m, 1.20 m, 1.50 m, and 2.0 m. The other three numerical models are established to simulate three single piles with equal diameters of 1.30 m and different lengths of 13.0 m, 19.0 m, and 26.00 m. The 12 models' results were compared with the reference case study (1.3 m diameter and 9.50 m length) to assess the change in pile behavior due to change either in pile diameter or length.

In contrast, for category B, pile geometrical dimensions are taken with the same values of the reference case study, and the soil properties are changed according to the factor considered for evaluation in each section. The effect of the soil geotechnical parameters is investigated through 17 numerical models, and no change in the pile geometry is considered for this part of the study. Thus, the same mesh size and geometry dimensions adopted previously in the reference calibrated model have also been adopted to investigate the effect of category B parameters on pile response. According to [39], the positions of the model boundaries did not affect the analysis results of stresses and displacement around the pile when a numerical model height of 40 m and width of 20 m was adopted. Furthermore, the pile was represented by a very fine mesh with a mesh size of 0.216 m (quadratic 8-node elements). Additionally, with the same mesh size, a fine soil mesh media was established around and below the pile element. Gradually, soil mesh size was increased to be 1.0 m at the boundary locations.

According to the calibration study [39], very good agreements were obtained between field measurements of the Alzey bridge case study and the results of the numerical analyses, even at the failure state, when an interface strength reduction factor (R) of 1.0 was adopted in the analyses. Consequently, the same value of the reduction factor was utilized in the 29 parametric models.

#### *4.3. Stages of Analysis*

Similarly to the analytical procedure utilized in the numerical calibration study, the analysis sequence followed in this parametric study consists of three stages of analysis. The first stage represents the initial stresses of the soil due to the overburden effect. This stage is essential to simulate the initial condition of the soil before installing the pile. The second stage started by changing the pile volume to concrete material instead of soil material to represent the pile installation phase. It is worth noting that a rigid interface element has been used in the second stage to provide the required numerical stability at this stage, as it is rigidly connecting the paired nodes of soil and pile to avoid singularity. In the third stage, the rigid interface has been replaced with the frictional interface, and the uniform load is applied to the pile head using incremental load steps (250 kN per load increment). The applied load in each case is defined with a value greater than the estimated value of ultimate pile capacity (using [42]) for each of the 30 pile cases to allow the numerical solver to achieve the failure according to the adopted convergence criteria.

#### *4.4. Load Transfer Mechanism and Failure Criteria*

The maximum load where convergence can be achieved in the numerical analysis is considered as the failure load. Failure is also indicated by the apparent large settlement that is expected to be induced under the application of the ultimate load. The pile load transfer mechanism is investigated by determining the pile stress in the vicinity of the pile base and multiplying it by the cross-sectional area of each pile case to obtain the pile bearing resistance at each load increment. Thus, pile friction resistance can be calculated by deducting the obtained pile bearing load from the total applied load at each loading increment. Then, the relations between pile load settlement, pile friction, and bearing capacities at each loading increment can be plotted. In contrast, the formation sequence, shape, and size of the formed plastic bulb at the failure status are explored at the obtained ultimate load of each pile model.

#### **5. Analysis Results and Discussion**

#### *5.1. Large Diameter Bored Pile Load-Settlement Relationship*

For pile geometrical parameters (category A), load-settlement results of the nine numerical pile models with different diameters and equal length (9.50 m) are presented in Figure 4a. Additionally, the results of the three numerical equal-diameter (1.30 m) pile models with different lengths are shown in Figure 4b. Those results were compared with the load–settlement relationship of the calibrated model (Figure 1) to assess the change in pile load settlement performance due to variation in either pile diameter D or pile length L.

Figure 4a,b highlight the effect of pile diameter and length on the ultimate pile capacity and the induced settlement at failure. Hence, the ultimate capacity increases with increases in pile diameter and length. Consequently, the induced settlement at the failure state is also increased. A significant increase in pile settlement is also observed at the last load increment in the results of the 12 finite element models (category A), as the induced pile settlement at the last load increment (failure load) ranged from about one and a half to more than two times the obtained settlement value at the pre-last loading increment.

Similarly, Figure 4c–e highlights that ultimate pile capacity is increased, and consequently, the induced settlement at failure is also increased with increases of the following soil geotechnical parameters from category B: clay effective cohesion (c'), effective friction angle (Ø'), and soil dilatancy angle (Ψ). Conversely, the induced pile settlement at failure is decreased with increases in the lateral earth pressure coefficient (K0), despite the observed increase in ultimate pile capacity, as shown in Figure 4f. This is attributed to the confinement of the surrounding soil due to stress increases according to lateral earth pressure coefficient increases. Conversely, the increase in soil Young's modulus (E) also caused a decrease in the induced pile settlement at failure, but without any observed change in ultimate pile capacity, as shown in Figure 4g.

#### *5.2. Pile Load Transfer Mechanism*

Based on the methodology explained before in Section 4, the ultimate total, bearing, and friction capacities are determined for each of the 29 numerical models performed in this study. The effect of each factor, either from category A or B, on the ultimate pile capacity is given in Figure 5a–g. Figure 5a–g highlights that the three ultimate capacities (total, bearing, and friction) are increasing with increases in pile diameter (D), length (L), effective clay cohesion (c'), effective friction angle (Ø'), and soil dilatancy angle (Ψ). However, Figure 5e reveals that the lateral earth pressure coefficient significantly affects the pile friction resistance; however, it has almost no effect on the ultimate bearing capacity. In contrast, Figure 5f indicates that soil Young's modulus has almost no effect on the ultimate friction, bearing, and total capacities of the large diameter bored pile. However, it has a major impact on the settlement of large diameter bored pile, especially at the failure state (see Figure 4g).

**Figure 4.** Effect of the pile geometrical and soil geotechnical parameters on the load settlement relationship. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Effective cohesion (c'). (**d**) Effective friction angle (Ø'). (**e**) Dilatancy angle (Ψ). (**f**) Lateral earth pressure coefficient (k0). (**g**) Young's modulus (Es).

**Figure 5.** Variation of the ultimate bearing (Qb), friction (Qf), and total capacities (Qt) with the different parameters. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Soil effective cohesion (c'). (**d**) Effective friction angle (Ø'). (**e**) Lateral earth pressure coefficient (K0). (**f**) Young's modulus (Es). (**g**) Soil dilatancy (Ψ).

In the same line, percentages of load transferred by bearing and by friction are calculated relative to total obtained ultimate capacity (QB/QT and QF/QT) for each of the 29 numerical models, and the variations of those percentages with each factor in both categories are given in Figure 6a–f.

**Figure 6.** Variation of (QB/QT) and (QF/QT) with different parameters. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Effective cohesion (c'). (**d**) Effective friction angle (Ø'). (**e**) Lateral earth pressure coefficient (k0) (**f**) Dilatancy angle (Ψ).

It can be seen from Figure 6a that the percentage of load transferred by bearing (QB/QT) increases from 10% to about 50%, with pile diameter increases from 0.40 m to 2.0 m. Conversely, the percentage of load transferred by friction decreased from about 90% for a pile with a diameter of 0.40 m to 47.3% for a pile with a diameter of 2.0 m.

It was also observed that the applied load was predominantly transferred by friction for the first nine cases (from 0.4 to 1.50 m in diameter). Especially in small diameter piles, the percentage of the load transferred by friction ranged from 90% (0.40 m small diameter pile) to 81% (0.60 m small diameter pile). An observed decrease in the percentage of the

load transferred by friction was noted after the pile diameter increase from 0.5 m to 0.6 m, as the Qf/QT decreased by about 10%.

For large-diameter piles (with diameters greater than 0.60 m), the load transferred by friction linearly decreased from 79% (a pile diameter of 0.70 m) to 58% (at a pile diameter of 1.50 m) with pile diameter increases. Thus, the total load friction and bearing percentages became nearly equal (52.6%:47.3%) in the last case (pile with 2.0 m diameter). These percentages indicate an apparent difference between the behavior of large and small diameter bored piles.

Figure 6b demonstrates that the applied load was predominantly transferred by friction (Qf) for the four cases with different lengths. The percentage of load transferred by bearing decreases from about 40% to about 20% with pile length increases. Consequently, the percentage Qf/QT is increased from about 60% (pile with length 9.50 m) to about 80% in the last case with the pile of 26 m length.

In summary, the percentages presented in Figure 6a,b indicate that the rate of load transferred by friction increases with increases in pile length. In contrast, this percentage (load transferred by friction) decreases with increases in pile diameter.

The effects from the category B parameters of effective clay cohesion (c'), effective friction angle (Ø'), lateral earth pressure coefficient (K), soil Young's modulus (E), and soil dilatancy angle (Ψ) on the percentage of load transferred by bearing or friction are given in Figure 6c–f.

It can be seen from Figure 6c that the percentage of load transferred by bearing (48%) nearly equals the percentage of load transferred by friction (52%) for a small value of soil cohesion of 5 kN/m2. The percentage QF/QT increased from about 52% (5 kN/m2) to near 63% for the case with an effective cohesion of 30 kN/m2. For the last three cases with effective cohesion higher than 20 kN/m2, the ratio between the load transferred by bearing and by friction of the total ultimate load tended to be constant (37% and 63%, respectively). Those percentages indicate that effective soil cohesion has the same impact on both bearing and friction pile resistances at the ultimate state.

Figure 6d demonstrates that the percentage of load transferred by friction decreases from about 70% (Ø' = 10◦) to about 53% in the fourth case with an effective friction angle of 30◦. In addition, the percentage QB/QT (47%) nearly equals the percentage of load transferred by friction (53%) in the last case with a friction angle of 30◦. An observed increase is noticed in the percentage of load transferred by bearing with effective friction angle increases, as it increases from about 30% at an effective friction angle of 10◦ to about 47% at an effective friction angle of 30◦. These percentages indicate that effective soil friction angle has a significant effect on the percentage of the load transferred by bearing, which is expected because of the related increase in soil bearing capacity with the increase in effective friction angle. Of note, in these trial analyses, the same value for the lateral earth pressure coefficient (K0) was adopted for all four of the cases with different effective friction angles.

It can be seen from Figure 6e that the percentage of load transferred by friction increased from about 52% at ko of 0.43 to about 62% in the third case with lateral earth pressure coefficient of 0.8; also, the percentage of load transferred by bearing (48%) nearly equals the percentage of load transferred by friction (52%) in the first case with ko of 0.43. These percentages indicate that the soil lateral earth pressure coefficient has a significant effect on the percentage of the load transferred by friction, which is attributed to the related increase in soil shear strength with increases in k0. On the other hand, Figure 6f highlights that the percentage of load transferred by friction and by bearing is almost unaffected by the change in dilatancy angle value, as about 63% of the applied load is carried by friction, and about 37% is transferred with bearing, in the four cases with different dilatancy angles.

In the next two sections, the results of the numerical models of this parametric study will be used to assess the calculated ultimate capacity of the LDBP using two different methods of both capacity-based and settlement-based design approaches. According to the available data in a field study [38], the Meyerhof 1976 capacity-based method [42] and Egyptian code settlement-based design (ECP 202/4) [19] approaches were chosen for this evaluation because other methods, such as AASHTO LRFD [28], require soil testing results such as SPT or CPT test results that were not provided in the field study.

#### *5.3. Average Ultimate Bearing Stress*

Results of the numerical models are used to explore the effect of several geometrical and geotechnical parameters on the bearing stress below the pile base at the failure state. Pile diameter (D), length (L), clay effective cohesion (c'), and soil effective friction angle (Ø') are the selected parameters to represent the two categories (A and B) in this evaluation. Numerical results are compared with estimated values of bearing stress using the chosen settlement-based design method of [19], as given in Figure 7a–d.

**Figure 7.** Influence in bearing stress with various parameters. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Effective cohesion (c'). (**d**) Effective friction angle (Ø').

As shown in Figure 7a the almost equal value of bearing stress (about 600 kN/m2) can be seen for piles with small diameters of 0.40 m and 0.50 m. A significant increase in bearing stress value was observed in the case of a pile with a diameter of 0.60 m, as the bearing stress increases to about 800 kN/m2, which is about 33% greater than the obtained bearing stress at 0.4 and 0.5 m pile diameters. Further, almost the same value of bearing stress was obtained below the bases of piles with diameters of 0.6, 0.7, 0.8, and 1.0 m. Bearing stress increases again with pile diameter increases for piles with diameters greater than 1.0 m, to achieve its maximum value of 1050 kN/m<sup>2</sup> below the 2.0 m diameter pile base.

It can be seen from Figure 7b that the bearing stress increases with increases in pile length, which is expected at the drained condition because of the increase in vertical soil stress at the pile base level. The bearing stress increases from 920 kN/m2 below the base of the pile with length 9.50 m, to nearly two times this value with a length of 26.0 m (1985 kN/m2).

On the other hand, Figure 7c presents the obtained bearing stress values below the pile bases of the six finite element models with different soil effective cohesion values (category B). It can be seen from this Figure that the bearing stress increases with increases in soil effective cohesion. The bearing stress increases from about 800 kN/m2 in the first case with an effective cohesion of 5 kN/m2, to about three times that value (2113 kN/m2) for the case with effective cohesion of 100 kN/m2. Similarly, Figure 7d shows an increase in bearing stress due to the increase in soil effective friction angle.

Although the results of the numerical models showed that the ultimate bearing stress is influenced by the pile diameter, pile length, soil effective cohesion, and soil effective friction angle, the bearing stress estimated by the settlement-based method (ECP 202/4) is independent of pile diameter, pile length, and soil effective cohesion. Most of the available settlement-based methods depend on settlement-based criteria to estimate the ultimate bearing resistance of the pile at a particular settlement value. Those methods generally defined the nominal ultimate pile capacity (Qtu) using different percentages of the pile head settlement to diameter ratio (s/d), such as 10%, according to ECP 202/4 [19].

#### *5.4. Average Ultimate Unit Skin Friction*

The effects of pile diameter (D), length (L), clay effective cohesion (c'), and soil effective friction angle (Ø') on the soil unit skin friction are investigated using the results of the numerical models. The obtained numerical results are compared with the calculated values of the interface shear strength (according to Mohr-Coulomb failure criteria (σ<sup>h</sup> tan Øi + ci)), to ensure that full mobilization has been achieved in each case.

As shown in Figure 8a–d, all numerical models reached the full friction mobilization. With that in mind, numerical results of the soil unit skin friction are also utilized to assess the ones determined using the chosen capacity-based method [42].

**Figure 8.** Variation of the average unit skin friction with various parameters. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Soil effective cohesion (c'). (**d**) Soil effective friction angle (Ø').

As demonstrated in Figure 8a, the obtained results for the average unit skin friction of the 10 pile cases with different diameters are almost equal. This indicates that the unit skin friction is not affected by the pile diameter (D). In contrast, Figure 8b highlights that the unit skin friction of the four pile cases with different lengths increases with increases in pile length. This is attributed to the increase in vertical stress due to the increase in overburden pressure at the pile base level. Additionally, Figure 8c,d highlights that the obtained average unit skin friction results increase with the increase in effective soil cohesion and effective friction angle.

On the other hand, the calculated soil skin friction values using the capacity-based method [42] are close to the numerical results, except for the three cases with different lengths (Figure 8b) and the last two cases with greater friction angles (Figure 8d). The difference between the numerical results and those obtained using the Meyerhof formula ranges from 10% to 21%. However, the Meyerhof method underestimated the unit skin friction of all cases. Of note, the capacity-based method takes the effect of pile diameter (D), length (L), clay effective cohesion (c'), and soil effective friction angle (Ø') into account, in contrast to the settlement-based methods that ignore the effect of several influencing factors, as presented in Section 4.3.

#### *5.5. Size of Plastic Bulb below Pile Base*

Figure 9a–d presents the shape of the formed plastic bulbs obtained from numerical analyses at the last load increment (failure load) for different pile model cases with various geometrical and geotechnical parameters. Plastic points are generally formed when the transferred stress achieves a value that almost equals the soil shear strength (according to Mohr-Coulomb's theory of rupture). Hence, the yielding occurs, and plastic strains are induced. The zone at which the plastic points are concentrated under the pile base is hereafter referred to as the plastic bulb. These plastic bulbs have been measured in diameter and length to assess the effect of each parameter considered in this study on the failure mechanism of large diameter bored piles.

Because the pile diameter is the main governing factor in this comparison, it was essential to investigate the differences between large and small diameter piles in terms of the size of the end bearing plastic bulb that formed under the pile base in the failure state. Figure 9a presents the change in width and length in the formed plastic bulb with pile diameter change. As shown, the plastic bulb is obviously increased in size and height with increases in pile diameter. The diameter of the plastic bulb (Dp) under the base of the first pile case with a diameter of 0.40 m is 2.14 m (~5D), and the greatest size of the plastic bulb is at 7.08 m (3.54D) when pile diameter increases to 2.00 m. The length of the plastic bulb below the pile bases (Lp) also increased with pile diameter increases, to about 5.14 m (~2.5D) below the 2.00 m diameter pile base. Of note, in the 10 cases with different diameters, the length of the plastic bulb below the pile base is nearly two and a half times the pile diameter.

In contrast, as shown in Figure 9b, plastic bulbs of almost identical diameter are obtained in all cases of piles with equal diameters and different lengths. A minor increase is observed with pile length increases. For the pile with a length of 9.50 m, the plastic bulb diameter measured at 3.68 times pile diameter and slightly increased to 4.05 D for the fourth case of the pile with a length of 26.0 m. Furthermore, the length of the plastic bulb slightly increased below the pile base from three times pile diameter at the pile with a length of 9.50 m to almost 3.10 of pile diameter at the pile with a length of 26.0 m.

**Figure 9.** Influence in the formed plastic bulbs due to changes in different parameters. (**a**) Pile diameter (D). (**b**) Pile length (L). (**c**) Soil effective cohesion (c'). (**d**) Soil effective friction angle (Ø').

Figure 9c indicates that the size of the formed plastic bulb is affected by effective soil cohesion. The formed plastic bulb increases in width and length as the soil effective cohesion increases, i.e., the plastic blub radius (rp) is measured as 2.66 m (~4D) in the first case (c' = 10 kN/m2) and increases to 3.6 m (5.54D) for the last case (c' = 100 kN/m2). The height of the plastic bulb below the pile base ranges from 3.11 to 3.62 times the pile diameter. Similarly, as presented in Figure 9d, the plastic bulb also increased in size and length with soil effective friction angle increases.

#### **6. Observations of the Ultimate Capacity of the Large Diameter Bored Piles**

Based on the results of the numerical models investigated within this parametric study, it was noted that full friction mobilization occurs when the transferred shear stresses from the shaft interface achieve a value that almost equals the soil shear strength (according to Mohr-Coulomb theory of rupture). This may be described as the beginning of the failure stage, as in this phase, the pile friction resistance tends to be constant or slightly decreases (See Figure 1). The full friction mobilization was also indicated in the numerical models by the formed plastic points that extended above the pile base along an interface length of more than three times the pile diameter (3D). Additionally, in many cases, the plastic points extended above the base to fully cover the whole length of the pile shaft interface.

After full mobilization, the pile friction resistance tends to be constant (inactive constant friction situation). However, the pile is still able to sustain higher applied loads safely through its bearing resistance. The applied load is predominantly transferred by bearing within the failure stage. Consequently, the transferred high bearing stresses at the pile base level resulted in a significant increase in both the size and length of the formed plastic bulb below the pile base. Finally, the transferred bearing stress value will achieve the value of the soil's ultimate soil bearing capacity. In this situation, the pile will transition from the inactive constant friction situation into the sliding friction situation, and apparent failure will be observed through the large induced pile settlement at the end of this stage. Hence, the applied load at which the failure is achieved could be described as the ultimate load of the large diameter bored pile.

Results of the numerical model simulation indicate that the induced settlement at the failure state (Sf) ranges from 1.5 to more than 2 times the obtained settlement (Sf-1) at the pre-last load increment (90% of the ultimate load), as shown in the results summary presented in Figure 10.

**Figure 10.** The ratio between the induced settlement at the failure load (Sf) and the settlement (Sf-1) under the load before the last load increment (0.9Qult).

To investigate the size and length of the end bearing plastic bulb at the failure stage, the ratio (Dp/D) between the plastic bulb diameter (Dp) and pile diameter (D) is calculated for (19) sample models with different pile proportions and soil properties as shown in Figure 11. In addition, the ratios (Lp/D) between the length of plastic bulb below the base level (Lp) to the pile diameter are obtained for the 19 models and shown in the same Figure. Results of the numerical models indicate that the plastic bulb diameter (Dp) ranges from 3 to 6 times the pile diameter (D). Moreover, the plastic bulb length (Lp) ranges from 2 to 4 times the pile diameter (D). These percentages are used to diagnose the expected plastic zone around the base of LDBP, as shown in Figure 12. Significant plastic deformations are highly expected to be induced in this area, which may cause the arching phenomenon that leads to the failure state and sudden excessive settlement of LDBP.

**Figure 11.** Relation between pile diameter and length and diameter of the formed plastic bulb below the pile base at the failure state.

**Figure 12.** Schematic diagram showing the formed plastic bulb's predicted size and the induced settlement at the failure state of large diameter bored piles.

#### **7. Application of the Study**

Based on the results of the numerical parametric study performed, the failure of LDBP can be diagnosed by the ratio between the induced settlement at the failure load (Sf) and the settlement (Sf-1) at the pre-last load increment (90% of the ultimate load). Moreover, it can be identified through the size and height of the plastic bulb formed around and below the base of LDBP at the failure state, as explained in the previous section. In addition, the achievement of full friction mobilization can be ensured using Mohr-Coulomb failure criteria (σ<sup>h</sup> tan Øi + ci). Thus, those three characteristics can be used to formulate a supervised machine learning algorithm. With that in mind, the parametric study results can extend the measurements of only one loaded to failure LDBP (Alzey case history) to provide the supervised algorithm with a data set for 29 cases with various pile geometrical and soil geotechnical parameters.

Figure 13 presents the decision tree that can be used to train the machine learning algorithm using the field measurements of the loading test and the parametric study results. The primary input in the tree will be the site investigation data needed to estimate the soil parameters, in addition to the pile geometry and field measurements of the available loading test. In the first layer of analysis, the algorithm should estimate both the ultimate pile capacity and settlement behavior. Secondly (layer 2), the estimated results shall be compared with the field measurements (layer 3) to assess the accuracy of the calculations. If the accuracy of the results exceeds the 95% benchmark (layer 4), the algorithm can open the data set of the parametric data set (Layer 5) and start performing the same procedure (forward process). In case the accuracy is below the benchmark, the algorithm shall re-analyze the data until hitting the benchmark (backward process).

**Figure 13.** The suggested decision tree of a supervised machine learning algorithm for predicting the ultimate capacity of LDBP.

Consequently, the algorithm can be trained using the loading test measurements and the supervised data set of the parametric models. Thus, the algorithm will be ready to start the prediction process for different cases after achieving the targeted accuracy benchmark for all cases provided in the data set.

#### **8. Conclusions**

In this study, 29 numerical models were established to investigate the effect of different geometrical and geotechnical parameters on the response of large diameter bored piles in clayey soils by comparing their results with the results of the calibrated model of the Alzey Bridge case study. The following conclusions are drawn:


**Author Contributions:** Conceptualization, M.E.A.-A. and A.M.H.; methodology, M.E.A.-A. and A.M.H.; experimental testing, validation M.E.A.-A., A.M.H. and T.M.S.; formal analysis, M.E.A.-A. and A.M.H.; investigation, M.E.A.-A., A.M.H. and T.M.S.; data curation, M.E.A.-A. and T.M.S.; writing—original draft preparation, M.E.A.-A.; writing—review and editing, M.E.A.-A., A.M.H. and T.M.S.; visualization, M.E.A.-A.; supervision A.M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The research is supported by the Structures and Materials (S&M) Research Lab of Prince Sultan University, Riyadh, Saudi Arabia.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**


## *Article* **Precomputation of Critical State Soil Plastic Models**

**Vicente Navarro \*, Virginia Cabrera, Gema De la Morena, Daniel González, Laura Asensio and Ángel Yustres**

Geoenvironmental Group, University of Castilla-La Mancha, Avda. Camilo José Cela s/n, 13071 Ciudad Real, Spain; virginia.cabrera@uclm.es (V.C.); gema.delaMorena@uclm.es (G.D.l.M.); Daniel.Gonzalez@uclm.es (D.G.); Laura.Asensio@uclm.es (L.A.); Angel.Yustres@uclm.es (Á.Y.)

**\*** Correspondence: vicente.navarro@uclm.es; Tel.: +34-926295300 (ext. 3264)

**Abstract:** In this paper, a simple precomputing procedure is proposed to improve the numerical performance of the technological application of critical state soil models. In these models, if associated plasticity is assumed, the normalization of the stress space allows both the yield surface and the plastic components of the elastoplastic matrix to be defined as a function of a single variable. This approach facilitates their parameterization and precomputation, preventing the repetition of calculations when the boundary value problems appear at the yield surface with the calculation of plastic strain. To illustrate the scope of the procedure, its application on a modified Cam Clay model is analysed, which shows that the method allows a significant reduction of about 50% (as compared with the conventional explicit integration algorithm) in the computational time without reducing the precision. Although it is intended for critical state models in soils, the approach can be applied to other materials and types of constitutive models provided that parameterization is possible. It is therefore a methodology of practical interest, especially when a large volume of calculations is required, for example when studying large-scale engineering systems, performing sensitivity analysis, or solving optimization problems.

**Keywords:** plasticity; critical state soil model; modified Cam Clay model; model normalization; precomputation

#### **1. Introduction**

Simulation of soil yield processes is often computationally expensive, especially when analysing active clays such as MX-80 bentonite [1], FEBEX bentonite [2], or GMZ bentonite [3]. These materials, which usually have a high smectite content, when hydrated unconfined from partially saturated conditions, can undergo deformations of 500%. Under confined conditions, they can lead to swelling pressures approaching 10 MPa. This makes the simulation of their hydromechanical behaviour always computationally demanding, particularly when simulating complex domains or when simulations have to be repeated a large number of times. Even considering the development experienced in recent years by microelectronic technology, which has allowed the multi-core and many-core hybrid heterogeneous parallel computing platform to facilitate a very important advance in computing power, the efficiency of the calculation algorithm continues to be a key issue in the application of massive calculation processes. This is shown, for example, in [4] where the application of the Generalized Likelihood Uncertainty Estimation method in the probabilistic estimation of parameters in hydrology is analysed. In the field of soil mechanics, the improvement of computational performance is especially important in design or parameter identification processes, where the computational time can determine the viability of the study [5].

In some constitutive models, such as the critical state soil models, both the stress space [6] and the plastic terms **D**<sup>p</sup> of the elastoplastic matrix **D**ep [7–9] can be normalized. Thus, different stress states are associated with the same point of the normalized yield surface and with the same **D**p. The result could be that when solving boundary value

**Citation:** Navarro, V.; Cabrera, V.; De la Morena, G.; González, D.; Asensio, L.; Yustres, Á. Precomputation of Critical State Soil Plastic Models. *Processes* **2021**, *9*, 2142. https:// doi.org/10.3390/pr9122142

Academic Editor: Avelino Núñez-Delgado

Received: 10 November 2021 Accepted: 25 November 2021 Published: 27 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

problems, for a nonnegligible number of times, the calculations associated with the plastic behaviour of the same point (or in the vicinity of the same point) of the normalized stress space are repeated. It is therefore interesting to precompute **D**<sup>p</sup> for a grid of normalised stresses and to interpolate any value within the points of that grid when solving boundary value problems.

This approach is the resolution strategy proposed in the present work. In the following sections, the normalisation of the yield surface and **D**<sup>p</sup> is first described. Then, the CPU time savings that the precomputation entails are evaluated. Later, the criteria are defined to ensure that the proposed method does not compromise the quality of the simulation, and finally, an inspection exercise is performed to evaluate its performance in solving boundary problems.

#### **2. Theoretical Background**

In addition to saturated conditions, this paper also assumes axisymmetric conditions. These hypotheses do not reduce the scope of the proposed formulation and yet allow it to be explained more clearly. Furthermore, in line with this quest for clarity, the modified Cam Clay model (MCCm) [6] is adopted as the reference critical state soil model. In this way, the complexity of the constitutive formulation is reduced to a minimum, and the present work is focused on the description of the proposed calculation method. However, the widespread use of the MCCm in computational geotechnics [10–12] means that the methodology presented here is not trivial.

Usually, the yield surface *f* of the MCCm is formulated as

$$f \equiv \left(\frac{q}{M}\right)^2 - p(p\_\bullet - p) = 0,\tag{1}$$

where *q* is the von Mises deviatoric stress, *M* is the slope of the critical state line, and the isotropic preconsolidation pressure *p*<sup>O</sup> is the hardening parameter. This parameter is a good normalisation parameter since by introducing the dimensionless stresses *P* = *p*/*p*<sup>O</sup> and *Q* = *q*/*p*O, *f* projects onto *F* (Figure 1), as defined by the expression

$$F \equiv \left(\frac{Q}{M}\right)^2 - P(1 - P) = 0.\tag{2}$$

**Figure 1.** Projection of the different yield surfaces *f* <sup>i</sup> on the normalized yield surface *F*.

All the points (*p*,*q*)i located on the same axis *η*, although they are associated with a different hardening parameters *p*Oi, define the same dimensionless stress (*P*,*Q*) (Figure 1). Moreover, in all of them, the vector **n**, normal to *f*, given by the expression

$$\mathbf{n} = \left(\frac{\partial f}{\partial p'} \frac{\partial f}{\partial q}\right) = \left(2p - p\_{\rm O} \frac{2q}{M^2}\right) = p\_{\rm O} \left(2P - 1 \prime \frac{2Q}{M^2}\right) = p\_{\rm O} \mathbf{N} \tag{3}$$

has the same direction of **N**, with its modulus varying according to the *p*<sup>O</sup> factor. Therefore, assuming an associated plasticity model, the plastic terms of the elastoplastic matrix will be the same for all yield stresses which project onto the same dimensionless stress. Indeed, if the elastoplastic matrix **D**ep defines the increase in the effective stress *d*σ associated with an increase in the strain *d*ε (the Voigt notation is used for both stress and strain tensors), then we have

$$d\sigma = \mathbf{D}^{\mathrm{ep}} \cdot d\varepsilon \tag{4}$$

If the consistency equation is applied to the definition of the yield surface *f* and both the flow rule and the hardening law are considered, then **D**ep can be expressed as (see, for example, [13,14])

$$\mathbf{D}^{\rm cp} = \mathbf{D}^{\rm c} \cdot \left( \mathbf{I} - \frac{\left(\frac{\partial f}{\partial \sigma}\right)^{\rm t} \cdot \frac{\partial f}{\partial \sigma} \cdot \mathbf{D}^{\rm c}}{\frac{\partial f}{\partial \sigma} \cdot \mathbf{D}^{\rm c} \cdot \left(\frac{\partial f}{\partial \sigma}\right)^{\rm t} - \frac{\partial f}{\partial \mathbf{x}} \cdot \frac{\partial \mathbf{x}}{\partial \varepsilon \mathbf{F}} \cdot \left(\frac{\partial f}{\partial \sigma}\right)^{\rm t}} \right) = \mathbf{D}^{\rm c} \cdot (\mathbf{I} - \mathbf{D}^{\rm p}) \tag{5}$$

where **D**<sup>e</sup> is the elastic strain matrix, **I** is the identity matrix with the dimension of **D**e, the symbol (·) <sup>t</sup> indicates the transpose operator, and **x** defines the vector of hardening parameters. As noted, in the MCCm, **x** is equal to the scalar *p*O. Its value is usually defined in critical state models by the hardening law

$$dp\_{\rm O} = \frac{1+\varepsilon}{\lambda-\kappa} \, p\_{\rm O} \, d\varepsilon\_{\rm V}^{\rm P} \tag{6}$$

where *λ* is the slope of the virgin compression line of the soil, and *dε*V<sup>p</sup> is the volumetric plastic strain. For the axisymmetric conditions adopted, where *d*ε = (*dε*V, *dε* S) in which *dε* <sup>S</sup> is the deviatoric strain, **D**<sup>p</sup> can be calculated as

$$\mathbf{D}^{\mathbb{P}} = \frac{\begin{pmatrix} \left(2p - p\_{\bullet}\right)^{2} & \left(2p - p\_{\bullet}\right)\frac{2q}{M^{2}}f\_{\mathbb{V}}\\ \left(2p - p\_{\bullet}\right)\frac{2q}{M^{2}} & \left(\frac{2q}{M^{2}}\right)^{2}f\_{\mathbb{V}} \end{pmatrix}}{\left(2p - p\_{\bullet}\right)^{2} + f\_{\mathbb{V}}\left(\frac{2q}{M^{2}}\right)^{2} + \frac{\kappa}{\lambda - \kappa}p\_{\mathbb{O}}(p - p\_{\mathbb{O}})} \tag{7}$$

where *fν* is a function of Poisson's ratio given by

$$f\_{\nu} = \frac{9(1 - 2\nu)}{2(1 + \nu)}\tag{8}$$

Dividing by (*p*O) <sup>2</sup> results in the dimensionless expression

$$\mathbf{D}^{\mathbb{P}} = \frac{\begin{pmatrix} \left(2P - 1\right)^2 & \left(2P - 1\right) \frac{2Q}{M^2} f\_{\mathbb{V}}\\ \left(2P - 1\right) \frac{2Q}{M^2} & \left(\frac{2Q}{M^2}\right)^2 f\_{\mathbb{V}} \end{pmatrix}}{\left(2P - 1\right)^2 + f\_{\mathbb{V}} \left(\frac{2Q}{M^2}\right)^2 + \frac{\mathbb{x}}{\lambda - \mathbb{x}} (2P - 1)}\tag{9}$$

Additionally, **D**<sup>e</sup> is defined as

$$\mathbf{D}^{\mathbf{c}} = \begin{pmatrix} K & 0 \\ 0 & 3G \end{pmatrix} \tag{10}$$

where *K* is the bulk modulus

$$K = \frac{1+\varepsilon}{\kappa}p\tag{11}$$

and *G* is the shear modulus

$$G = \frac{3(1 - 2\nu)}{2(1 + \nu)}K\tag{12}$$

These moduli depend on the void ratio and are updated as the soil is loaded and deformed. Therefore, the elastic stiffness matrix depends on the strain path followed, as well as the preconsolidation stress, which must be updated according to Equation (6) from the volumetric plastic strain. However, Equation (9) shows that **D**<sup>p</sup> depends on *P* and *Q* only. Actually, it depends on only one parameter, since both *P* and *Q* are on the surface *F* = 0 (Figure 1), which, like **D**p, can be parametrised by the angle *θ* shown in Figure 2. Furthermore, in this figure, it is proposed to select *N* equispaced points (with an Δ*θ* angle between them) at *F* = 0, calculating **D**<sup>p</sup> <sup>i</sup> for each of them. Subsequently, when solving a boundary value problem, the value of **D**<sup>p</sup> is approximated by the interpolation **D**p\* for angle values between two previously calculated values

$$\mathbf{D}^{\mathbf{P}^\*} (\theta) = \frac{(\theta\_{\mathbf{i}+1} - \theta\_{\mathbf{i}})}{(\theta\_{\mathbf{i}+1} - \theta\_{\mathbf{i}})} \mathbf{D}^{\mathbf{P}} \mathbf{i} + \frac{(\theta - \theta\_{\mathbf{i}})}{(\theta\_{\mathbf{i}+1} - \theta\_{\mathbf{i}})} \mathbf{D}^{\mathbf{P}} \mathbf{i} + 1 \; \text{, } \theta\_{\mathbf{i}} \le \theta < \theta\_{\mathbf{i}+1} \tag{13}$$

**Figure 2.** (**a**) Discretization of the dimensionless yield surface as a function of the angle *θ*, including the position of the tension line. (**b**) Diagram of the interpolation strategy. Note that for simplicity, the interpolation of the plastic matrices is represented in one-dimensional space.

The calculation procedure proposed in this work is based on precomputing the vector of matrices {**D**<sup>p</sup> 1, **D**<sup>p</sup> 2, ... , **D**<sup>p</sup> N}, thus avoiding repeating the calculation of **D**<sup>p</sup> during the resolution of the boundary value problems.

#### **3. Evaluation of the Precomputation Efficiency**

The efficiency of the proposed procedure depends on the cost of calculating **D**p. If this cost is not an important part of the total cost of the calculation process, the performance associated with using **D**p\* might not be of interest. Figure 3a summarises the process followed to solve a computational step when integrating the constitutive model. Table 1 lists the number of operations associated with the calculation of **D**<sup>e</sup> and **D**p. Although these numbers can vary depending on the way in which the algorithm is programmed, the changes are not significant in terms of the relative computational cost of each magnitude. Using the results of the benchmark defined in Appendix A (where it is indicated that the estimated cost of each multiplication is equivalent to 1.101 sums, with the cost of division being 3.698 sums), the costs shown in the last column of Table 1 were obtained. While it is advisable to be cautious about the reliability in absolute terms of the representativeness of these numbers, they do characterise the weight that each calculation carries in relative terms. Figure 3b shows the strategy followed when applying the precom-

obtaining **D**ep\* by

$$\mathbf{D}^{\mathbf{c}\mathbf{p}\*}\left(\theta\right) = \mathbf{D}^{\mathbf{c}} \cdot \left(\mathbf{I} - \mathbf{D}^{\mathbf{p}\*}\left(\theta\right)\right) \tag{14}$$

putational procedure. As noted, **D**<sup>e</sup> must be calculated as in the conventional method. In addition, interpolation of the **D**<sup>p</sup> <sup>i</sup> matrices must be made to obtain **D**p\* by performing the operations indicated in Table 1. Using again the values of Appendix A, the computational cost of Table 1 is obtained, which is approximately 42% of the cost of calculating **D**p. The computational time savings are remarkable, and the application of the approach proposed here is therefore of interest.

**Figure 3.** Integration of a computational step using (**a**) a conventional explicit algorithm and (**b**) the precomputational approach.


**Table 1.** Equivalent sum costs of the main stages that conform to a conventional explicit integration algorithm and an algorithm with the proposed procedure.

#### **4. Precomputation Density**

To define the precomputation strategy, the minimum number of points N to be taken at *F* = 0 must be determined in such a way that the accuracy of the calculation when using **D**p\* is comparable to the accuracy obtained when calculating **D**p. To accomplish this goal, a set of constant strain rate paths with a final strain of 20% and with different ratios between volumetric and shear strains, as in Figure 4, were inspected. In all of these inspection exercises, initial spherical conditions were assumed with an effective mean stress of *p* = *p*O/2 = 200 kPa. The analysis was performed for the three soils, namely Weald clay, Klein Belt Ton, and kaolin, which are characterised in Table 2. The differences between their properties make the conclusions reached not limited to a particular material.

**Figure 4.** Strain paths *ε*α as a function of the angle *α* considered to define the number of *N* subdivisions of the yield surface.

**Table 2.** Parameters used to characterise the behaviour of the clays assumed in the models (adapted from Schofield and Wroth [15]).


As a reference, all paths were simulated using an explicit Euler integration method. An intensive substepping process was applied, i.e., the number of calculation steps was divided into smaller substeps successively until the solution differed from the previous solution by only the sixth significant digit, which is a strategy similar to that used by Sołowski and Gallipoli [16]. These results were compared with the precomputed–interpolated solution using the same number of substeps during the explicit integration.

For each path, for each value of *α* in Figure 4, the number of points N to be computed on the yield surface was determined in such a way that the normalised root-mean-square deviation between the stress solution obtained by explicit Euler integration and the solution reached by calculation with the precomputation–interpolation procedure was less than 100, 10−1, 10−2, 10−3, and 10−4. Figure 5 shows the values of N obtained for each of these tolerances and for each value of *α*.

It is important to note that, as illustrated in Figure 6, when imposing a total strain of 20%, in all cases, a significant yield was attained. In other words, trajectories that were highly sensitive to the quality of the precomputation results were analysed.

In all paths, as the tolerance decreased, N increased. This trend can be clearly seen in Figure 7, where for each tolerance value, the maximum value of N obtained for any *α* was plotted. The behaviour of N was well defined and not erratic; additionally, Figure 7 can be used to define the number of points N to be taken to avoid exceeding a certain value of tolerance.

**Figure 5.** Minimum number *N* of subdivisions of the yield surface necessary to ensure a certain level of relative tolerance for each of the *α* strain paths defined in Figure 4 for (**a**) Weald clay, (**b**) Klein Belt Ton, and (**c**) kaolin.

**Figure 6.** Ratio of elastic and plastic strain increments for the *α* paths defined in Figure 4 with respect to the number of computational steps for (**a**) Weald clay, (**b**) Klein Belt Ton, and (**c**) kaolin.

**Figure 7.** Fitting of the minimum number *N* of subdivisions of the yield surface is necessary to ensure a certain level of relative tolerance for any stress–strain path for (**a**) Weald Clay, (**b**) Klein Belt Ton and (**c**) Kaolin.

#### **5. An Inspection on Solving Boundary Value Problems**

Although the comparison in Table 1 highlights the interest of the precomputational approach, it does not illustrate its scope clearly enough, because it is associated with the simulation of individual stress–strain paths. Considering a finite element model in which conventionally the stresses are a state function, those paths would illustrate the behaviour of a single Gaussian point. It is reasonable to question what will happen when, in a real problem, the number of Gaussian points is in the hundreds or thousands.

It is not easy to answer this question. Savings will depend on the structure of the finite element code used and the calculation and information management strategies implemented. Each code will have a different response. However, in general, it is reasonable to assume that if *N*GP is the number of Gaussian points considered in the mesh, the cost of a conventional explicit integration, *χ*C, can be estimated as

$$\chi\_{\rm C} = N\_{\rm GP} \left( \eta\_{\rm Axce} \, \mathrm{G}\_{\rm De} + \eta\_{\rm Acep} \left( \mathrm{G}\_{\rm De} + \mathrm{G}\_{\rm Dp} \right) \right) \tag{15}$$

where the computational costs of **D**<sup>e</sup> and **D**p, and *C*De and *C*Dp have been defined (as well as the computational cost of **D**p\* and *C*Dp\*) in Table 1, and for the total number of computational time steps needed to solve the boundary problem, *η*Δε<sup>e</sup> defines the weighted value of the elastic increments, and *η*Δεep defines one of the elastoplastic increments. In turn, the cost of the precomputation–interpolation strategy, *χ*PI, is estimated as

$$\chi\_{\rm PI} = N\_{\rm GP} \left( \eta\_{\Delta \varepsilon \rm e\*} \, \mathcal{C}\_{\rm De} + \eta\_{\Delta \varepsilon \rm e\*} \left( \mathcal{C}\_{\rm De} + \mathcal{C}\_{\rm Dp\*} \right) \right) \tag{16}$$

where *η*Δεe\* and *η*Δεep\* are the values analogous to *η*Δε<sup>e</sup> and *η*Δεep, respectively, when the integration of the constitutive model is done using the precomputation and interpolation procedure. Their values can, in principle, be different. However, since the discretization of the yield surface is proposed with a very demanding tolerance (10−<sup>4</sup> in the analyses performed above), it is to be expected that *η*Δεe\* ≈ *η*Δε<sup>e</sup> and *η*Δεep\* ≈ *η*Δ<sup>ε</sup>ep. In this case, *ρ*, the ratio between the cost of a traditional explicit integration and the cost of the precomputation–interpolation strategy, is given by the expression

$$\rho = \frac{\chi\_{\rm Pl}}{\chi\_{\rm C}} = \frac{\eta\_{\rm Acc} \,\mathrm{C\_{De}} + \eta\_{\rm Accep} \left(\mathrm{C\_{De}} + \mathrm{C\_{Dp\*}}\right)}{\eta\_{\rm Acc} \,\mathrm{C\_{De}} + \eta\_{\rm Accep} \left(\mathrm{C\_{De}} + \mathrm{C\_{Dp}}\right)} = \frac{1 + \frac{\eta\_{\rm Laser}}{\eta\_{\rm Laser}} \left(1 + \frac{\mathrm{C\_{Dp\*}}}{\mathrm{C\_{Dp}}}\right)}{1 + \frac{\eta\_{\rm Laser}}{\eta\_{\rm Laser}} \left(1 + \frac{\mathrm{C\_{Dp}}}{\mathrm{C\_{Dp}}}\right)}\tag{17}$$

The relative cost is independent of the mesh size. Furthermore, it is necessary to add to the computational costs considered in Equations (15) and (16) the costs associated with the generation of the database {**D**<sup>p</sup> 1, **D**<sup>p</sup> 2, ... , **D**<sup>p</sup> N}. Consequently, it is to be expected that if the mesh is small, the relative importance of the precomputation cost is not negligible. However, the proposed methodology is oriented to large-scale computation, and thus it is a consistent estimator of the efficiency, and for a large computational volume, as shown in Equation (17), it is independent of the mesh size. The ratio *η* = *η*Δ<sup>ε</sup>ep/*η*Δε<sup>e</sup> expresses the relative importance of the computation of plastic steps in a problem. As expected, given that the precomputation is aimed at optimizing the calculation of the plastic steps, its efficiency is conditioned by the relevance of the plastic strains over the elastic strains, in other words, by the value of *η*. Using the values of *C*De, *C*Dp, and *C*Dp\* from Table 1, the change in *ρ* with respect to *η* was plotted and is shown in Figure 8. If *η* = 0, when the yield is not reached, precomputation is not applied. Only elastic steps are calculated, and both methods are equal (*ρ* = 1). However, as seen at the end of Section 3, as soon as some plastic steps are required, the precomputation is more efficient (*ρ* < 1). The larger the volume of the calculation of the plastic steps, the higher the efficiency. When the entire calculation is in the plastic regime (i.e., in normally consolidated conditions), the efficiency tends to a maximum value of 54%.

**Figure 8.** Dependence of *ρ* (ratio between the cost of a traditional explicit integration and the cost of the precomputation–interpolation strategy) on the value of *η* (relative importance of the computation of the plastic steps over the elastic steps).

#### **6. Conclusions**

A methodology was presented to parameterise, precompute, and approximate the plastic behaviour of a soil by interpolation. The procedure can be applied to materials whose constitutive models allow the parametrisation and subsequent normalisation of their yield surface as well as the plastic components of the elastoplastic matrix. In this work, the modified Cam Clay model was used as a reference model. After describing the precomputation and interpolation procedure, the efficiency of the method was evaluated by estimating the computational cost of one computational step. Then, a criterion was defined to determine the density of the precomputed database so that the quality of the approximations obtained with the proposed method is like that of solutions based on conventional explicit integration methods. Finally, an inspection was conducted to assess the efficiency of the proposed method in the solution of boundary value problems. It was concluded that the proposed method can reduce the computational burden to as much as 54% of the conventional cost, a significant savings when solving problems in technological applications. This finding is of particular interest in the case of modelling large-scale engineering systems, in sensitivity analysis, or in solving optimisation problems, where a large volume of calculations is required.

**Author Contributions:** Conceptualization, V.N.; methodology, V.N.; software, D.G.; validation, D.G. and G.D.l.M.; formal analysis, Á.Y.; investigation, L.A.; resources, L.A.; data curation, G.D.l.M. and V.N.; writing—original draft preparation, D.G. and G.D.l.M.; writing—review and editing, V.N., V.C., and Á.Y.; visualization, V.C.; supervision, V.N.; project administration, V.N.; funding acquisition, V.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by Junta de Comunidades de Castilla-La Mancha and the European Regional Development Fund (European Union) through the project SBPLY/19/180501/000222.

**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.

#### **Appendix A. Characteristics of Computers**

To obtain a reasonable estimation of the computational cost of the operations (multiplications and divisions) indicated in Table 1, a benchmark was conducted for each type of operation with six computers (see Table A1) in such a way that the cost of the operation was converted into equivalent sums. Calculation loops were performed until the CPU time was stationary in relation to a similar addition loop for each computer and for both multiplication and division. Thus, calculating the average of the results obtained with the computers used, the cost of multiplication was estimated to be equivalent to 1.101 sums, and the cost of the division was equal to 3.698 sums.

**Table A1.** Characteristics of the computers used in the benchmarks to quantify the different computational costs.


#### **References**

