An Efficient Model for Numerical Simulation of the Mechanical Behavior of Soils. Part 2: Applications

. A numerical model to simulate the mechanical behavior of soils was introduced in Part 1 of this paper (also pub-lished in this issue). Detailed information about the analytical model were presented, where the critical state theory for soil mechanics was considered in the context of the elastoplastic formulation. Moreover, an efficient numerical formulation to deal with nonlinear applications was also presented, featuring important characteristics such as reduced integration techniques, explicit integration of the constitutive equation and a corotational formulation for the kinematical description of the continuum. In this second part of the present work, the numerical model proposed in the previous paper is applied to some classical examples of soil mechanics to demonstrate the applicability of the present formulation. Effects of a geometrically nonlinear approach over the numerical predictions are investigated and comparisons are performed taking into account results obtained by using a geometrically linear model. In addition, some comparisons are also carried out considering evalu-ations performed with different constitutive formulations in order to observe the mechanical behavior of the soil mass under different constitutive assumptions.


Introduction
In the first part of the present work (Braun & Awruch, 2013) a numerical model to simulate the mechanical behavior of soil masses was presented. An analytical model based on critical state soil mechanics (CSSM) was formulated utilizing the elastoplastic framework, which is entirely justified by experimental evidences. A corotational formulation was also presented in order to improve the kinematical description of the continuum. In addition, a finite element model was developed considering eight-node hexahedral elements with one-point quadrature and special techniques for the integration procedure of the constitutive equation. It was observed that the formulation presented previously leads to a highly efficient numerical model, especially for nonlinear analysis, where the computational effort is significantly reduced when compared with other formulations employing standard approaches.
In the present paper some classical problems, including three-dimensional examples, are numerically simulated to demonstrate the applicability of the present formulation when geotechnical analyses are carried out. Effects of a geometrically nonlinear approach over the numerical predictions are investigated considering comparisons performed with respect to results obtained by using a geometrically linear approach for the numerical model proposed here. Furthermore, investigations are also performed taking into account different constitutive formulations in order to compare responses obtained under different constitutive assumptions for the applications analyzed in this work.
Whenever possible, predictions obtained here are compared to results obtained by other authors.

One element triaxial test
A triaxial compression test is numerically simulated here using a single hexahedral element in order to perform a simple and reliable verification over numerical routines of the present model referring to elastoplastic analysis and Cam-Clay formulation. Information about geometrical parameters and boundary conditions employed in the present analysis are shown in Fig. 1 and material properties related to Cam-Clay characterization of the soil mass are found in Table 1. In order to reproduce the experimental procedure, vertical displacements are imposed gradually at the top of the element and horizontal forces are applied on the vertical faces and calculated based on an initial hydrostatic stress state. The initial conditions are mainly characterized according to the overconsolidation ratio (OCR), which may be defined as follows: where p 0 and p c0 are initial values for effective normal stress p and preconsolidation pressure p c . Two different initial conditions are analyzed here: a lightly overconsolidated state, characterized by OCR = 1.2, and a strongly overconsolidated state, characterized by OCR = 6.0. In addition, comparisons between the geometrically linear and geomet-rically nonlinear approaches are also performed to determine effects of geometrical nonlinearities over the computed results. Results obtained in this work are presented in Fig. 2 and compared with numerical predictions performed by Sheng et al. (2000). A good agreement can be observed between the reference work and predictions carried out with the geometrically linear approach. Hardening and softening behaviors are clearly reproduced, where hardening is associated to the lightly overconsolidated state as well as softening is related to the strongly overconsolidated state. For the lightly overconsolidated state, the stress path intersects the yield surface and then returns towards the critical state line following approximately the same path. The critical state is obtained at e a @ 20% and e a @ 50% for OCR = 6.0 and OCR = 1.2, respectively, which is in accordance with the reference results. On the other hand, when the geometrically nonlinear model is utilized, the critical state is not ob-tained. The stress path related to the lightly overconsolidated state returns from the point of intersection on the yield surface using a little different path, which follows an asymptotic curve towards the critical state line (CSL). The same behavior is observed for the strongly overconsolidated state after yielding. That particular behavior leads to a continuous increase of the deviator stress, which violates one of the basic critical state assumptions. However, it is worth to mention that geometrically nonlinear models tend to increase stiffness when deformation takes place, which is originated from stress rate terms added to the elastic constitutive matrix (Braun & Awruch, 2013 -Eq. 25) and changes in the reference configuration.

