Introduction

Anchors and anchoring system usage techniques are widely used in the world in different fields of civil engineering to resist vertical or horizontal uplift forces or overturning moments acting on the structures, such as transmission towers, sheet piles, retaining walls, deepwater offshore developments, airport hangars, wind loads on tall structures, buoyancy forces on buried underwater pipelines, earthquake, ice forces at different embedment, sizes and shapes [18]. Since the pullout capacity of anchors is greatly influenced by anchor geometry and local soil conditions, many studies have focused on these important factors [16, 17, 19, 20, 22, 35, 51].

However, the settlement behavior of anchor foundations is overlooked by most of the researchers due to the shortcomings of the models to predict the real soil behavior. Furthermore, traditional design approaches are based solely on the peak friction angle (\(\upphi_{\text{p}}\)) and principle of associated flow (AF) rule, which overestimate the foundation capacity in soils [10, 13, 17, 24, 30, 34, 36, 46, 49, 55]. According to the study of Bolton [5], the \(\upphi_{\text{p}}\) is controlled by the critical state friction angle (\(\upphi_{\text{c}}\)) and dilation angle (ψ). Besides, several studies have shown that the value of \(\upphi_{\text{c}}\) and the relationship between (\(\upphi_{\text{p}} - \upphi_{\text{c}}\)) and ψ depend on the sand type [8, 29, 48]. Therefore, \(\upphi_{\text{p}}\) alone cannot uniquely predict the uplift capacity of anchors in sand. Moreover, Kumar and Kouzer [26] have studied the effect of embedment ratio (H/D) and \(\upphi_{\text{p}}\) on the pullout capacity of plate anchors in the sand and found that the pullout capacity increases with the rise of H/D and \(\upphi_{\text{p}}\). Besides, Rowe and Davis [41] used FEM to analyze a strip anchor in sand and found that the dilatancy could significantly influence the ultimate pullout resistance.

Despite the availability of complex soil models for advanced applications, Mohr–Coulomb (MC) yield criterion is widely used in the elasto-plastic geotechnical analysis due to its simplicity in application. In particular, some limitations create a barrier in its extensive application. Most of the cases, the intermediate principal stresses are neglected, and the MC criterion is expressed as a function of major (normal stress) and minor (shear stress) principal stress. However, in practical cases, the influence of intermediate principal stresses are predominant to determine the yield strength of the material [1, 9, 25, 37]. Another problem is that there are six sharp corners and tip which create singularity problems in the numerical implementation, particularly at axial symmetry conditions, whereas, the yield curve is in the form of an irregular hexagon on deviatoric (π) plane with fastigium and wracking [25, 50]. For that reason, to eliminate the computational difficulties, the apex and edges of the MC yield surface should be smoothed or rounded or special algorithmic treatment is needed [23]. Consequently, different types of techniques are proposed by numerous researchers to solve this problem for the smooth approximation of MC yield surface, for example, Drucker–Prager model, Mogi criterion, Twin Shear Unified criterion, SMP criterion, Zienkiewicz–Pande criterion, Williams–Warnke model, Gudehus–Argyris model, Rubin smooth function [4, 14, 27, 31,32,33, 37, 42, 58,59,60,61,62].

A numerical study by Walters and Thomas [56], considering the sandbox experiment, have shown that non-associated strain-softening and mobilized dilatancy are important parameters to determine the propagation of localized failure surfaces above trapdoors in the sand. In this regard, Sakai and Tanaka [43, 44] simulated the shear band proliferation precisely, in their two-dimensional (2D) anchor using a non-associated strain hardening–softening model with a shear band effect and nonlinear mobilized dilatancy in the elasto-plastic framework.

However, Meyerhof and Adams [36] proposed an empirical correction factor (“Shape factor”) to the peak uplift load for the rectangular anchor, calculated based on plane strain assumptions, as a function of embedment ratio only. Also, Dickin [11] and Frydman and Shaham [15] proposed the same as a function of embedment and aspect ratio. Besides, Murray and Geddes [38] and Vermeer and Sutjiadi [54] recommended simplified design equations and Ovesen and Strømann [40] proposed empirical design charts for isolated anchors.

