Effect of using different approximation models to the exact Mohr–Coulomb material model in the FE simulation of Anchor Foundations in sand

Singularities in the Mohr–Coulomb (MC) yield criterion create numerical difficulties when implemented in the finite element model (FEM). In this research, a two-dimensional (2D) and a three-dimensional (3D) elasto-plastic FEM coupled with isotropic strain hardening–softening or simple strain-softening law incorporating non-associated flow rule including shear band, is used to investigate the effect of singularities, embedment, shape, soil frictional as well as dilation angle on vertically uploaded shallow anchor foundations buried in Toyoura sand. Mohr–Coulomb–Drucker Prager (MC–DP) model provides a more appropriate solution than Drucker Prager–Drucker Prager (DP–DP) model. Besides, the peak resistance factor and the settlement at the ultimate capacity have been found as the functions of embedment, shape, and sand density. Also, the numerical model gives satisfactory agreement with the previous studies. In particular, a simple strain-softening model cannot capture the settlement behavior of anchor foundations in sand efficiently.

exact MC model considering anchor foundations buried in the sand of different densities ranging from dense to loose. Drucker-Prager (DP) model is used as a potential surface for fitting the MC law because of the availability of different cone parameters in the DP model to approximate the MC hexagon on the π-plane. The MC material model can be fitted into the inner and outer apices [39]. Also, design graphs are developed by some parametric studies on shapes, embedment, frictional, and dilatation angles. The results of numerical predictions are also compared with the classical bearing capacity equations.

Experimental procedures
The experimental setup used for this study is shown in Fig. 1. A cylindrical container having a diameter of 590 mm is used as a soil bin to neglect the boundary effects. At the time of experiments, two types of flat anchor foundations (circular and rectangular) are considered. The diameter (D) of the circular anchor for this study is 10 cm at an embedment ratio (H/D, where H = depth of sand mass) of 2 having the plate thickness of 5 mm. However, the width (B) of the rectangular anchor is 5 cm at an aspect ratio (L/B, where L = length) of 1 to 6 and embedment ratio (H/D) of 2. Each anchor is installed at the bottom of the soil bin. The data of vertical pullout load at pulling out speeds of 0.003 cm/ min and upward displacement are recorded by a load cell, connected with an anchor rod and a displacement transducer is installed on the top of the anchor rod respectively with the help of a computerized data acquisition system. All movement of the anchor is operated by a direct current (DC) motor until the residual conditions are achieved. Airdried uniform Toyoura sand (specific gravity ( G s ) = 2.64, effective diameter of particle ( D 50 ) = 0.016 cm, maximum void ratio ( e max ) = 0.98, minimum void ratio ( e min ) = 0.61, and no fine content less than 0.00075 cm) is used in this test. The air pluviation method is used to build up the ground above anchor foundations with the help of two sieves through the raining device. The experiment is conducted on the sand of different uniform densities which are achieved by differing falling height and deposition density. To evaluate the uniformity of the sand mass, the density of the ground is measured by using six 100 cc cylindrical samplers (D = 5 cm), which are placed on the bottom of the soil bin. In this study, three different types of sand (dense having a density of 1630 kg/ m 3 ; medium having a density of 1480 kg/m 3 ; loose having a density of 1350 kg/m 3 , error allowance = ± 10 kg/m 3 ) having a relative density of 95% to 5% are used. The experimental failure mechanisms are presented in Fig. 2a, considering the circular anchor [3,12] and Fig. 3a, considering the rectangular anchor.