Embankment analysis
An embankment under the action of gravity is analyzed in this section. In order to verify the present numerical scheme, the Mohr-Coulomb constitutive model is first considered to compare results with numerical predictions performed by Zienkiewicz et al. (1978). Simulations are then carried out considering the Cam-Clay formulation utilized in this paper and a geometrically linear model, whose results are then compared with the previous predictions.
Geometrical characteristics of the computational domain employed in this example are shown in Fig. 3, where boundary conditions for the displacement field are also illustrated. Displacements in the orthogonal direction (axis z) are locked to reproduce the plane strain condition. The physical parameters adopted in the present simulations are found in Table 2 and all analyses are performed using a finite element mesh with 459 eight-node hexahedral elements. It is important to notice that material properties for  the Cam-Clay constitutive model are not available in the reference work. Therefore, a set of physical constants usually employed by other authors to characterize a Cam-Clay modeling of the soil is considered here. In order to reproduce stiffness characteristics of the Mohr-Coulomb simulations, a hydrostatic stress state p is imposed over the computational domain, which is calculated using Eq. 3 from the first part of the present work (see Braun & Awruch, 2013) and the classical elastic relation: where E is the Young's modulus, K is the compressibility modulus and n is the Poisson's ratio. Consequently, the preconsolidation pressure can be defined based on the overconsolidation ratio OCR (see Eq. 1), which assumes the following values in this study: OCR = 1.025, OCR = 1.05, OCR = 1.075, OCR = 1.1 and OCR = 2.0.
All results obtained in the present work are summarized in Table 3 together with predictions presented by Zienkiewicz et al. (1978), where comparisons in terms of vertical (u y ) and horizontal displacements (u x ) measured at point A (see Fig. 3) are performed. In addition, these results are plotted against load ratio Qn/Qt in Fig. 4, where Qn and Qt are the current and total load applied on the embankment, respectively. All simulations carried out using the Mohr-Coulomb constitutive model resulted in predictions very similar to those obtained by Zienkiewicz et al. (1978). By using the Modified Cam-Clay formulation, a reasonable agreement with the Mohr-Coulomb results can be observed for analyses performed with OCR ³ 1.075, for which identical predictions are produced owing to the fully elastic state established in this range of overconsolidation ratios. On the other hand, as OCR decreases the horizontal displacements u x tend to follow a distinct curve when compared with that provided by using the Mohr-Coulomb formulation with smaller cohesion values, where negative Soils    displacements are observed if c < 7 kN/m 2 is considered. Nevertheless, all curves referring to vertical displacements u y obtained with the present Cam-Clay model and decreasing OCR values show a similar trend to that verified in the Mohr-Coulomb simulations when the cohesion is reduced.
In Fig. 5, final patterns of flow vectors and stress fields obtained from the Mohr-Coulomb and Modified Cam-Clay models are presented. In order to confront simulations performed with different constitutive models but with similar characteristics, comparisons are carried out considering that responses obtained with the Cam-Clay formulation are referred to OCR = 1.075 and responses obtained with the Mohr-Coulomb model are referred to c = 20 kN/m 2 . As can be noticed, a good agreement can be verified in all comparisons performed here. It is worth to mention that stress measures from the Cam-Clay responses are referred to sign conventions usually adopted in soil mechanics, where compression is considered as positive as well as tension is taken as negative. Moreover, Cam-Clay axial stresses s xx and s zz are shifted by an initial value related to the hydrostatic stress state imposed at the start of the simulation.

