Finite element implementation and application of a sand model in micropolar theory

In traditional finite element failure analyses of geotechnical structures, the micro grain rotations cannot be modelled and numerical solutions are mesh dependent. In this study, a user element including rotational degree of freedom has been developed based on micropolar theory (Cosserat theory), then an enhanced non-associated sand model is calibrated with laboratory data and used to model the plane strain tests. The simulated results demonstrate the polarized model is able to model reasonably the sand behavior as well as the grain rotations in the localized region. What’s more, with this enhanced model, the mesh independent numerical solutions in terms of mechanical responses, shear bands thickness and orientations have been obtained. In failure analysis of geostructures, significant rotations of soil grains have been observed to occur in the strain localized regions, but the current commercial Finite Element tools cannot model the micro rotations. Therefore, a user defined element must be developed to include the rotational degree of freedom. The micropolar approach is proven to be effective to model the grain rations in present paper. More suitable than other classical soil or sand constitutive models, the selected non-associated sand model in present paper is capable of describing well the contraction and shear dilatancy behaviors of sand. Then the model has been enhanced by means of micropolar technique, in this way, the reasonable strain localization phenomena in laboratory tests could be predicted well. There are always the mesh dependent problems for traditional simulations of the strain localization phenomena in finite element analysis. It can be found in present paper that the mesh independent numerical solutions are obtained by means of micropolar technique. Furthermore, the micropolar approach can obviously improve convergence difficulties in finite element analyses. In failure analysis of geostructures, significant rotations of soil grains have been observed to occur in the strain localized regions, but the current commercial Finite Element tools cannot model the micro rotations. Therefore, a user defined element must be developed to include the rotational degree of freedom. The micropolar approach is proven to be effective to model the grain rations in present paper. More suitable than other classical soil or sand constitutive models, the selected non-associated sand model in present paper is capable of describing well the contraction and shear dilatancy behaviors of sand. Then the model has been enhanced by means of micropolar technique, in this way, the reasonable strain localization phenomena in laboratory tests could be predicted well. There are always the mesh dependent problems for traditional simulations of the strain localization phenomena in finite element analysis. It can be found in present paper that the mesh independent numerical solutions are obtained by means of micropolar technique. Furthermore, the micropolar approach can obviously improve convergence difficulties in finite element analyses.


Article highlights
(1) In failure analysis of geostructures, significant rotations of soil grains have been observed to occur in the strain localized regions, but the current commercial Finite Element tools cannot model the micro rotations. Therefore, a user defined element must be developed to include the rotational degree of freedom.The micropolar approach is proven to be effective to model the grain rations in present paper. (2) More suitable than other classical soil or sand constitutive models, the selected non-associated sand model inpresent paper is capable of describing well the contraction and shear dilatancy behaviors of sand. Then the modelhas been enhanced by means of micropolar technique, in this way, the reasonable strain localization phenomenain laboratory tests could be predicted well. (3) There are always the mesh dependent problems for traditional simulations of the strain localization phenomena in finite element analysis. It can be found in present paper that the mesh independent numerical solutions are obtained by means of micropolar technique. Furthermore, the micropolar approach can obviously improve convergence difficulties in finite element analyses.

Introduction
For the failures of geotechnical structures related to granular assemblies, strain localization always occurs accompanied with the reduction of bear capacity. In this process, significant grain rotations inside the localized regions play important roles [1]. Moreover, researchers have confirmed grains' translations and rotations in shear strain localized region with the digital correction of particle-scale volume, allowing the complete kinematics of particles to be characterized in laboratory tests [2]. Many constitutive models have been developed to describe the behavior of granular soils. However, most models are formulated based on classical continuum mechanics theory in finite element analysis, which cannot consider rotational degrees of freedom (dofs), but can only take into account translational dofs. Consequently, when modeling the strain localization phenomena, there are mesh dependent problems and pathological solutions. One of the potential reasons for this defect is that when the simulated localized region is over-developed for static problems using classical models, the property of field control equations will change from elliptical to hyperbolic. Consequently, we cannot obtain the unique numerical solution, especially in the post-failure stage.
In recent years, various regularization techniques, e.g. micropolar theory, viscosity, non-local theory, and highgradient theory, have been used to alleviate the mesh dependency or sensitive problems in simulating strain localizations by finite element method (FEM). Compared with other approaches, micropolar theory is easy to be implemented and physically meaningful. Moreover, it is also suitable for capturing shear-dominated failure mode. Accordingly, extending a classical model to the micropolar model has been embraced a lot. Numerical simulations using micropolar Lade hardening model have been carried out by Alsaleh et al. [3,4] to model the strain localization phenomena. Arslan and Sture [5] used the micropolar Drucker-Prager model to study the micro and macro relations in materials. de Borst and Sluys [6] investigated the strain localization phenomena and strain softening behaviors under both static and dynamic condtions with the micropolar J 2 model. The micropolar approach was also used to model the behaviors of hypoplastic materials by Huang and Bauer [7]. Li and Tang [8] used a pressure dependent micropolar model to simulate the strain localizations of a soil layer and the failure of a slope. Different from other researchers' attempts on conventional models, such as J 2 model, Drucker-Prager model among others, a non-associated critical state based model for granular materials has been first formulated within the micropolar framework in present study. The polarized sand model demonstrates great ability to reasonably model the micro grain rations and the localized failure mode in plane strain tests. Furthermore, the improved model is able to guarantee the accuracy of solutions and efficiently overcome the mesh dependency problems in numerical simulations. In other words, the polarized model is more reasonable than the model based on classical theory to model shear bands in laboratory tests, in terms of mechanical responses, grain rotations, shear band patterns, etc.
In this paper, the formulations of elastic and plastic models based on micropolar theory have been illustrated and derived in detail at first. Then the enhanced non-associated sand model by the micropolar approach has been implemented into a user subroutine with finite element method. After that, the polarized model has been calibrated and verified by the laboratory plane strain tests. At last, the applications of the polarized sand model has been conducted by modelling shear bands in different meshes.

