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

Jiangxin Liu, Zhen-Yu Yin, Lijian Wu, Pierre-Yves Hicher

► To cite this version:

Jiangxin Liu, Zhen-Yu Yin, Lijian Wu, Pierre-Yves Hicher. Finite element implementation and application of a sand model in micropolar theory. SN Applied Sciences, 2021, 3 (8), pp.725. 10.1007/s42452-021-04708-z. hal-04737844

HAL Id: hal-04737844

https://hal.science/hal-04737844v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

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

Jiangxin Liu¹ID · Zhen-Yu Yin² · Lijian Wu¹ · Pierre-Yves Hicher³

Received: 2 April 2021 / Accepted: 22 June 2021

Published online: 02 July 2021

© The Author(s) 2021 OPEN

## Abstract

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.

## 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 rotations in present paper.

(2) 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 phenomenon in 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.

**Keywords** Micropolar technique · Numerical simulations · Strain localization · Mesh independency · Granular material

<figure><img src="image_3.png" /></figure>

<footnote>Jiangxin Liu, liu.jx@rioh.cn | ¹Research Institute of Highway Ministry of Transport, Beijing 100088, China. ²Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong. ³Institut de Recherche en Génie Civil Et Mécanique (GeM), Ecole Centrale de Nantes, Université de Nantes, UMR CNRS 6183, Nantes, France.</footnote>
# 1 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 high-gradient 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 conditions with the micropolar J₂ 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₂ 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.

# 2 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.

$$ \mathbf{u} = [u_x \ u_y \ \omega_z]^T \qquad (1) $$

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,

$$ \epsilon = [\epsilon_{xx} \ \epsilon_{yy} \ \epsilon_{zz} \ \epsilon_{xy} \ \epsilon_{yx} \ \kappa_{zx} l_c \ \kappa_{zy} l_c]^T \qquad (2) $$

$$ \sigma = [\sigma_{xx} \ \sigma_{yy} \ \sigma_{zz} \ \sigma_{xy} \ \sigma_{yx} \ m_{zx}/l_c m_{zy}/l_c]^T \qquad (3) $$
in which $m_{zx}$ and $m_{zy}$ are couple stresses corresponding to the micro curvatures $\kappa_{zx}$ and $\kappa_{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].

## 2.1 Equilibrium equations

The unified equilibrium governing equations for both static and dynamic problems in micropolar theory can be expressed

$$ \begin{cases} \sigma_{ij,i} + f_j = 0 \\ m_{kj,j} + e_{kij} \sigma_{ij} + c_k = 0 \end{cases} \qquad (4) $$

in which $\sigma_{ij}$ and $m_{kj,j}$ denote Cauchy stress and micro-couple 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 $\tau_S = (\sigma_{xy} + \sigma_{yx})/2$ produces shear strain whereas the skew part $\tau_A = (\sigma_{xy} - \sigma_{yx})/2$ only causes rotations.

## 2.2 Kinematics equations

In micropolar theory, the formulations of normal strain components $\epsilon_{xx}$ and $\epsilon_{yy}$ are identical to the forms in classical continuum theory, whereas, the shear strain components $\epsilon_{xy}$ and $\epsilon_{yx}$ are linked not only to the translations but also to the rotation. The curvatures $\kappa_{xz}$ and $\kappa_{yz}$ are calculated from the gradients of rotational quantity $\omega_z$.

$$ \epsilon_{xy} = \frac{\partial u_y}{\partial x} - \omega_z, \quad \epsilon_{yx} = \frac{\partial u_x}{\partial y} + \omega_z \qquad (5) $$

$$ \kappa_{xz} = \frac{\partial \omega_z}{\partial x}, \quad \kappa_{yz} = \frac{\partial \omega_z}{\partial y} \qquad (6) $$

## 2.3 Elastic constitutive laws

Summarize the equilibrium and kinematic equations in matrix-vector forms,

$$ \epsilon = \mathbf{L}u \qquad (7) $$

