Axisymmetric adaptive upper-bound finite element limit analysis formulation based on second-order cone programming for bearing capacity of circular footing

TL;DR

This paper introduces a new computational method to accurately estimate the bearing capacity of circular footings using advanced finite element analysis.

Contribution

A novel adaptive axisymmetric upper-bound finite element limit analysis formulation is proposed using second-order cone programming.

Findings

The proposed method improves computational efficiency by using second-order cone programming.

Mesh adaptation based on plastic dissipation enhances accuracy with fewer elements.

Results show the method provides accurate upper-bound solutions for circular footing bearing capacity.

Abstract

An adaptive axisymmetric upper bound finite element limit analysis (UB-FELA) formulation has been presented for Mohr-Coulomb (M-C) materials in this paper. The computational domain is discretized using quadratic velocity elements. For the sake of computational efficiency, the axisymmetric UB finite element problem is recast into a model of second-order cone programming (SOCP). To enhance the precision of the proposed UB finite element method using a reduced element count, this study implements a mesh adaptation algorithm grounded in plastic dissipation. The collapse loads for determining the circular footings are then estimated by application of the proposed axisymmetric UB limit analysis formulas. By comparing the results to those reported in the literature, the analysis indicates that the method presented in this paper yields an accurate UB solution.

Figures

Quadratic velocity elements.

General layout and typical mesh of the problem (rough interface) for UB method.

The impact of the number of iterations on the load carrying capacity factor and the element count.

Mesh refinement in the calculation of Nc for smooth interface withφ = 30°. (a) fifth adaptive mesh refinement. (b) tenth adaptive mesh refinement. (c) fifteenth adaptive mesh refinement. (d) twentieth adaptive mesh refinement.

Mesh refinement in the calculation of Nc for rough interface withφ = 30°. (a) fifth adaptive mesh refinement. (b) tenth adaptive mesh refinement. (c) fifteenth adaptive mesh refinement. (d) twentieth adaptive mesh refinement.

Equations

Full text

1 Introduction

The efficacy of the finite element limit analysis (FELA) method in solving geotechnical stability problems has been validated by numerous researchers. According to different situations, the problems can be classified into plane strain problems, axisymmetric and 3D problems. In practical engineering, the circular excavation and bearing capacity of circular foundation can be considered as axisymmetric problems.

Pastor and Turgeman [1] initially proposed an axisymmetric UB-FELA method, grounded in linear optimization techniques, to examine this problem. This method need enforce a total of 3p inequality constraints at each node, where p represents the number of sides of the polygon employed for the linearization of the yield surface. Nevertheless, the methodology presented by Pastor and Turgeman [1] appears to be computationally less efficient. Kumar and Khatri [2] and Kumar and Chakraborty [3] proposed a novel linearization method for the yield criterion that reduces computational complexity by requiring each node to be subject to only three additional inequality constraints. However, the substitution of the nonlinear yield function with a multitude of linear constraints results in an increase in computational cost that is prohibitive for large-scale problems. Recently, SOCP has been utilised by numerous researchers as a means of addressing significant optimisation challenges [4–9]. The formulation with the usage of SOCP or Semidefinite programming (SDP) permits the treatment of the yield function in its intrinsic form, this approach addresses the limitation of the linearisation of the yield function inherent to the linear programming technique. SOCP has been used in 2D M-C problems [4–7], axisymmetric M-C problems [8–9], and 3D Drucker-Prager (D-P) problems [10], where has shown the high computational efficiency in FELA method. Tang et al. [8] were the initial researchers to apply SOCP to the axisymmetric LB problem for M-C materials, adopting the formulation proposed by Pastor and Turgeman [1].The computational results reported from Tang et al. [8] have demonstrated that SOCP can effectively solve the lower bound problems than linear programming. Consequently, the application of SOCP is rendered feasible for enhancing the computational efficiency of the previously mentioned limit analysis method for M-C materials. This article provides detailed expressions for the flow rule constraints and element plastic dissipation for the axisymmetric UB limit analysis method based on M-C materials. This article presents a comprehensive account of the expressions utilized in the flow rule constraints and element plastic dissipation within the FELA method. The implementation of the axisymmetric UB-FELA method encompasses these key research areas.