General formulations of model in micropolar theory
In addition to the translational field, extra rotational dofs are endowed to the element in micropolar theory (also called Cosserat theory [9]). Rotational dofs, independent from the translational dofs, are easily added in the balance equations at constitutive level. Furthermore, micro rotations and the corresponding couple stresses can be explained from the viewpoint of material mechanics. The rotations in a micropolar model are considered as the rotations of grains or aggregates. In contrast, the classical model considers only the translations or displacements of material points. In order to illustrate the micropolar theory clearly, the plain strain problems are first investigated. After adding a rotational dof to the micro-element, each element has three generalized dofs.
Contrary to the stress and strain in traditional continuum mechanics, the strain and stress vectors in micropolar media can also consider microscopic curvatures and moments, in which m zx and m zy are couple stresses corresponding to the micro curvatures κ zx and κd zy , and l c , is a internal size scale like mean particle size d 50 .
For a model based on micropolar theory, not only the equilibrium and compatibility equations, but also the constitutive laws have been generalized by including the additional components of strain and stress vectors [10].

Equilibrium equations
The unified equilibrium governing equations for both static and dynamic problems in micropolar theory can be expressed in which σ ij and m kj,j denote Cauchy stress and microcouple stress, and f j and c k represent body force and micro moment, respectively. The first sub-equation is the same with that in a classical model. The second one is a special form in micropolar theory. For the shear stress components in micropolar theory, the symmetrical part S = xy + yx 2 produces shear strain whereas the skew part A = xy − yx 2 only causes rotations.

Kinematics equations
In micropolar theory, the formulations of normal strain components xx and yy are identical to the forms in classical continuum theory, whereas, the shear strain components xy and yx are linked not only to the translations but also to the rotation. The curvatures xz and yz are calculated from the gradients of rotational quantity ω z .

Elastic constitutive laws
Summarize the equilibrium and kinematic equations in matrix-vector forms, where u (including rotation) and f (including body moment) denote displacement vector and body force vector, respectively, and L is named strain operator matrix, (4) ij,i + f j = 0 m kj,j + e kij ij + c k = 0 The strain rate is related to stress rate through elastic stiffness matrix, in which λ = 2 Gυ/(1-2υ) υ is Poisson's ratio, G and G c denote shear modulus and micropolar shear modulus, respectively.

Elastoplastic constitutive laws
For small strain micropolar continuum, it is supposed that total strain includes elastic part and plastic part. Substitute Eqs. (12) into (10), we get The plastic strain is formulated in which ̇ is named plastic multiplier, m denotes the direction of plastic flow and Q is the potential function. The consistency condition controls plastic flow. In plasticity theory, the Kuhn-Tucker must be guaranteed where κ contains hardening variable vector, and = and h = −̇̇ . Then, the stress-strain relation in a micropolar continuum can be obtained by combining Eqs. (13) and (17) 3 Extension of classical model to a micropolar one

Introduction of the classical model
The models for granular materials have been improved from linear elasticity to ideal plasticity and then to nonlinear hardening during the past decades. In present study, a different model based on critical state proposed by Yin and Hicher [11] has been selected for polarization. The constitutive laws in p'-q plane are derived The forms of bulk modulus K and shear modulus G are expressed [12] in which e denotes void ratio, υ is Poisson's ratio, p at and p' represent the atmospheric pressure and mean effective stress, respectively. ζ controls the nonlinear effect. K 0 , G 0 are the initial values of K and G, which can be linked by υ,