$$ \mathbf{L}^{\mathrm{T}} \boldsymbol{\sigma} + \mathbf{f} = 0 \qquad (8) $$

where $\mathbf{u}$ (including rotation) and $\mathbf{f}$ (including body moment) denote displacement vector and body force vector, respectively, and $\mathbf{L}$ is named strain operator matrix,

$$ \mathbf{L}^{\mathrm{T}} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 & 0 & \frac{\partial}{\partial y} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 & \frac{\partial}{\partial x} & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 1 & l_c \frac{\partial}{\partial x} & l_c \frac{\partial}{\partial y} \end{bmatrix} \qquad (9) $$

The strain rate is related to stress rate through elastic stiffness matrix,

$$ \dot{\boldsymbol{\epsilon}} = \mathbf{D}^{\mathrm{e}} \dot{\boldsymbol{\epsilon}}^{\mathrm{e}} \qquad (10) $$

$$ \mathbf{D}^{\mathrm{e}} = \begin{bmatrix} \lambda + 2G & \lambda & \lambda & 0 & 0 & 0 & 0 \\ \lambda & \lambda + 2G & \lambda & 0 & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2G & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & G + G_c & G - G_c & 0 & 0 \\ 0 & 0 & 0 & G - G_c & G + G_c & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2G & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2G \end{bmatrix} \qquad (11) $$

in which $\lambda = 2 Gu/(1-2\nu)$ $\nu$ is Poisson's ratio, $G$ and $G_c$ denote shear modulus and micropolar shear modulus, respectively.

## 2.4 Elastoplastic constitutive laws

For small strain micropolar continuum, it is supposed that total strain includes elastic part and plastic part.

$$ \dot{\boldsymbol{\epsilon}} = \dot{\boldsymbol{\epsilon}}^{\mathrm{e}} + \dot{\boldsymbol{\epsilon}}^{\mathrm{p}} \qquad (12) $$

Substitute Eqs. (12) into (10), we get

$$ \dot{\boldsymbol{\sigma}} = \mathbf{D}^{\mathrm{e}}(\dot{\boldsymbol{\epsilon}} - \dot{\boldsymbol{\epsilon}}^{\mathrm{p}}) \qquad (13) $$

The plastic strain is formulated

$$ \dot{\boldsymbol{\epsilon}}^{\mathrm{p}} = \dot{\boldsymbol{\lambda}}\mathbf{m}, \quad \text{with} \quad \mathbf{m} = \frac{\partial \mathbf{Q}}{\partial \boldsymbol{\sigma}} \qquad (14) $$

in which $\dot{\boldsymbol{\lambda}}$ is named plastic multiplier, $\mathbf{m}$ denotes the direction of plastic flow and $\mathbf{Q}$ is the potential function. The consistency condition controls plastic flow. In plasticity theory, the Kuhn-Tucker must be guaranteed

$$ F(\boldsymbol{\sigma}, \mathbf{\kappa}) = \phi(\boldsymbol{\sigma}) - \bar{\sigma}(\mathbf{\kappa}) = 0 \qquad (15) $$

$$ \begin{cases} \dot{\lambda} \ge 0 \\ F \le 0 \\ \dot{\lambda}F = 0 \end{cases} \qquad (16) $$

$$ \dot{\mathbf{F}} = \mathbf{n}^{\mathrm{T}} \dot{\boldsymbol{\sigma}} - h \dot{\boldsymbol{\lambda}} = 0 \qquad (17) $$