Two-dimensional analysis of a rigid strip footing
A two-dimensional analysis of a rigid strip footing resting on a soil layer is carried out in this section. Footing tests have been extensively used to validate numerical algorithms based on elastoplastic formulations applied to soil mechanics. This problem may be considered as a very difficult task for most of the numerical schemes due to strong rotations verified in the principal stresses and the singularity observed near the edge of the footing.
All simulations are performed using the Modified Cam-Clay model and results are compared with numerical predictions obtained by Sheng et al. (2000), who also utilized a Cam-Clay formulation. Comparisons are also carried out with respect to the geometrically linear and nonlinear approaches adopted in this problem. Geometry and boundary conditions utilized in the present study are shown in Fig. 6 and material parameters adopted in the Cam-Clay characterization of the soil mass are listed in Table 4. The computational domain is divided into 2703 eight-node hexahedral elements and the plane strain condition is imposed by displacement restrictions applied in the orthogonal direction (axis z). In order to simulate the action of the rigid foundation on the soil mass, uniform vertical displacements are imposed along the footing-soil interface. It is important to notice that critical state models cannot mobilize any strength if the effective normal stress p is zero or under tensile stress. Therefore, a compressive initial stress state must be established over the soil, which is accomplished in this example by considering a hydrostatic stress state generated from the soil specific weight. Moreover, the soil is assumed to be overconsolidated to 50 kPa at the ground surface, which implies that, at the ground surface, the preconsoliation pressure is 50 kPa.
The load-displacement curves obtained in this work are presented in Fig. 7, where a prediction computed by Sheng et al. (2000) is also shown. Displacements and footing loads are measured at point A underneath the footing, as indicated in Fig. 6. A good agreement is verified when predictions provided by the geometrically linear model is taken into account. On the other hand, results obtained us- 174 Soils and Rocks, São Paulo, 36 (2)  ing the geometrically nonlinear approach lead to an overestimation of loads owing to an increase of stiffness, which is a characteristic usually observed in geometrically nonlinear models. Yielding can be observed in all predictions when the imposed displacements reach almost 20% of the footing width.
Stress paths computed at points A and B underneath the footing (see Fig. 6) are found in Fig. 8   sons are performed with results presented by Sheng et al. (2000). In the present work, those points are referred to the Gauss points of elements localized underneath the footing borders. A good correlation can be observed when both the linear and nonlinear responses obtained at point A are compared with the respective result presented by the reference work. As can be noticed, predictions obtained with geometrical linear and nonlinear models lead to very similar stress paths in the elastic range. However, some mayor differences between linear and nonlinear responses can be observed after yielding. At point B, a reasonable agreement is verified between the geometrically linear stress path obtained here and the corresponding stress path predicted by Sheng et al. (2000), which is not verified when the geometrically nonlinear response is considered after the onset of the plastic yielding. Unlike point A, where the stress distribution is more regular, the stress distribution near point B presents strong gradients, which demands a precise localization of the measure point. Since the measure points are not precisely defined in the reference work, the measure points utilized in this work are probably localized at different positions if compared with those adopted by Sheng et al. (2000).
In Fig. 9, mesh configurations and stress distributions over the computational domain are shown, which are referred to the final configuration of the footing settlement. Comparisons are performed with respect to predictions obtained using the geometrically linear and geometrically nonlinear approaches, where a reasonable agreement is verified for s zz .

Two-dimensional analysis of flexible strip footings
In the present analysis, numerical simulations of soil layers under the action of two-dimensional flexible strip footings are performed, where two distinct configurations referred to load and geometrical characteristics adopted in the numerical modeling are considered. Flexible strip footings may be represented applying uniform vertical loads along the footing-soil interface. The first configuration considered here is analyzed using the Mohr-Coulomb and Modified Cam-Clay models in order to compare predictions performed by Zienkiewicz et al. (1978), who utilized the Mohr-Coulomb formulation to determine settlements for a strip load on a soil layer with an ideal weightless material. The second configuration was proposed by Borja & Tamagnini (1998) to verify their critical state formulation, which is based on an alternative Cam-Clay model. In this second problem, comparisons of results obtained with the geometrically linear and nonlinear approaches utilized in the present scheme are also carried out.
Geometrical characteristics referring to the configurations employed in the present study are found in Fig. 10, where boundary conditions are also shown. Physical properties of the soil are presented in Table 5 according to the constitutive models and configurations utilized here. The computational domain is divided into 480 and 540 eightnode hexahedral elements, considering the first and second configurations respectively. For simulations referred to the first configuration, an initial hydrostatic stress state p 0 is imposed over the computational domain using the same procedure adopted in Section 2.2, whereas physical constants must be utilized based on the properties indicated in Table 5. Just as in the previous analysis, the initial stress state for the second configuration is generated using the soil specific weight. Plane strain conditions are imposed for both configurations.
Load-displacement curves computed in the present work are presented in Fig. 11 considering both configurations simulated here, where results obtained from the reference authors are also shown. Several overconsolidation  ratios (OCR = p c0 /p 0 ) are simulated here in order to reproduce the reference predictions and investigate its influence over the mechanical response. All displacements are measured at a point corresponding to the point of contact between the soil layer and the intermediate point underneath the footing, which correspond to the top left node of the computational domain.
Results obtained for the first configuration show that simulations with a conventional Mohr-Coulomb formula-tion cannot converge in this problem. It is worth to notice that Zienkiewicz et al. (1978) only obtained a stable solution by using an alternative Mohr-Coulomb model. Comparisons are also performed using the Modified Cam-Clay model with various OCR's, where the elastic responses are well correlated to the Mohr-Coulomb results. Nevertheless, the plastic response is quite different, as expected, since each formulation considers different assumptions to describe the plastic phenomena. Even so, some simulations Soils   approximated the Mohr-Coulomb responses, such as the curves referred to OCR = 1.1 and OCR = 1.2. Regarding results obtained for the second configuration, the best agreements between the reference work and simulations carried out here are referred to OCR = 2.5 and OCR = 2.25, although the plastic behavior immediately after the onset of the plastic yielding in these cases are distinct to that predicted by Borja & Tamagnini (1998), which is explained by some particular characteristics observed in the constitutive model adopted in the reference work. All simulations performed here reproduced well the mechanical response in the elastic range, but only that referred to OCR = 1.5 obtained the yielding load indicated by the reference result. As expected, the prediction obtained with the geometrically nonlinear model resulted in smaller vertical displacements when compared with that obtained with the linear model.
Final mesh configurations and stress components computed with the numerical simulation carried out over the second configuration are presented in Fig. 12, where comparisons are performed taking into account results obtained with the geometrically linear and nonlinear models and OCR = 2.5. The mesh deformation referred to the nonlinear analysis is clearly smaller than that obtained with the linear model, although both predictions show expressive lateral displacements in the region beneath the footing location.