The plastic strain is formulated
The yield function is expressed in which k p , controlling the initial slope of the curve q∕p � − p d , affects the plastic shear modulus, and M p = 6sin(ϕ p )/(3-sin(ϕ p ), representing the peak stress ratio, can be obtained from the peak friction angle ϕ p . The potential surface, considering dilation and contraction, can be formulated as where A d is the input parameter controlling dilatancy, M pt = 6sin(ϕ pt )/(3-sin(ϕ pt )) represents the stress ratio at phase transformation state in p'-q plane. ϕ p and phase transformation friction angle ϕ pt be calculated by ϕ u and e c (on the critical state line in p'-q plane and in e-log p' plane in Fig. 1).
And we calculate void ratio e c on the critical state line Input parameters n p and n d control the peak and phase transition point, λ represents the CSL slope in e-log p' plane, and e ref means the reference critical void ratio on critical state line at p at . From above equations and the p'-q-e plane Fig. 1, we can find that for both loose and dense materials, it guarantees the void ratio and the stress reach the critical state line simultaneously.
The plasticity multiplier dλ can be derived as the manner in conventional plasticity theory, then the constitutive laws can be solved.

Enhanced model in micropolar theory
By considering curvatures and coupled stresses, the classical model aforementioned has been extended to the micropolar model. Thus, more advanced than the strain and stress invariants of a classical model, the modified invariants of the micropolar model were formulated.
In the formula, the plastic micro-curvature rate tensor and the plastic deviatoric strain rate tensor are represented by ̇p ij and ̇e p ij , and m ij and s ij denote the micro moment tensor the deviatoric stress tensor, respectively. Thereafter, the deviatoric stress q can be expressed in a matrix-vector form where P matrix is defined Similarly, the equivalent plastic strain p d is newly updated in matrix-vector form with the definition of matrix Q It can be observed that after the strain and stress vectors are expanded to include micro curvatures and micro moments, the strain and stress invariants must be also redefined, then the traditional model can be extended to a micropolar model. If the rotational dof is restricted, the micropolar model can be also retrieved to a traditional model.

User defined element in micropolar theory
As aforementioned, a micropolar element should have three dofs for plane strain problem. Through the user interface in software ABAQUS, a user element with 8 nodes and 4 reduced integration points, considering the independent rotation has been developed as shown in Fig. 2. According to FEM, the displacement and rotation u of the user element can be derived from the counterparts δ e at each node via the interpolation function N.
The interpolation functions are defined The global coordinates of the element can also be obtained by interpolating nodes' coordinates.
After the partial differential operation, we get The Jacobian matrix [J], used to realize the map between (ξ, η) and (x, y), is formulated

Finite element discretization
According to the virtual displacement principle for plane strain problem, the potential energy of a structure can be expressed as the weak form in which f and T denote the interior body force vector and the external surface force vector, respectively. t represents the thickness. Substitute Eqs. (40) into (44), and the total structural system potential energy can be obtained by adding all the discretized elements' potential energy.
For a random virtual displacement, it must be required the partial differential Π p = 0 . In this way, the equilibrium equations of an element can be formulated K e and P e represent the stiffness matrix of the element and the load vector acting on the nodes, respectively. Then the governing field equations for a structure can be expressed in which a denotes the array containing all the node displacements, and K and P represents the structure stiffness matrix and the equivalent node load array of the total system, respectively. For static problems, the Newton-Raphson iteration is adopted to conduct the equilibrium, and the state related variables can be calculated and updated by implicit integration in ABAQUS. In finite element analysis, random discretized irregular elements can be transformed to regular ones via Jacobian matrix. Then, the integration operation can be replaced by sum operation at the integration point. Taking the element stiffness matrix as an example

Calibration of model parameters
The polarized model includes 12 parameters, which are classified into four groups: (a) parameters about elasticity ζ, υ, and K 0 , (b) parameters defining critical state λ, ϕ u and e ref , (c) parameters interlocking plasticity n p , n d , A d and k p , (d) and the particular parameters l c and G c in the micropolar model.
The elastic parameters ζ, K 0 , can be obtained from an isotropic compression test, and the parameters interlocking plasticity n p , n d , A d and k p , can be identified with one To determine the critical state parameters λ, ϕ u as well as e ref , mechanical curves of no less than three set of triaxial tests with different initial densities and confining pressures need to be fitted. Poisson's ratio υ can be assigned based on materials.
For the determination of the micropolar parameter l c and G c , it has been found that the shear band thickness is linearly proportional to d 50 in laboratory biaxial tests or linearly proportional to the micro scale l c in numerical simulations [3,[13][14][15][16]. Therefore, l c , representing the scale of microstructure, is set to be d 50 , and G c is set to half the shear modulus G as others do [17,18].
Considering the low efficiency of conventional curvefitting method, The genetic optimization method [19], finding model parameters more quickly and accurately, has been adopted to conduct the parameters calibration.
Laboratory data of F-75 Ottawa sand shown in Fig  selected to identify the model parameters. The experimental data shown in Figs. 5 and 6 of one isotropic compression test [20] as well as five triaxial tests [21] has been fitted. In experimental tests, with a constant confining pressure applied to the specimens using a cell pressure reservoir, a constant axial displacement rate is applied on the specimens. The bottom end platen is restrained from movement and the top end platen is rigidly connected to the loading ram. During the testing process, the specimen's deformation is monitored by noting the displacements of the grid imprinted on the membrane surface covering the specimen. The void ratio can be obtained according to the volume of drained water. The elastic parameters obtained from the isotropic compression test are supposed to be universal for triaxial tests. There are 9 model parameters should be identified by fitting the laboratory data, and the rest 3 parameters (υ, l c and G c ) can be set in advance. It requires at least three drained triaxial tests with different initial confining pressure and void ratio to inversely search for the critical state parameters using the genetic optimization approach. The other two parallel triaxial tests can be used to verify the calibrated parameter. The calibrated parameters are listed in Table 1.
From the comparisons between laboratory data and simulations in Figs. 5 and 6, it can be found that the behavior of F-75 dry sand in isotropic compression and triaxial tests can be very well predicted by the micropolar model.

Experimental verification by biaxial tests
In order to study the effects of particle size, void ratio, and confining pressure on the performances of granular soils, researchers have carried out a lot of biaxial tests on F-75 Ottawa sand [22]. In present study, four experimental data set of biaxial tests of dense and medium dense sand under high and low confining pressure, have been selected to conduct model verification (i.e. D r = 55%, σ c = 15 kPa; D r = 47%, σ c = 100 kPa; D r = 97%, σ c = 15 kPa; D r = 87%, σ c = 100 kPa). During testing process, the bottom and top end platens were constrained from lateral movement. With the confining pressure imposed on testing specimen, a continuous displacement increment was applied on the top of the specimen. And deformation configuration was finally recorded on the membrane grid covering on the specimen.
When modeling strain localization, the numerical simulations by the classical model were proven to be pathological and mesh dependent in post-bifurcation regime. Therefore, the enhanced micropolar model was adopted herein to describe the granular soil behavior. l c was regarded as  Table 1. Notably, the critical friction angle ϕ u has been slightly increased as shown in Table 2. The experimental findings of Alshibli et al. [21] can support the improvement of ϕ u , that the strength of specimen in plane strain test is slightly higher than that in triaxial test. The restrained lateral movement of the platen end may lead to a higher residual stress. What's more, the difference of ϕ u may be also explained by the differences of Lode angles in different cases.
In numerical analysis, a weak element was assigned to top left corner of the model to induce a single shear band as the shear band patterns in the laboratory tests. Because obvious grain rotations have been observed inside shear strain localized zones, accordingly, the simulated shear bands characterized by micro rotations were compared with the shear bands in laboratory specimen. Figure 7, UR3 representing the rotations, displays the comparison between the shear band in a laboratory specimen and a numerical model of dense sand with high confining pressure, which shows that the shear band in a testing specimen can be well modelled by the micropolar model, in terms of its location, inclination, and thickness. Figures 8   and 9 plot the simulated and experimental data of F-75 sand in biaxial tests in the same plane, from which it can be found that the principal stress ratio vs. axial strain curves of experimental data are modelled well and the volumetric strains can also be moderately predicted.

Application in mesh dependency problems
Admittedly, the numerical solutions to strain localization by the classical model may have numerical difficulties or suffer from ill-posedness in finite element analysis. As a regularization approach, the performances of the polarized model in improving convergences and dealing with mesh dependencies have been investigated. Numerical simulations of plane strain test, with model dimensions of 20 cm height, 10 cm width, and 1 m thickness, have been conducted by the classical model as well as the micropolar model. Lateral movements of the bottom and top ends of the model were locked. There were two stages in the simulations. The first stage was the 100 kPa confining pressure applied on the model, and

Shear band and mechanical behavior
Firstly, the numerical simulations of shear bands were performed by the classical model. The distributions of plastic strain of four different meshes are displayed in Fig. 10, and four corresponding curves of displacement vs. load are plotted in Fig. 10. From the shear bands in Fig. 10 and the mechanical responses in Fig. 10, it can be found that the simulations with mesh 10 × 20 and mesh 15 × 30 have been completed, whereas, there are convergence difficulties around the peak point for mesh 20 × 40 and mesh 30 × 60. The bifurcations of numerical calculations in modeling strain softening lead to solutions mesh dependent, resulting in shear bands concentrate in the domain of a single element. When the element size is extremely refined to zero, the consumed energy by failure also converges to zero, which is not the feature of real materials [23]. What's more, the acoustic tensors of many Gauss points in localized regions become singular, which is responsible for the convergence difficulties  Figures 10 and 11 show that numerical simulations are severely sensitive to mesh sizes, that is to say, the coarser the mesh is, the thicker the shear band will be, and the peak load capacity of coarse mesh is a little higher and delayed than fine mesh. Furthermore, it can be found the residual strength of a model with coarse mesh is greater than that with fine mesh.

Inclination of shear band by the classical model
When it came to the investigations of the mesh dependency problems in the past, many discussions were focused on shear band thickness and the strength in softening regime, however, the differences of shear band inclinations  have never been studied. In present study, we also considered the shear band inclinations of different meshes. The inclined angle between the horizontal line and the centerline of the shear band can be measured to define shear band orientations. The centerline may be regarded as the sliding surface. Because the computation of mesh 30 × 60 faces numerical difficulties earlier, shear band inclined angles of other three meshes are measured in Fig. 12 (i.e., β 1 = 52.69°, β 2 = 57.65°, and β 3 = 60.15°), from which we can find a finer mesh produces a steeper shear band inclination.
To summarize the numerical mesh dependent problems by the classical model. First, upon the occurrence of bifurcation, the acoustic tensors of the elements inside shear band will become singular, resulting in numerical problems and convergence difficulties. Second, shear bands patterns, specially the thickness and inclination, are severely sensitive to the element size. Third, load peak as well as the residual strength also significantly rely on the element size.

Shear band and mechanical behavior
To verify the performances of micropolar technique in solving mesh dependency problems, the shear bands in biaxial tests were simulated again using the micropolar model. The shear bands of four different meshes are displayed in Fig. 13. Better than the calculations by the classical model, all four simulations by the polarized one could be completely finished without convergence difficulties or numerical problems. And the very similar shear bands of the four different meshes can be easily observed at first glimpse. Furthermore, the curves of displacement vs. load of different meshes are presented in Fig. 14, showing that the mesh dependency problems have been significantly relieved even in the softening regime. Except the too coarse mesh 10 × 20, the mechanical curves of other three meshes are completely consist with each other.

Inclination of shear band by the micropolar model
The regularization efficiency of the micropolar approach in relieving mesh dependency problems can be also demonstrated by the the shear band inclinations with different meshes in Fig. 15, which shows a well unified shear band inclinations of the three fine meshes β 2 = β 3 = β 4 = 53.22° except for a slightly different β 1 = 53.10° in coarse mesh 10 × 20.
To summarize, the micropolar model demonstrates a good convergence property and is capable of solving the mesh dependency problems significantly in simulating strain softening behavior and strain localization phenomena.

Conclusion
In present study, a traditional elastoplastic model for granular material based on critical state has been improved by micropolar approach to take into account the micro-rotations. The finite element implementation of the polarized model has been formulated in detail.
After the illustrations of model, parameters calibration are carried out by fitting the experimental data of F-75 Ottawa sand. With the calibrated parameters, it has been verified that the micropolar model is capable of predicting well the results in laboratory biaxial tests. In contrast, if the classical model is used to model strain softening or strain localization phenomena, the solutions in postbifurcation regime are pathological and mesh dependent. At last, by discussing the simulated shear bands in biaxial tests with the classical model as well as with the polarized model, the great performance of micropolar approach in regularizing numerical solutions or relieving mesh dependency problems has been verified. Not only the shear band patterns in terms of the thickness and orientation but also the mechanical responses are independent on the discretization in finite element analysis. Most notably, the physically meaningful grain rotations considered in a micropolar model is also one of the most attracting factors in present study. With the outstanding features of physical meaning, convergence efficiency and finite element mesh independency, there is no doubt the polarized model can be used to model more strain localized failure modes in geotechnical structures, such as slope, shallow foundation, passive retaining walls, among others.

Conflicts of interests
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the Research Article SN Applied Sciences (2021) 3:725 | https://doi.org/10.1007/s42452-021-04708-z source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.