Numerical modeling
Four node quadrilateral isoparametric finite elements are used for analyzing the circular anchor problems as illustrated in Fig. 2b for 2D FE analysis. The mesh extends at least at a distance of four times the diameter of the circular anchor from the edges to neglect the boundary effects. In the center part of the model, the discretization is finer with quadrilateral elements than at the edges, where finite deformation is expected. To obstruct out-of-plane displacements of the central vertical symmetrical boundaries, the mesh base is fixed in all three coordinate directions, except in the anchor plate area, and zero-displacement boundary conditions are introduced. Besides, along the anchor boundary until failure, differential quasi-elastic displacement is applied at small sequential increments and the corresponding nodal forces are averaged concerning the displacement control nodes to ascertain the pullout capacity (P). In case of 3D FE simulations, first-order cubic hexahedral finite elements are used for analyzing the problems of a rectangular anchor. By taking the benefits of symmetry, one-quarter of the domain is discretized into finite elements as illustrated in Fig. 3b. The discretization is finer in the central zone (zone above the anchor and extending to H from the anchor edges, using 5 mm sided cubic hexahedral elements) than the edges due to the same reason as described for 2D FE analysis. The mesh extends at least at a distance of two times the depth of sand mass from the edges of the rectangular anchor and one time above the anchor to neglect the boundary effects [7]. More elements are used for longer anchors in the longitudinal direction to maintain the uniformity of element size across the models, for example, square anchor with L/B = 1 consists of 12,288 elements, and the rectangular anchor with L/B = 6 consists of 21,888 elements (which is approximately 1.5 times more than that of the square anchor). The normalized displacement factor, δ/D or δ/B, and the pullout resistance factor, N p (= δ d HA, δ d is the effective unit weight of dry sand and A is the anchor surface area) are ascertained by following the recommendations of Lindsay [28]. The numerical models use the nonlinear elasticity [21], non-associated flow rules (DP potential and MC yield surface), sophisticated strain hardening-softening [44,52] or simple strain softening constitutive law [56] in an elasto-plastic framework including the shear band [44]. Both types of material models are available for 2D analysis while only the simple strain-softening model option is available in the personally developed FORTRAN program for 3D analysis. An explicit dynamic relaxation (DR) method [45,52] devised with the generalized return mapping algorithm [47] is used due to its robust ability to handle bifurcation problems of extremely nonlinear material; due to its explicit nature needs very small data storage and easy to be coded; use of approximate critical damping for linear isoparametric elements with reduced integration can inherently inhibit the hourglass modes at failure. This study focuses on the effect of different approximate models to the exact MC material model in the finite element method, such as DP approximation (Table 1); Argyris et al. [2]; and rounding parameter (lode angle) on the anchor (Fig. 4). The lode angle, θ can be introduced as

Table 1 Drucker-Prager approximation to Mohr-Coulomb model
If Drucker-Prager function is chosen as the potential function for its differentiability with other or Mohr-Coulomb yield criteria (generalized elasto-plastic model), in that case, φ must be replaced by dilatation angle, ψ

Parameters Compression
Plane strain cone Compromise cone where J 3D = third deviatoric stress invariants and J 2D = second deviatoric stress invariants.
The MC and DP yield functions can be written in a general form using a function, g(θ), which determines the shape of the yield surface on the deviatoric plane as and the function, g(θ) , is such that g π 6 = 1 , like shape function, and this function determines the shape of the π-section. For instance, for MC, it becomes If now the values of α, γ, and g(θ) are put in Eq. (2), it will regenerate the invariant form of the MC yield function. For smoothing the corners at θ = ± π 6 in the MC's π-section as defined by Eq. (3), this function can be written as [2], graphic is shown in Fig. 5.
where g − π 6 = K and to match with MC failure criteria, the value of K can be calculated using Eq. (3) as If K = 1, Eq. (4) will generate g(θ) = 1, which leads the circular π-section and no dependence of third stress invariant (i.e. DP model).
To find the inscribed circle, the value of θ at which this circle is tangential to the MC hexagon can be found by differentiating Eq. (3) for θ and setting the result to zero to have Substituting this value in Eq. (3), we can have the value of g(θ) of the inscribed circle as shown in Fig. 5. According to Fig. 5, extension cone and compression cone are defined to match the MC model either in triaxial extension or triaxial compression. Compromise cone is the average of compression cone and extension cone. However, an internal cone is inscribed into the MC [57]. Besides, another popular DP approximation to the MC criterion is obtained by forcing both criteria to predict identical collapse loads under plane strain conditions.

Different approximations to the Mohr-Coulomb material model
The material parameters used for the analysis are summarized in Table 2. The effect of intermediate principal stresses is neglected in the MC criterion, which contradicts the fact that uniaxial compressive strength is always smaller than biaxial compressive     shows the different approximations to the exact MC material model in the finite element method using a sophisticated strain hardening-softening model at H/D = 2 and D = 10 cm. It is observed that MC-DP model, where MC represents yield surface and DP is the potential failure surface, gives better results compared with the Mohr-Coulomb-Mohr-Coulomb (MC-MC), Argyris-Argyris, Drucker Prager-Drucker Prager (DP-DP) (different matching) material model (in all cases, the preceding portion indicates yield surface and the latter part points out potential failure surface) which proves its appropriateness to use it in this study.