Three-dimensional analysis of compaction and footing problems
The numerical model proposed in this work is verified in this section using numerical tests that demonstrate its ability to predict results for three-dimensional stress states. Two different problems are analyzed here, which are characterized with different geometrical properties and distinct load characteristics. The first problem is based on the footing test performed by Lee et al. (2005), where a cubical model under the action of a flexible square footing is inves- 178 Soils and Rocks, São Paulo, 36 (2)  tigated. The second problem is referred to the compaction test proposed by Xia & Masud (2009), who analyzed spatial density changes induced during compaction of a threedimensional soil model. The Modified Cam-Clay formulation proposed in this work is utilized in both the problems analyzed here. A comparison between results obtained with the geometrically linear and nonlinear models is carried out for the compaction test. Geometrical and load characteristics related to the problems studied in this section are shown in Fig. 13. For both geometrical models, symmetry boundary conditions are imposed on the planes defined by (x, y, z) = (0, 0, 0). The computational domain is composed of 24x24x14 and 14x14x8 eight-node hexahedral elements, considering the footing and compaction tests, respectively. Physical constants adopted in the present simulations are listed in Table 6. Gravity and external loads are applied simultaneously in the footing test and initial stress states utilized in both analyses carried out here are based on the procedure explained in Section 2.2, using the Young's modulus indi-cated in the respective reference works and other physical constants presented in Table 6.
Load-settlement curves obtained in the present simulations are shown in Fig. 14 together with predictions taken from the respective reference works, where displacements are measured at (x,y,z) = (0,0,2) and (x,y,z) = (0,0,11) for the footing and compaction test, respectively. In the elastic range, the response obtained in the footing test with the present formulation agrees very well with the numerical predictions observed in Lee et al. (2005). On the other hand, after plastic yielding a slightly different response is verified, although a reasonable agreement is still observed. For the compaction test, several overconsolidation ratios (OCR) were analyzed. Results show that when OCR = 3.0 is considered, a very good agreement is verified in the range 0 £ q £ 16 between the present predictions and those referred to Xia & Masud (2009). However, the plastic behavior is quite different in this case as well as in the remaining cases studied here. It is worth to mention that Xia & Masud Soils  (2009) adopted a constitutive formulation based on the smooth surface cap model. In the comparison performed with respect to the influence of geometrical nonlinearity, predictions obtained with the geometrically nonlinear model and OCR = 3.0 indicate an anticipation of the yielding load and a modification of the mechanical response in the plastic range when compared with the geometrically linear response under the same conditions. Figure 15 illustrates the final mesh configuration and final stress distribution of s zz over the computational domain referred to the footing test. An intermediate stress distribution of s zz and the corresponding mesh related to the compaction test at the fictitious time t = 60% are also   shown. As expected, the characteristic compression bulb under the load area is easily identified in the stress field s zz for both results. The stress distribution and the mesh configuration obtained in the compaction test show a good agreement with the corresponding results obtained by Xia & Masud (2009).

Conclusions
A numerical model based on one-point quadrature and critical state formulation was proposed in a previous paper to simulate the mechanical behavior of soil masses. Results obtained here demonstrated good agreement with the reference works whenever comparisons were possible to be made, considering that experimental predictions obtained by other authors were also included. Investigations performed with respect to geometrical nonlinearity showed that an approach considering large displacements and rotations may be not important for the cases analyzed here. On the other hand, it is well known that many geotechnical problems are characterized by large strains in the plastic range. Comparisons considering different constitutive models were also carried out, where responses with different behavior were demonstrated. Although results concerning the numerical efficiency of the present formulation were not presented here, simulations indicated that the stress integration scheme utilized in this work presents high sensibility to the tolerance values adopted in the iterative procedures. These conclusions are similar to those obtained by Sloan et al. (2001) and, therefore, tolerance values should be chosen carefully (see Sloan et al., 2001 for suitable tolerance values). In order to improve the present model a finite strain formulation to solve problems involving large deformation should be implemented in future works. An alternative to be considered may be found in Nazem et al. (2008), where an arbitrary Lagranian-Eulerian algorithm is proposed.