On the other hand, many research results have shown that the use of mesh adaptive refinement strategy and high-order element in the UB-FELA method is an effective means to improve computational accuracy[6,7,10–12]. Makrodimopoulos and Martin [6] provided an UB-FELA method with 6-node high-order triangle elements for M-C Materials. Sun et al. [12] provided an UB-FELA method based on the quadratic velocity fields and mesh adaptive refinement strategy. However, the utilization of mesh adaptive refinement strategy and high-order element on the axisymmetric UB method still needs further research.

In this manuscript, a detailed theoretical basis for an axisymmetric upper limit finite element calculation method based on the SOCP technique and quadratic velocity element are presented. Additionally, A mesh adaptive algorithm driven by plastic dissipation is also proposed with the aim of obtaining highly accurate UB solutions. Numerous studies have confirmed the ability of the mesh adaptive strategy to improve the accuracy of plastic limit analysis calculations [5]. The study also reveals that the proposed method yields accurate UB solutions, as evidenced by a comparison with previous research findings on the stability of circular foundations involving both smooth and rough substrate-soil interfaces.

2 Axisymmetric UB-FELA with SOCP

2.1 Quadratic velocity finite element

To discretize domain of computation, quadratic velocity elements are utilized. It has been demonstrated that quadratic velocity elements (6-node triangular element) can provide much better accuracy solutions but do not reduce the computation efficiency of solving UB optimization model [6]. Fig 1 shows the type of quadratic velocity element without velocities discontinuities between the adjacent elements. It is possible to express the velocity tensor at any specified point within a triangular element in terms of the velocity tensors at its six vertices.

Quadratic velocity elements.

$$[eqn]$$

In this context, ui and vi denote the velocities at a specified point in the radial (r) and axial (z) directions, respectively.

The term Ni represents the shape function associated with the quadratic velocity element and can be calculated using the following expressions.

$$[eqn]$$

Where L1, L2, L3 is the shape function of the 3-nodes triangular element,

$$[eqn]$$

where,

$$[eqn]$$

And $[eqn]$ is the area of the element.

The necessary element strains for the axisymmetric problem can be determined through,

$$[eqn]$$

In this context, the symbols εr, εz, and εθ are used to denote the plastic normal strain rates in the radial (r), axial (z), and circumferential (θ) directions, respectively. The symbol γrz is used to represent the plastic shear strain rate.

2.2 M-C Criterion for axisymmetric problem

The M-C Criterion in the context of an axisymmetric condition is formulated as follows [1,8].

$$[eqn]$$

The preceding three criteria can also be expressed in the form of three distinct types of 3D second-order cones.

$$[eqn]$$

Where $[eqn]$ , $[eqn]$ is nodal stresses, $[eqn]$ represents second-order cone constraint.

$$[eqn]$$ $$[eqn]$$

2.3 The plastic strain dissipation function and associated flow rule constraints

This section comprehensively introduces the application of SOCP within the axisymmetric UB-FELA methodology for M-C materials. The initial requirement for employing the UB limit analysis method is to derive an expression for the dissipation function, as well as the delineation of the set of plastically admissible strains[13]. Under axisymmetric conditions, the M-C criterion can be reformulated as a SOCP, leveraging the duality principles inherent in conic optimization based on Eq. (6).

In fact, it is difficult to directly establish the associated flow rule constraints without stress components based on axisymmetric M-C Criterion conditions. According to the research results of Makrodimopoulos [14], the analytical expressions for the associated flow rule constraints and plastic dissipation can be obtained by solving Eq. (9a).

$$[eqn]$$

where $[eqn]$ , $[eqn]$ . By solving Eq. (9a), the analytical expressions for the associated flow rule constraints and plastic dissipation shown in Eq. (9b) can be obtained.

$$[eqn]$$

where $[eqn]$ can be expressed by Eq.(10).

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Utilizing this method, According to Eq. (8), Eq. (9b), and Eq. (10), the expression for the dissipation function is reformulated as shown in Eq. (11).