In this study, a rigorous numerical model incorporating strain hardening–softening or strain-softening law with the shear-band effect is used to examine the implications of the implementation of different approximate models and rounding parameter to the 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. Air-dried uniform Toyoura sand (specific gravity (\({\text{G}}_{\text{s}}\)) = 2.64, effective diameter of particle (\({\text{D}}_{50}\)) = 0.016 cm, maximum void ratio (\({\text{e}}_{ \text{max} }\)) = 0.98, minimum void ratio (\({\text{e}}_{ \text{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/m3; medium having a density of 1480 kg/m3; loose having a density of 1350 kg/m3, error allowance = ± 10 kg/m3) 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.

Fig. 1
figure 1

Testing apparatus

Fig. 2
figure 2

Experimental failure mechanism and FE mesh for circular anchor: a Experimental failure mechanism, and b FE mesh (H/D = 2, D = 10 cm, and the number of elements = 2400)

Fig. 3
figure 3

Experimental failure mechanism and FE mesh in 3D for rectangular anchor: a Experimental failure mechanism, and b FE mesh (H/B = 2, L/B = 2, B = 50 mm, elements = 14,208)

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, Np (= δ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

$${{\uptheta }}= \frac{1}{3}\sin^{ - 1} \left( {\frac{{3\sqrt 3 {\text{J}}_{{3{\text{D}}}} }}{{2{\text{J}}_{{2{\text{D}}}}^{3/2} }}} \right)$$
(1)

where J3D = third deviatoric stress invariants and J2D = second deviatoric stress invariants.

Table 1 Drucker–Prager approximation to Mohr–Coulomb model
Fig. 4
figure 4

Rounding of MC yield surface in π-plane

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

$${\text{f}}\left( {{{\upsigma }},{{\upkappa }}} \right) = {{\upalpha I}}_{1} + \frac{{\sqrt {{\text{J}}_{2} } }}{{{\text{g}}\left( {{\uptheta }} \right)}} - {{\upgamma }} = 0$$
(2)

where \({{\upalpha }} = \frac{2\sin \upphi }{{\sqrt 3 \left( {3 - \sin \upphi } \right)}}\), \({{\upgamma }} = \frac{{6{\text{c}}\cos \upphi }}{{\sqrt 3 \left( {3 - \sin \upphi } \right)}}\) and the function, \({\text{g}}\left( {{\uptheta }} \right)\), is such that \({\text{g}}\left( {\frac{{{\uppi }}}{6}} \right) = 1\), like shape function, and this function determines the shape of the π-section. For instance, for MC, it becomes

$${\text{g}}\left( {{\uptheta }} \right) = \frac{{\frac{{\left[ {\left( {\cos \frac{{{\uppi }}}{6} - \sin \frac{{{\uppi }}}{6}} \right)\sin \upphi } \right]}}{\sqrt 3 }}}{{\left( {\cos {{\uptheta }} - \frac{{\sin {{\uptheta }}\sin \upphi }}{\sqrt 3 }} \right)}} = \frac{3 - \sin \upphi }{{2\sqrt 3 \left( {\cos {{\uptheta }} - \frac{{\sin {{\uptheta }}\sin \upphi }}{\sqrt 3 }} \right)}}$$
(3)

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 \({{\uptheta }} = \pm \frac{{{\uppi }}}{6}\) in the MC’s π-section as defined by Eq. (3), this function can be written as [2], graphic is shown in Fig. 5.

$${\text{g}}\left( {{\uptheta }} \right) = \frac{{2{\text{K}}}}{{\left( {1 + {\text{K}}} \right) - \left( {1 - {\text{K}}} \right)\sin 3{{\uptheta }}}}$$
(4)

where \({\text{g}}\left( { - \frac{{{\uppi }}}{6}} \right) = {\text{K}}\) and to match with MC failure criteria, the value of K can be calculated using Eq. (3) as

$${\text{K}} = \frac{{\sqrt {{{\text{J}}_2}} \left( {{\text{at }}\uptheta = - \frac{\uppi }{6}} \right)}}{{\sqrt {{{\text{J}}_2}} \left( {{\text{at }}\uptheta = + \frac{\uppi }{6}} \right)}} = \frac{{3 - \sin \upphi }}{{3 + \sin \upphi }}$$
(5)
Fig. 5
figure 5

Effect of different approximations in the π-plane

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

$${{\uptheta }} = \tan^{ - 1} \left( {\frac{\sin \upphi }{\sqrt 3 }} \right)$$
(6)

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 strength for soils. As the strains predicted by the MC model are conservative, this model does not capture the strain-softening behavior which is found in natural soils. Conversely, the DP model considers the effect of intermediate principal stresses. Figure 6 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.

Table 2 Parameters for hardening and softening material model
Fig. 6
figure 6

Approximation to Mohr–Coulomb material model using 2D FEM: a dense sand, b medium sand, and c loose sand

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 (\({{\uptheta }}_{\text{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.

Fig. 7
figure 7

Effect of rounding (H/D = 2, D = 10 cm) using 2D FEM: a dense sand, b medium sand, and c loose sand

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.

Effet of embedment ratio

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

Fig. 8
figure 8

The effect of embedment: a dense sand, b medium sand, and c loose sand

Fig. 9
figure 9

Plastic deformation zones above the anchor in dense sand (figures are not to scale): a H/D = 1, and b H/D = 6

Fig. 10
figure 10

Relationship between a Peak resistance factor, and b displacement factor as a function of H/D

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 parameters used for the analysis are summarized in 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 (\({\text{N}}_{\text{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.

Table 3 Parameters for simple strain softening material model
Fig. 11
figure 11

Uplift resistance–displacement factor relationships at H/B = 2 in dense sand

Fig. 12
figure 12

Shape effect on different types of sand

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 analyzing the influence of dilatancy. The uplift capacity of vertically uploaded circular anchor foundations with H/D varying from 1 to 6 has been simulated where \(\upphi\) ranges from 30° to 45° and ψ from 0° to 45°. Figures 13 and 14 show the values of \({\text{N}}_{\text{pu}}\) to peak frictional and dilation angle at H/D = 1 to 6, respectively. The uplift capacity increases significantly with the peak frictional and dilation angle and the effect of dilatancy becomes greater with the increase of frictional angle and embedment ratio. Greater uplift capacity is obtained for a greater dilation angle when the embedment ratio and frictional angle are kept constant. When H/D = 2, \(\upphi\) = 30° and 45°, the uplift capacity in the case of ψ = 30° is about 1.5 and 4 times the values in the case of ψ = 0°, respectively. When H/D = 6, \(\upphi\) = 30° and 45°, the uplift capacity in the case of ψ = 30° is about 1.7 and 8.8 times of the values in the case of ψ = 0°, respectively.

Fig. 13
figure 13

Effect of peak frictional angle on breakout factors for various H/D (D = 10 cm): a H/D = 1, b H/D = 2, c H/D = 3, d H/D = 4, e H/D = 5, and f H/D = 6

Fig. 14
figure 14

Effect of dilation angle on breakout factors for various H/D (D = 10 cm): a H/D = 1, b H/D = 2, c H/D = 3, d H/D = 4, e H/D = 5, and f H/D = 6

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

$${\text{S}} = \frac{{{\text{N}}_{\text{pu}} \left( {\text{Rectangular Anchor}} \right)}}{{{\text{N}}_{\text{pu}} \left( {\text{Continuous Anchor}} \right)}}$$
(7)
Fig. 15
figure 15

Peak resistance factor as a function of H/B (B = 10 cm) by different studies

Fig. 16
figure 16

Shape effect by different studies (H/B = 2, B = 5 cm)

The effect of shape diminishes with the increase of aspect ratio, gradually decreased to reach the strip anchor condition at \({\raise0.7ex\hbox{${\text{L}}$} \!\mathord{\left/ {\vphantom {{\text{L}} {\text{B}}}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\text{B}}$}} = \infty\). 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. It can be reported that the numerical model presents a very close agreement with past studies.

Fig. 17
figure 17

Effect of frictional angle by different studies at H/B = 2 (B = 10 cm)

Conclusions

The numerical models used in this research work demonstrate their appropriateness to closely predict the experimental uplift resistance–displacement relationships of vertical anchor problems in the sand. MC–DP model gives better results compared to the MC–MC, Argyris–Argyris, DP–DP (different matching) which proves its appropriateness to use it in this study. Conversely, it is required to verify the MC–DP model for other boundary value problems and soil type. The bearing capacity of the anchor foundation decreases with the increase of aspect ratio and increases with the increase of embedment. Furthermore, the settlement at the ultimate bearing capacity of the anchor foundations always increases with the increase of embedment and it cannot be well predicted by using a simple strain-softening model. Also, it is reported that the numerical model provides a very close agreement with past studies.