Effect of rounding of Mohr-Coulomb material model
As the singularities or the corners and the apex of the MC yield surface caused problems in the numerical implementation, an approximate yield surface with smoothed or rounded corners has to be used (Fig. 4). The singularities in the yield surfaces, where the gradient concerning the stresses is undefined, occur at θ = ± 30°. According to Borst [6], the solution is difficult to obtain in MC condition. To verify this, MC condition with pseudo-equilibrium element is used. The material parameters for this study are already shown in Table 2. Figure 7 that presents the effect of the rounding parameter ( θ T ) considering the anchor problem at the sand of different densities. In the MC-DP model, the solution is obtained even at the time of the Lode angle limit being 29.0°, and even the result is better than that obtained using the DP-DP model. However, the success of using the popular DP yield surface highly depends on the proper selection of a matching parameter to agree with the MC yield surface condition.

Parametric study
The effect of embedment, shape, frictional, and dilatation angle on the pullout capacity of circular or rectangular anchor foundations buried in the sand of different densities have been studied through model tests and extensive FEM analysis. Figure 8 presents the numerical relationship between the pullout resistance and displacement factor for different H/D in different layers of sand using the material parameters shown in Table 2. All curves show three distinct phases: the initial phase (the uplift resistance increased with the displacement of the anchor), followed by the decreasing tendency of uplift resistance with the displacement of anchor, and finally, the uplift resistance remains unaltered with the further uplifting of the anchor (residual state). It is observed that the uplift capacity and displacement factor depend on H/D and H. However, the numerical passive plastic zone is more developed at peak load at the higher embedment (Fig. 9). Besides, peak load mobilization at lower displacement (stiffer response) is noticed at the smaller embedment. Furthermore, peak resistance and displacement factor increase with the increment in embedment ratio, whereas peak resistance factor increases and displacement factor decreases with the increase in density of sand, which is attributed to the increase in confining pressure (Fig. 10).

Shape effect (3D modeling)
The width (B) of the foundation buried in dense Toyoura sand is used in this study is 5 cm at an aspect ratio (L/B) of 1 to 6 and embedment ratio (H/B) of 2. The material  Table 3. However, a satisfactory agreement is found between the experimental and numerical pullout resistance-displacement factor relationships. The rate of softening after the peak is decreasing with the increase in aspect ratio (Fig. 11). The relationship between mobilized peak resistance factor ( N pu ) and L/B considering rectangular anchor buried in the sand of different densities is presented in Fig. 12. The shear resistance along vertical planes plays an important role in increasing resistance factors of anchor foundations with a lower aspect ratio. The peak resistance factor is a combined function of L/B, B, and H/B as well as relative density, which violates the non-robust empirical shape factor used in traditional design. However, it decreases rapidly with the increase of L∕ B. Moreover, displacement is not well predicted due to the use of a simple strain-softening model.

Effect of peak frictional and dilation angle (using simple strain softening model)
Dilatancy of sand depends on density and stress level. When the density of sand is dense with low-stress level, it demonstrates shear dilation and when loose with the high-stress level it often exhibits shear contraction. The dilation angle perspective to zero degrees relates to a soil distorts plastically with zero volume change, which is a reasonable supposition for loose sands [41]. Here, the dilation angle is assumed from 0° to 45° in  Dilation angle, ψ 20 10 0 Residual frictional angle, φ r (°) 33 33 33 Poisson's ratio 0. Comparison with past studies Figure 15 shows the different studies to determine the peak resistance factor as a function of H/B in dense sand. Besides, Fig. 16 shows different studies on the shape effect of rectangular anchor in dense sand at H/B = 2. The dimensionless shape factor S can be defined as 28    The effect of shape diminishes with the increase of aspect ratio, gradually decreased to reach the strip anchor condition at L B = ∞ . The experimental result [53] of the analogous trapdoor with the same width is used for the strip anchor. According to Fig. 17, with the increase of peak frictional angle, the peak resistance factor is increased significantly.