$$[eqn]$$

Additionally, the associated flow rule constraints for each element can be obtained as follows,

$$[eqn]$$

Where $[eqn]$ are the auxiliary variables, $[eqn]$ are strain variables.

3 Numerical formulation

3.1 Calculation of power dissipation.

To mitigate the likelihood of errors during the computational process, the current section initially provides an analytical expression that describes the rate of external work related to body forces and the dissipation of power resulting from plastic deformation within the framework of continuous quadratic velocity fields. The further details regarding the derivation of the formula are presented in Appendix. In an axisymmetric problem solved by means of SOCP optimization, the power dissipated due to plastic deformation can be expressed through the subsequent formula:

$$[eqn]$$

where the rate of external work performed by body forces is given by the following analytical expression:

$$[eqn]$$

In this expression, A signifies the area of the element, γ denotes the unit weight of the soil mass, vi corresponds to the vertical velocity at the node, and r1, r2, r3 represent the radial coordinates for the three vertices of the triangular element.

3.2 Numerical implementation.

For a six-node triangular element, the plastic strain rate field (εr, εz, εrz) exhibits linear variation, whereas εθ remains a constant, corresponding to the strain at the centroid of the triangular element. The quadratic velocity fields is applied in the UB limit analysis method to help obtain high-precision UB solution. The flow rule constraint for axisymmetric UB formulation can be applied in the three vertices. This paper formulates the axisymmetric UB problem, incorporating finite elements, as a conventional SOCP problem. In fact, the large scale mathematical optimization calculation model based on second-order cone programming technology has been widely proven to be efficiently solved[15–20].Constructing an upper bound finite element limit analysis calculation model that satisfies the associated flow rule and velocity boundary conditions everywhere can ensure obtaining strict upper bound solutions. The axisymmetric upper bound finite element method proposed in this article assumes that the circumferential strain εθ inside each element is a constant. This assumption means that Eq. (15) does not fully comply with the associated flow rule, and this present method cannot obtain a strict upper bound solution. However, based on existing research [3,11] and the results of this study, this assumption does not lead to computational errors.

$$[eqn]$$

In which NE denotes the element number, NP represents the flow rule point number, NB is the boundary point number, and NP = 3NE. The strain components ε_r_, ε_z_, γ_rz_, ε_θ_ can be converted into the expression of velocity components based on Eq.(1) and Eq.(5). ω represents the angle between the velocity boundary and the horizontal coordinate axis. The analytical expression of flow rule constraints and plastic dissipation is of great significance for the constructing upper bound finite element limit analysis calculation programs. The upper bound finite element limit analysis model based on second-order cone programming can be obtained through Eq. (15), which is only related to node velocity variables and second-order cone optimization variables.

Moreover, the computational program for the UB finite element method, founded on SOCP and utilizing continuous secondary velocity field, is coded in MATLAB. All computations were executed on a Lenovo system equipped with a 3.2 GHz CPU and 8GB RAM, operating under the Windows 7 operating system.

4 Mesh adaptive strategy

Numerous scholars have adapted the adaptive meshing technique for application within the FELA method [11,21]. Nonetheless, research on the application of adaptive meshing techniques to axisymmetric problems remains nascent. Zhang et al. [21] introduced an adaptive meshing strategy based on plastic dissipation for axisymmetric UB limit analysis problems. This paper presents a plastic-dissipation-based adaptive strategy to enhance the precision of determining UB solutions and associated failure mechanisms. The proposed plastic-dissipation-based adaptive meshing strategy is analogous to that introduced by Zhang et al. [21], employing a 3-node triangular element. The mesh adjacent to failure zones should be refined, whereas the mesh located farther from plasticized zones may be coarsened incrementally. It appears that elements adjacent to the failure regions are likely to exhibit high plastic-dissipation Pe values. The plastic-dissipation Pe values for each element are sorted in descending order. Consequently, elements requiring refinement can be identified using the local refinement indicator η.

$$[eqn]$$