where $\mathbf{\kappa}$ contains hardening variable vector, and $\mathbf{n} = \frac{\partial \mathbf{F}}{\partial \boldsymbol{\sigma}}$ and $h = -\frac{\dot{\boldsymbol{\kappa}}}{\dot{\boldsymbol{\lambda}} \frac{\partial \mathbf{F}}{\partial \boldsymbol{\kappa}}}$. Then, the stress-strain relation in a micropolar continuum can be obtained by combining Eqs. (13) and (17)
$$ \dot{\sigma} = \left[ \mathbf{D}^e - \frac{\mathbf{D}^e \mathbf{m} \mathbf{n}^T \mathbf{D}^e}{h + \mathbf{n}^T \mathbf{D}^e \mathbf{m}} \right] \dot{\mathbf{e}} \qquad (18) $$
$$ \begin{gathered} \frac{\partial g}{\partial \sigma'_{ij}} = \frac{\partial g}{\partial p'} \frac{\partial p'}{\partial \sigma'_{ij}} + \frac{\partial g}{\partial q} \frac{\partial q}{\partial \sigma'_{ij}} \\ \text{with } \frac{\partial g}{\partial p'} = A_d \left( M_{pt} - \frac{q}{p'} \right); \frac{\partial g}{\partial q} = 1 \end{gathered} \qquad (25) $$

# 3 Extension of classical model to a micropolar one

## 3.1 Introduction of the classical model

The models for granular materials have been improved from linear elasticity to ideal plasticity and then to non-linear 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

$$ \dot{\epsilon}_v^e = \frac{\dot{p}'}{K}, \quad \dot{\epsilon}_d^e = \frac{\dot{q}}{3G} \qquad (19) $$

The forms of bulk modulus K and shear modulus G are expressed [12]

$$ K = K_0 \cdot p_{at} \frac{(2.97 - e)^2}{1 + e} \left( \frac{p'}{p_{at}} \right)^{\zeta} \qquad (20) $$

$$ G = G_0 \cdot p_{at} \frac{(2.97 - e)^2}{1 + e} \left( \frac{p'}{p_{at}} \right)^{\zeta} \qquad (21) $$

in which e denotes void ratio, u 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 v,

$$ \frac{K}{G} = \frac{2(1 + v)}{3(1 - 2v)} \qquad (22) $$

The plastic strain is formulated

$$ \dot{\epsilon}_{ij}^p = d\lambda \frac{\partial g}{\partial \sigma'_{ij}} \qquad (23) $$

The yield function is expressed

$$ f = \frac{q}{p'} - \frac{M_p \epsilon_d^p}{k_p + \epsilon_d^p} \qquad (24) $$

in which k_p', controlling the initial slope of the curve q/p' - ε_d^p, 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).

$$ \tan \phi_p = \left( \frac{e_c}{e} \right)^{n_p} \tan \phi_u, \quad \tan \phi_{pt} = \left( \frac{e_c}{e} \right)^{-n_d} \tan \phi_u \qquad (26) $$

And we calculate void ratio e_c on the critical state line

$$ e_c = e_{ref} - \lambda \ln \left( \frac{p'}{p_{at}} \right) \qquad (27) $$

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.

<figure><img src="image_4.png" /><figcaption>Fig. 1 Principle of the critical state based nonlinear sand model</figcaption></figure>
$$d\lambda = \frac{\left[ \frac{\partial f}{\partial \sigma_{ij}} \right]^T D_{ijkl} d\varepsilon_{kl}}{\left[ \frac{\partial f}{\partial \sigma_{ij}} \right]^T D_{ijkl} \frac{\partial g}{\partial \sigma_{kl}} - \frac{\partial f}{\partial \varepsilon_d^p} \frac{\partial g}{\partial q}} \quad (28)$$

## 3.2 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.

$$\varepsilon_d^p = \left[ \frac{1}{3} \dot{e}_{ij}^p \dot{e}_{ij}^p + \frac{1}{3} \dot{e}_{ij}^p \dot{e}_{ji}^p + \frac{2}{3} \dot{\kappa}_{ij}^p \dot{\kappa}_{ij}^p l_c^2 \right]^{\frac{1}{2}} \quad (29)$$

$$J_2 = \frac{1}{4} s_{ij} s_{ij} + \frac{1}{4} s_{ij} s_{ji} + \frac{1}{2} m_{ij} m_{ji} / l_c^2 \quad (30)$$

In the formula, the plastic micro-curvature rate tensor and the plastic deviatoric strain rate tensor are represented by $\dot{\kappa}_{ij}^p$ and $\dot{e}_{ij}^p$, 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

$$q = \sqrt{\frac{1}{2} \boldsymbol{\sigma}^T \mathbf{P} \boldsymbol{\sigma}} \quad (31)$$

where **P** matrix is defined

$$\mathbf{P} = \begin{bmatrix} 2 & -1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 & 0 \\ -1 & -1 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3/2 & 3/2 & 0 & 0 \\ 0 & 0 & 0 & 3/2 & 3/2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 \end{bmatrix} \quad (32)$$

Similarly, the equivalent plastic strain $\varepsilon_d^p$ is newly updated in matrix-vector form

$$\varepsilon_d^p = \sqrt{\frac{2}{3} (\boldsymbol{\varepsilon}^p)^T \mathbf{Q} \boldsymbol{\varepsilon}^p} \quad (33)$$

with the definition of matrix **Q**

$$\mathbf{Q} = \begin{bmatrix} 2/3 & -1/3 & -1/3 & 0 & 0 & 0 & 0 \\ -1/3 & 2/3 & -1/3 & 0 & 0 & 0 & 0 \\ -1/3 & -1/3 & 2/3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/2 & 1/2 & 0 & 0 \\ 0 & 0 & 0 & 1/2 & 1/2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \quad (34)$$

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.

# 4 Finite element implementation of the polarized model

## 4.1 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 $\boldsymbol{\delta}^e$ at each node via the interpolation function **N**.

$$\mathbf{u} = \mathbf{N} \boldsymbol{\delta}^e \quad (35)$$

$$\mathbf{u} = \begin{Bmatrix} u \\ v \\ \omega_z \end{Bmatrix}, \quad \boldsymbol{\delta}^e = \begin{Bmatrix} \delta_1 \\ \delta_2 \\ \delta_3 \\ \delta_4 \\ \delta_5 \\ \delta_6 \\ \delta_7 \\ \delta_8 \end{Bmatrix}, \quad \delta_i = \begin{Bmatrix} u_i \\ v_i \\ \omega_{zi} \end{Bmatrix} \quad (i = 1, 2, \dots, 8) \quad (36)$$

$$\mathbf{N} = [\mathbf{I}_1 \ \mathbf{I}_2 \ \mathbf{I}_3 \ \mathbf{I}_4 \ \mathbf{I}_5 \ \mathbf{I}_6 \ \mathbf{I}_7 \ \mathbf{I}_8], \quad \text{with } \mathbf{I} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (37)$$

The interpolation functions are defined

<figure><img src="image_5.png" /><figcaption>Fig. 2 Plane strain element of micropolar continuum: a 8-node element; b reduced integration</figcaption></figure>
$$ \left\{ \begin{array}{ll} N_1 = \frac{1}{4}(1 - \xi)(1 - \eta) - \frac{1}{2}(N_5 + N_8), & N_5 = \frac{1}{2}(1 - \xi^2)(1 - \eta) \\ N_2 = \frac{1}{4}(1 + \xi)(1 - \eta) - \frac{1}{2}(N_5 + N_6), & N_6 = \frac{1}{2}(1 + \xi)(1 - \eta^2) \\ N_3 = \frac{1}{4}(1 + \xi)(1 + \eta) - \frac{1}{2}(N_6 + N_7), & N_7 = \frac{1}{2}(1 - \xi^2)(1 + \eta) \\ N_4 = \frac{1}{4}(1 - \xi)(1 + \eta) - \frac{1}{2}(N_7 + N_8), & N_8 = \frac{1}{2}(1 - \xi)(1 - \eta^2) \end{array} \right. \qquad (38) $$

Then, the strain components can be derived from the displacement and rotation.

$$ \varepsilon = \mathbf{L}\mathbf{u} = \mathbf{L}\mathbf{N}\delta^{\mathbf{e}} = \mathbf{B}\delta^{\mathbf{e}} \qquad (39) $$

B = LN is the strain matrix. Then the stress can be obtained according to constitutive law

$$ \sigma = \mathbf{D}^{\mathrm{ep}}\varepsilon = \mathbf{D}^{\mathrm{ep}}\mathbf{B}\delta^{\mathbf{e}} \qquad (40) $$

The global coordinates of the element can also be obtained by interpolating nodes' coordinates.

$$ \left\{ \begin{array}{l} x = \sum_{i=1}^{8} N_i x_i = N_1 x_1 + N_2 x_2 + \cdots + N_8 x_8 \\ y = \sum_{i=1}^{8} N_i y_i = N_1 y_1 + N_2 y_2 + \cdots + N_8 y_8 \end{array} \right. \qquad (41) $$

After the partial differential operation, we get

$$ \left\{ \begin{array}{l} \left[ \frac{\partial N_i}{\partial \xi} \quad \frac{\partial x}{\partial \xi} \quad \frac{\partial y}{\partial \xi} \right] \left\{ \left[ \frac{\partial N_i}{\partial x} \quad \frac{\partial N_i}{\partial y} \right] \right\} = [\mathbf{J}] \left\{ \frac{\partial N_i}{\partial y} \right\} \\ \left[ \frac{\partial N_i}{\partial \eta} \quad \frac{\partial x}{\partial \eta} \quad \frac{\partial y}{\partial \eta} \right] = [\mathbf{J}] \left\{ \frac{\partial N_i}{\partial x} \right\} \end{array} \right. \qquad (42) $$

The Jacobian matrix [J], used to realize the map between (ξ, η) and (x, y), is formulated

$$ [\mathbf{J}] = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial \sum N_i x_i}{\partial \xi} & \frac{\partial \sum N_i y_i}{\partial \xi} \\ \frac{\partial \sum N_i x_i}{\partial \eta} & \frac{\partial \sum N_i y_i}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial N_1}{\partial \eta} & \frac{\partial N_2}{\partial \eta} & \frac{\partial N_8}{\partial \eta} \\ \frac{\partial N_1}{\partial \xi} & \frac{\partial N_2}{\partial \xi} & \frac{\partial N_8}{\partial \xi} \\ \frac{\partial N_1}{\partial \eta} & \frac{\partial N_2}{\partial \eta} & \frac{\partial N_8}{\partial \eta} \end{bmatrix} \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ \vdots & \vdots \\ x_8 & y_8 \end{bmatrix} \qquad (43) $$

## 4.2 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

$$ \prod_p = \int_{\Omega} \frac{1}{2} \varepsilon^{\mathrm{T}} \mathbf{D}^{\mathrm{ep}} \varepsilon \mathrm{tdxdy} - \int_{\Omega} \mathbf{u}^{\mathrm{T}} \mathrm{ftdxdy} - \int_{S} \mathbf{u}^{\mathrm{T}} \mathrm{Tds} \qquad (44) $$

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.

$$ \prod_p = \sum_e \prod_p^e = \sum_e \left( (\delta^e)^T \int_{\Omega_e} \frac{1}{2} \mathbf{B}^T \mathbf{D}^{\mathrm{ep}} \mathbf{B} \mathrm{tdxdy} \delta^e \right) - \sum_e \left( \int_{\Omega_e} (\delta^e)^T \mathbf{u}^T \mathbf{f} \mathrm{tdxdy} \right) - \sum_e \left( \int_{S_e} (\delta^e)^T \mathbf{N}^T \mathbf{T} \mathrm{ds} \right) \qquad (45) $$

For a random virtual displacement, it must be required the partial differential $\frac{\partial \prod_p}{\partial \delta^e} = 0$. In this way, the equilibrium equations of an element can be formulated

$$ \mathbf{K}^{\mathbf{e}} \delta^{\mathbf{e}} = \mathbf{P}^{\mathbf{e}} \qquad (46) $$

$$ \mathbf{K}^{\mathbf{e}} = \int_{\Omega_e} \mathbf{B}^{\mathrm{T}} \mathbf{D}^{\mathrm{ep}} \mathbf{B} \mathrm{tdxdy} \qquad (47) $$

$$ \mathbf{P}^{\mathbf{e}} = \int_{\Omega_e} \mathbf{u}^{\mathrm{T}} \mathbf{f} \mathrm{tdxdy} + \int_{S_e} \mathbf{N}^{\mathrm{T}} \mathbf{T} \mathrm{ds} \qquad (48) $$

$\mathbf{K}^{\mathbf{e}}$ and $\mathbf{P}^{\mathbf{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

$$ \mathbf{K} \mathbf{a} = \mathbf{P} \qquad (49) $$

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 New-

ton-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
$$ \mathbf{K}^e = \int_{\Omega_e} \mathbf{B}^T \mathbf{D}^{ep} \mathbf{B} t dxdy = \int_{-1}^{1} \int_{-1}^{1} \mathbf{B}^T \mathbf{D}^{ep} \mathbf{B} |\mathbf{J}| d\xi d\eta = \sum_{i=1}^{m} \sum_{j=1}^{n} H_i H_j \mathbf{B}^T \mathbf{D}^{ep} \mathbf{B} t |\mathbf{J}| \quad (50) $$

in which $H_i$ or $H_j$ represents the weight factors for interpolation.

All the formulations of the polarized sand model are programed in the User defined Element (UEL) procedure. The computation flow chart can be found in Fig. 3.

## 5 Calibration of model parameters

The polarized model includes 12 parameters, which are classified into four groups: (a) parameters about elasticity $\zeta$, $\nu$, and $K_0$, (b) parameters defining critical state $\lambda$, $\phi_u$ and $e_{ref}$, (c) parameters interlocking plasticity $n_p$, $n_d$, $A_d$ and $k_p'$ (d) and the particular parameters $I_c$ and $G_c$ in the micropolar model.

The elastic parameters $\zeta$, $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

triaxial drained test. To determine the critical state parameters $\lambda$, $\phi_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 $\nu$ can be assigned based on materials.

For the determination of the micropolar parameter $I_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 $I_c$ in numerical simulations [3, 13–16]. Therefore, $I_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 curve-fitting 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. 4 ($e_{max}$ = 0.805, $e_{min}$ = 0.486, $G_s$ = 2.65, $d_{50}$ = 0.22 mm) was

<figure><img src="image_6.png" /><figcaption>Fig. 3 Computation flow chart of UEL</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4 Particle size distribution of F-75 sand</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5 Calibration of elastic parameters with an isotropic compression test on F-75 sand</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6 Calibration of plastic parameters with at least three triaxial
drained tests on F-75 sand: **a** vertical axial strain vs. deviatoric
stress; **b** mean stress versus deviatoric stress; **c** vertical axial strain
versus void ratio; **d** mean stress versus void ratio</figcaption></figure>

selected to identify the model parameters. The experi-
mental data shown in Figs. 5 and 6 of one isotropic com-
pression 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 speci-
men's deformation is monitored by noting the displace-
ments 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
(u, Ic and Gc) can be set in advance. It requires at least three
drained triaxial tests with different initial confining pres-
sure 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 behav-
ior of F-75 dry sand in isotropic compression and triaxial
tests can be very well predicted by the micropolar model.

**6 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 simu-
lations by the classical model were proven to be pathologi-
cal and mesh dependent in post-bifurcation regime. There-
fore, the enhanced micropolar model was adopted herein
to describe the granular soil behavior. I_c was regarded as
$d_{50}$ of F-75 sand (0.22 mm) and $G_c$ was half of G. Other calibrated model parameters are listed in Table 1. Notably, the critical friction angle $\phi_u$ has been slightly increased as shown in Table 2. The experimental findings of Alshibli et al. [21] can support the improvement of $\phi_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 $\phi_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.

## 7 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

<table><caption>**Table 1** Calibrated constitutive parameters</caption><thead><tr><th>Parameters</th><th>K<sub>0</sub></th><th>ζ</th><th>φ<sub>u</sub></th><th>e<sub>ref</sub></th><th>λ</th><th>k<sub>p</sub></th><th>A<sub>d</sub></th><th>n<sub>p</sub></th><th>n<sub>d</sub></th></tr></thead><tbody><tr><td>Values</td><td>60</td><td>0.63</td><td>35.8</td><td>0.776</td><td>0.015</td><td>0.004</td><td>0.4</td><td>1</td><td>2</td></tr></tbody></table>

<table><caption>**Table 2** Model parameters used to simulate biaxial tests</caption><thead><tr><th>Parameters</th><th>K<sub>0</sub></th><th>ζ</th><th>φ<sub>u</sub></th><th>e<sub>ref</sub></th><th>λ</th><th>k<sub>p</sub></th><th>A<sub>d</sub></th><th>n<sub>p</sub></th><th>n<sub>d</sub></th></tr></thead><tbody><tr><td>Values</td><td>60</td><td>0.63</td><td>38.5</td><td>0.776</td><td>0.015</td><td>0.004</td><td>0.4</td><td>1</td><td>2</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Comparisons between laboratory test and simulation for very dense sand under 100 kPa confining pressure at 10% axial strain: a laboratory test; b simulation by the micropolar model</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Comparisons between laboratory tests and numerical simulations for medium dense F-75 sand: **a** axial strain versus principle stress ratio; **b** axial strain versus volumetric strain</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 9** Comparisons between laboratory tests and numerical simulations for very dense F-75 sand: **a** axial strain versus principle stress ratio; **b** axial strain versus volumetric strain</figcaption></figure>

then imposed shear load on the model top by continuous incremental displacement in the second stage. To investigate the mesh dependency problems, four mesh densities from coarse to fine, i.e. mesh 10×20, mesh 15×30, mesh 20×40 and mesh 30×60, have been discussed. Mesh 10×20 means the model is uniformly discretized to 20 and 10 elements in length and width direction, respectively.

## 7.1 Mesh dependent problems by the classical model

### 7.1.1 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
<figure><img src="image_15.png" /><figcaption>Fig. 10 Shear bands in different meshes by the classical model: a mesh 10×20; b mesh 15×30; c mesh 20×40; d mesh 30×60</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 11 Displacement versus load of different meshes by the classical model</figcaption></figure>

and pathological solutions. 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.

### 7.1.2 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

<figure><img src="image_17.png" /><figcaption>Fig. 12 Shear band inclined angles in different meshes by the classical model: a mesh 10×20 β₁ = 52.69°; b mesh 15×30 β₂ = 57.65°; c mesh 20×40 β₃ = 60.15°</figcaption></figure>
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.

## 7.2 Mesh independent performance by the micropolar model

### 7.2.1 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

<figure><img src="image_18.png" /><figcaption>Fig. 14 Displacement versus load of different meshes by the micropolar model</figcaption></figure>

relieved even in the softening regime. Except the too coarse mesh 10×20, the mechanical curves of other three meshes are completely consistent with each other.

### 7.2.2 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 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.

<figure><img src="image_19.png" /><figcaption>Fig. 13 Shear bands in different meshes by the micropolar model: a mesh 10×20; b mesh 15×30; c mesh 20×40; d mesh 30×60</figcaption></figure>
<caption>Fig. 15 Shear band inclined angles in different meshes:</caption>

a mesh 10 × 20 $β_1$ = 53.10°;
b mesh 15 × 30 $β_2$ = 53.22°;
c mesh 20 × 40 $β_3$ = 53.22°; d mesh 30 × 60 $β_4$ = 53.22°

<figure><img src="image_20.png" /></figure>

**8 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.

**Funding** Dr. Jiangxin Liu was funded by Research Institute of High-way Ministry of Transport (CN).

**Declarations**

**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
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://creativecommons.org/licenses/by/4.0/.

## References

1. Oda M, Kazama H (1998) Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils. Géotechnique 48(4):465-481. https://doi.org/10.1680/geot.1998.48.4.465

2. Viggiani C, Bésuelle P, Hall S, Desrues J (2010) Imaging sand deformation at the grain scale. In: EPJ Web of Conferences: EDP Sciences. p. 22021

3. Alsaleh MI (2004) Numerical modeling of strain localization in granular materials using Cosserat theory enhanced with micro-fabric properties. Dissertation, Jordan University of Science and Technology

4. Alshibli KA, Alsaleh MI, Voyiadjis GZ (2006) Modelling strain localization in granular materials using micropolar theory: numerical implementation and verification. Int J Numer Anal Meth Geomech 30(15):1525-1544. https://doi.org/10.1002/nag.534

5. Arslan H, Sture S (2008) Finite element analysis of localization and micro-macro structure relation in granular materials. Part II: Implementation and simulations. Acta Mech 197(3):153-171. https://doi.org/10.1007/s00707-007-0514-0

6. de Borst R, Sluys L (1991) Localisation in a Cosserat continuum under static and dynamic loading conditions. Comput Methods Appl Mech Eng 90(1-3):805-827. https://doi.org/10.1016/0045-7825(91)90185-9

7. Huang W, Bauer E (2003) Numerical investigations of shear localization in a micro-polar hypoplastic material. Int J Numer Anal Meth Geomech 27(4):325-352. https://doi.org/10.1002/nag.275

8. Li X, Tang H (2005) A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modeling of strain localisation. Comput Struct 83(1):1-10. https://doi.org/10.1016/j.compstruc.2004.08.009

9. Eringen AC (1976) Nonlocal polar field theories. Contin Phys 4(Part III):205-264

10. Liu J (2018) Numerical investigations of the strain localization in geotechnical engineering within the framework of micropolar theory. Dissertation, Ecole centrale de Nantes

11. Yin Z-Y, Hicher P-Y, Jin Y-F (eds) (2020) Elastoplastic Modeling of Soils: From Mohr-Coulomb to SIMSAND. In: Practice of constitutive modelling for saturated soils, pp 121-176

12. Richard FE, Hall JR, Woods RD (1970) Vibrations of soils and foundations

13. Alshibli KA, Batiste SN, Sture S, Curreri PA (2002) Experimental Observations of Localization Phenomena in Sands: Plane Strain Versus Triaxial Compression Conditions

14. Desrues J (1990) Shear band initiation in granular materials: experimentation and theory. Geomate Const Equ Model 283-310

15. Desrues J, Viggiani G (2004) Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry. Int J Numer Anal Meth Geomech 28(4):279-321. https://doi.org/10.1002/nag.338

16. Roscoe KH (1970) The influence of strains in soil mechanics. Géotechnique 20(2):129-170. https://doi.org/10.1680/geot.1970.20.2.129

17. de Borst R, Mühlhaus H (1991) Continuum models for discontinuous media. In: Fracture Processes in Concrete, Rock and Ceramics Proceedings of the International rilem/esis Conference, Noordwijk, the Netherlands, June 19-21 2, 601-618(1991)

18. Tang H, Li X (2008) Numerical analysis for the effects of constitutive parameters in Cosserat continuum model on the simulation results of the strain localization. Chin J Comput Mech 25:676-681

19. Jin Y-F, Yin Z-Y, Shen S-L, Hicher P-Y (2016) Investigation into MOGA for identifying parameters of a critical-state-based sand model and parameters correlation by factor analysis. Acta Geotech 11(5):1131-1145

20. Omar, T (2010) Specimen size effect on shear behavior of loose sand in triaxial testing. Dissertation, The University of Western Ontario

21. Alshibli KA, Batiste SN, Sture S (2003) Strain localization in sand: plane strain versus triaxial compression. J Geotech Geoenviron Eng 129(6):483-494. https://doi.org/10.1061/(ASCE)1090-0241(2003)129:6(483)

22. Alshibli KA, Sture S (2000) Shear band formation in plane strain experiments of sand. J Geotech Geoenviron Eng 126(6):495-503. https://doi.org/10.1061/(ASCE)1090-0241(2000)126:6(495)

23. Bažant ZP, Belytschko TB, Chang T-P (1984) Continuum theory for strain-softening. J Eng Mech 110(12):1666-1692. https://doi.org/10.1061/(ASCE)0733-9399(1984)110:12(1666)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.