In which ne denotes the total count of marked elements, NE represents the total count of all elements. An appropriate value for η can be determined through trial calculations, typically lying within the range of 0.4 to 0.5. A value of η equal to 0 signifies that no elements should be refined, whereas η equal to 1 indicates that all elements within the computational domain are to undergo a refinement process. The adaptive limit analysis program is implemented following the methodology outlined below.

Address the UB problem as a SOCP optimized method, which based on using the initial set of elements.Initially, the elements are sorted in descending order based on the computed plastic damage value Pe, followed by the identification of elements for refinement based on the refinement factor η. The plastic-dissipation Pe values for each element can be obtained based on Eq. (21), in which the necessary data can be extracted from the solution results.Two new subtriangles are generated from each parent triangle by connecting the vertices of the parent triangle to the midpoint of its longest side.Based on the refined mesh data, rebuild the upper bound finite element limit analysis calculation model. Return to step 1 and start a new calculation.Each adaptive step is iterated through the predetermined number of iterations to achieve the optimal UB solution and the associated failure mechanism.

5 Numerical examples

The bearing capacity of circular foundations is a classic axisymmetric stability problem[21–26]. A procedure was conducted to compute the bearing capacity of circular foundations subjected to both smooth and rough substrate-soil interfaces, and the accuracy and efficacy of this procedure were evaluated based on the resultant outcomes. The accuracy, computational efficiency, and failure mechanisms associated with the structures were assessed in the context of the UB solution.

Only half of the total construction shown in Fig.2 is considered and $[eqn]$ is diameter of circular footing. The soil mass is characterized as a M-C material, possessing soil unit weight (γ), soil cohesion (c), and internal friction angle (φ).

General layout and typical mesh of the problem (rough interface) for UB method.

The bearing capacity for a rigid circular footing devoid of surcharge is expressed in [8] as:

$$[eqn]$$

Wherein Nc, which is a function of the internal friction angle φ, denotes the coefficient of load-bearing capacity. Qu is the ultimate bearing capacity of circular foundation, which can be obtained by the upper bound model shown as Eq. (15). The perfectly rough footing-soil interface condition can be simulated by specifying the node velocity in the u and v directions for the footing nodes, u = 0, v = -1. For a smooth footing, the interface condition can be simulated by specifying the node velocity in the v direction for the footing nodes, v = -1.

5.1 Method validation

Table 1 presents the bearing capacity coefficients Nc for both smooth (bracketed) and rough (unbracketed) substrate-soil interfaces, with friction angles φ ranging from 5° to 30°. The results are juxtaposed with (b) the UB solutions of Mohapatra and Kumar [11], which employ discontinuous quadratic velocity fields, (c) the UB solutions of Kumar and Chakraborty [3], utilizing a three-node element with velocity discontinuity, (d) the UB solutions of Zhang et al. [21], which are based on a three-node element with velocity discontinuity and employ a mesh adaptive strategy, (e) the LB solutions of Kumar and Khatri [2], obtained using a three-node element, (f) solutions derived from the stress characteristics method [23], (g) solutions derived from the stress characteristics method [22], (h) numerical solutions based on the FLAC 2D method [24]. Table 1 indicates that the results derived from the current analysis are in reasonable agreement with the various solutions documented in the literature. It is significant to observe that the bearing capacity values for the smooth footing-soil interface are positioned between the UB and LB solutions as documented in the literature. This implies that the current method is capable of yielding UB solutions with satisfactory accuracy. Additionally, it is observed that the UB solution obtained from the present method is marginally higher than that derived from the stress characteristics method.

Table 1: Comparison between the Nc values of the circular footings.

5.2 Discussion on the effectiveness of mesh adaptive strategy

The computation accuracy can be also improved by using more elements, but it seems to be a large time consumption. It has been proved that an adaptive refinement method can help obtain high-precision results with fewer elements in axisymmetric UB problems [21]. A plastic dissipation-driven mesh adaptation strategy, akin to the approach presented in literature [21], is employed in this paper to enhance the UB solutions.

Fig 3 indicates that the adaptive strategy adopted by this paper can improve the calculation accuracy of UB solution. However, note that the element number increases gradually with the adaptive number increases, but the calculation accuracy is limitedly improved after a few times mesh refinement. It means that if the numeric results are approximate to actual value, improving the computational results by using mesh refinement becomes very difficult.

The impact of the number of iterations on the load carrying capacity factor and the element count.

Fig 4 and Fig 5 illustrates the alterations in the refined mesh and the bearing capacity coefficient Nc for both smooth and rough foundation-soil interfaces at a friction angle of φ = 30°. It is important to note that only the local region is depicted to facilitate a clearer explanation of the issue. The area of damage associated with the rough footing is markedly greater than that associated with the smooth footing, and a distinct disparity exists between the damage modes of the two footings. The observed failure pattern is in close agreement with the characteristic meshes presented by Martin [23]. Fig 4 and Fig 5 further indicates that the arrangement of the refined mesh parallels the patterns of power dissipation, consequently, the refined mesh configuration proposed herein may indirectly reveal the failure pattern. Additionally, Fig 4 and Fig 5 demonstrates that the meshes located distant from the failure zone are significantly coarser compared to those within the failure zone, this inference suggests why the adaptive refinement strategy is effective in achieving a superior UB solution with a reduced number of elements.

Mesh refinement in the calculation of Nc for smooth interface withφ = 30°. (a) fifth adaptive mesh refinement. (b) tenth adaptive mesh refinement. (c) fifteenth adaptive mesh refinement. (d) twentieth adaptive mesh refinement.

Mesh refinement in the calculation of Nc for rough interface withφ = 30°. (a) fifth adaptive mesh refinement. (b) tenth adaptive mesh refinement. (c) fifteenth adaptive mesh refinement. (d) twentieth adaptive mesh refinement.

5.3 Detailed Nc values for the rigid circular footing

The precise UB and LB solutions for the bearing capacity coefficients Nc of smooth (bracketed) and rough (unbracketed) substrate-soil interfaces are presented in Table 2, with the friction angle φ ranging from 0° to 39°. The UB solutions listed in Table 2 are derived from the presented UB method, whereas the LB solutions are obtained through the LB method proposed by Sun and Yang [27]. It should be noted that the UB and LB solutions for bearing capacity, corresponding to smooth and rough footing-soil interfaces, are closely aligned. This implies that the current method is capable of yielding UB solutions with satisfactory accuracy.

Table 2: The UB and LB of the bearing capacity factors Nc.

6 Conclusions

An upper bound limit analysis (UB-LA) formulation has been proposed to address axisymmetric geomechanically stability challenges, integrating finite element method and SOCP optimization techniques. This study also incorporates a mesh adaptive refinement strategy and employs high-order elements to facilitate the acquisition of high-precision UB solutions with a reduced number of elements. The results of the current analysis are demonstrated to be in close alignment with the UB and LB solutions documented in the literature. Consequently, the developed method is employed to examine the bearing capacity of a circular foundation in relation to both smooth and rough footing-soil interfaces. The method holds considerable promise as a technique for deriving highly precise UB solutions and delineating distinct failure patterns in axisymmetric geomechanically scenarios.

Appendix

The rate at which external work is performed by body forces is defined as,

$$[eqn]$$

According to the finite element shape function, the Jacobian function can be expressed as,

$$[eqn]$$

The power dissipation $[eqn]$ can be calculated by,

$$[eqn]$$

The power dissipation due to the plastic deformation with continuous quadratic velocity field is defined as

$$[eqn]$$

The same approach can be used to deal with the integral.

$$[eqn]$$

Where,

$$[eqn]$$

Supporting information

S1 FigQuadratic velocity elements.(TIF)

S2 FigGeneral layout and typical mesh of the problem (rough interface) for UB method.(TIF)

S3 FigThe impact of the number of iterations on the load carrying capacity factor and the element count.(TIF)

S4 FigMesh refinement in the calculation of N_c_ for smooth interface with φ = 30°.(TIF)

S5 FigMesh refinement in the calculation of N_c_ for rough interface with φ = 30°.(TIF)