On the observation of extended microstructures and localized deformations in granular materials - A feature of excat quasi-convex energy envelope

Abstract

In this paper a most promising scope of the exact quasiconvex energy envelope in modeling the granular materials with microstructures is presented. This study shows that it is possible to observe both the extended microstructures and localized deformations in granular materials using a variational model based on the mathematical relaxation theory. The variational model is derived within the framework of Cosserat continuum. The computational algorithm based on finite element method is used to carry out numerical computations. The features of the proposed model are studied for three representative examples: the Couette shear cell, the rectangular specimen in compression and the indentation of a granular structure. The obtained results demonstrate on the possible applications and features of exact quasiconvex energy envelops in modeling the granular materials.

Keywords:
Granular materials; energy methods; extended microstructures; localized deformations; shear bands; relaxed Cosserat energy; and Finite element method

1 INTRODUCTION

Energy relaxation methods in analysis of materials with microstructures have been extensively studied for a large number of engineering problems in mechanics in the context of elasticity (Ball, 1976; Ball and James, 1987; DeSimone and Dolzmann, 2002; Govindjee et al., 2007; Khan and Hackl, 2018; Kohn, 1991; Kohn and Strang, 1986a,b; Raoult, 2010), elasto-plasticity (Conti and Theil, 2005; Dret and Raoult, 1995; Gürses and Miehe, 2011; Hackl and Heinen, 2008) and plasticity (Carstensen etal., 2008, 2002; Conti et al., 2009, 2007; Conti and Ortiz, 2005). In modeling the mechanical behavior of such materials it isalways important to construct a relaxed energy envelope of the corresponding non-quasiconvex energy potential. So far, there are very few studies (Conti et al., 2007; Conti and Ortiz, 2005; Conti and Theil, 2005; DeSimone and Dolzmann, 2002; Dret and Raoult, 1995; Khan and Hackl, 2018; Kohn, 1991; Kohn and Strang, 1986a), (Kohn and Strang, 1986b; Kohn and Vogelius, 1987; Raoult, 2010) in literature where the computation of exact quasiconvex envelope was possible.

Quasiconvex energy envelopes in energy minimization problems are of interest in every physical situation where there arefine scale oscillations of gradients of the infimizing deformations. Fluctuations of these infimizing deformations at fine scale leads to the development of material microstructures. Development and formation of these microstructures in granular materials is a debatable subject and have attracted a great interest within the framework of Cosserat continuum Alsaleh et al. (2006); Alshibli et al. (2006); Bardet(1994); deBorst (1991); Lambrecht et al. (2003); Miehe et al. (2004); Oda and Kazama (1998), (Sulem and Derrolaza, 2002; Tejchman and Bauer, 1996; Tejchman and Gudehus (1993); Tordesillas et al., 2004), (Tordesillas et al., 2005; Trinch and Hackl (2014)) in the past.

Cosserat elasticity Forest and Sievert (2006); Kafadar and Eringen (1976); Khan and Hackl (2012); Schaefer (1967); Steinmann (1994) in comparison to Boltzmann continuum theory incorporates independent rotational degrees of freedom at eachmaterial point. Forthis reason it can adequately describes the microstructural behavior of granular materials which are discrete in nature whereeach grain has micro-rotations independent of its macro-rotations. It is therefore widely acceptable and successfully usedtheory for modeling the microstructural behavior of granular materials in literature Altenbach et al. (2012); Sawada et al. (2006); Tejchman and Bauer (1996); Tejchman and Gudehus (1993); Trinch and Hackl (2014). In this work, the description of granular material deformation is based on the Cosserat continuum theory. For a review on the applications of Cosserat continuum theory the reader is referred to Altenbach et al. (2012) and the references therein.

Research has shown an agreement on the fact that formation of these microstructures can be extended microstructures Gajo et al. (2004); Gudehus and Nübel (2004) or it can appear in the form of localized deformations deBorst (1991); Trinch and Hackl (2014). In both the cases the development of these microstructural zones are highly influenced by grain rotations at the microscale Sawada et al. (2006); Tordesillas et al. (2005). In the zones of localized deformations and extended microstructures the particles organize themselves in non-permanent meso-scale structures Oda and Kazama (1998). To understand different features of particle rotations in the deformation mechanisms of granular materials the reader is referred to Bardet (1994); Oda and Kazama (1998) and the references therein. Although efforts have been made to observe and predicton the extended microstructures and localized deformations in granular materials both by experimental Aranson and Tsimring (2009); Debregeas et al. (2001); Kaus and Podladchikov (2006); Latzel et al. (2003); Savage and Sayed (1984); Sawada et al. (2006), (Howell et al., 1999; Tejchman and Gudehus, 1993); Utter and Behringer, 2009) and numerical simulations Sawada et al. (2006); Tejchman and Bauer (1996); Tordesillas et al. (2004), there is no unified variational approach available in literature that allows to compute and observe both the extended microstructure and localized deformations with the same model. However, in past the phenomena of extended microstructure formation and localized deformations have been studied as bifurcation phenomena. This study presents a relaxed computational algorithm using the finite element method that enables to predict on the formation and development of both the granular material microstructural phenomena with the same variational model in the framework of Cosserat elasticity. Moreover, this approach allows to elaborate on the different features of exact quasiconvex energy envelop in a number of physical situations.

The rest of the paper is organized as follows. In Section 2, a relaxed variational model for granular materials is presented where non-microstructural and microstructural material regimes are explicitly characterized and a formula for the relaxed energy is given. In section 3, a two field variational formulation of the proposed method is presented. In section 4, finite element descretization of the proposed method is presented and the solution algorithm to compute the deformation microstructures of granular materials is given. In Section 5, numerical results and discussion is presented demonstrating on the possible applications of the proposed methodology and the characterization of different material microstructural phases. Finally, conclusions are drawn in Section 6.

2 A RELAXED COSSERAT CONTINUUM MODEL FOR GRANULAR MATERIALS

Let$u:\Omega \subset {ℝ}^d\mapsto {ℝ}^d$be the displacement and$\phi =axl(\Phi ):\Omega \subset {ℝ}^d\mapsto {ℝ}^d$where $\Phi :\Omega \subset {ℝ}^d\mapsto so(d):=\left\{R\in {{\rm M}}^{d\times d}|R^T=-R\right\}$be the micro rotation vector field in dimension$d,$of the granular particles of the material bounded by$\Omega$with Lipschitz boundary$\partial \Omega .$Here$M^{d\times d}$is the set of$d\times d$ second order tensors. Further assume that an external force potential is defined with the considerations of external body force b, couple m, traction force t _u and traction moments t _ϕ such that

$\mathcal{l}(u,\phi )={\displaystyle {\int }_{\Omega }(b.u+m.\phi )}dV+{\displaystyle {\int }_{\partial {\Omega }_u}{t}_u.udS}+{\displaystyle {\int }_{\partial {\Omega }_{\phi }}{t}_{\phi }.\phi dS,}$ (1)

then the deformed configuration of materials can be completely determined from the following minimization problem

$\underset{u,\Phi ,\phi }{\inf }\left\{{\displaystyle {\int }_{\Omega }W(\nabla u,\Phi ,\nabla \phi )dV-\mathcal{l}(u,\phi ),{\text{u}\text{=}\text{u}}_{\circ }\text{at}\partial {\Omega }_u\text{and}\phi \text{=}{\phi }_{\circ }\text{at}\partial {\Omega }_{\phi }}\right\}\text{,}$ (2)

where $(u,\Phi ,\phi )\in {\text{W}}^{\mathrm{1,}p}(\Omega ,{ℝ}^d)\times {\text{W}}^{\mathrm{1,}p}(\Omega ,so(d)\times {\text{W}}^{\mathrm{1,}p}(\Omega ,{ℝ}^d).$ Here ${\text{W}}^{\mathrm{1,}p}$is the Sobolev space of admissible deformations with$p\in (\mathrm{1,}\infty )$related to growth of the energy function$\text{W}\text{.}$ Within the framework of generalized elasticity the ansatz for the energy potential$\text{W}$in (2) is taken as

$W(\nabla u,\Phi ,\nabla \phi )=\frac{1}{2}e:ℂ:e+\frac{1}{2}\kappa :\overline{ℂ}:\kappa +i.e.p,$ (3)

Where $i.e.p$ is energy potential$\alpha {\left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)}^2$derived in Khan and Hackl (2018) and due to counter rotations of granular particles at the continuum scale. Here, $\alpha$and $\beta$are non-negative constants. Where, the macroscopic deformation tensor and the rotational strain tensor are defined as$e=\nabla u-\Phi$and$\kappa =\nabla \phi$respectively. The fourth order constitutive tensors C and C of elastic constants are determined as

$ℂ(\nabla u,\nabla \phi )=\frac{{\partial }^{2}W(\nabla u,\Phi ,\nabla \phi )}{\partial \nabla u\otimes \partial \nabla u}\text{and}\overline{ℂ}(\nabla u,\nabla \phi )=\frac{{\partial }^{2}W(\nabla u,\Phi ,\nabla \phi )}{\partial \nabla \phi \otimes \partial \nabla \phi }.$ (4)

Within the framework of an isotropic elastic Cosserat medium the enhanced energy function in (3) takes the form

$W(\nabla u,\Phi ,\nabla \phi )=\{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+\mu {\Vert \text{dev}\epsilon \Vert }^2+{\mu }_c{\Vert asy\nabla \text{u-}\Phi \Vert }^2+\frac{\overline{\lambda }}{2}{\left(tr\kappa \right)}^2\\ +\overline{\mu }{\Vert sym\kappa \Vert }^2+{\overline{\mu }}_c{\Vert asy\kappa \Vert }^2+\alpha {\left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)}^2\end{array}$ (5)

where$\lambda$is the classical dilatancy parameter, $\mu$is the classical elastic shear modulus, ${\mu }_c$is the Cosserat shear modulus,$\overline{\lambda }$is the Cosserat material dilatancy parameter, $\overline{\mu }$ is the bending modulus,${\overline{\mu }}_c$is the coupled shear modulus, $\epsilon$is the Cauchy strain tensor, and$\kappa$is the curvature strain tensor.

The energy potential in equation (5) is non-quasiconvex, hence leading to non-attainment of the minimizers in the energy minimization problem (2). These non-attainments of the minimizers are essentially due the possible fluctuations in the displacement and microrotation fields at the fine scales. Physically, these field fluctuations can be seen Gürses and Miehe (2011); Miehe et al. (2004) as distortions of the finite element meshes. The fine scale oscillations of the minimizing displacement and microrotation field variables will lead to the development of micro-structures in the material. Formation of such internal structures can be extended microstructures Gajo et al. (2004); Gudehus and Nübel (2004) which is distributed through the material domain or the localized deformation microstructures deBorst (1991); Trinch and Hackl (2014) that appear in the form of narrow shearing bands. Hence to guarantee the existence of the unique minimizing translational and microrotational deformations in this case an exact relaxed energy envelope of energy potential (5) is computed in our previous papers Khan and Hackl (2013, 2018). The result is summarized below

Relaxed Energy

The relaxed energy corresponding to energy function in (5) is divided into three different material phases each of which is characterized as follows

$\left\{u,\phi \right\}\text{is}\text{in}\{\begin{array}{l}\text{phase}\text{1}\text{if}{\Vert \kappa \Vert }^2\ge {\beta }^2{\Vert \text{dev}\epsilon \Vert }^2+\frac{\mu }{2\alpha {\beta }^2}\\ \text{phase}\text{2}\text{if}\frac{-{\mu }_{\circ }}{2\alpha }\le {\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\le \frac{\mu }{2\alpha {\beta }^2}\\ \text{phase}\text{3}\text{if}{\Vert \kappa \Vert }^2\le {\beta }^2{\Vert \text{dev}\epsilon \Vert }^2-\frac{{\mu }_{\circ }}{2\alpha }\end{array}$

Phase 1 is corresponding to the material regime where there is a microstructure in micro-rotations. Phase 2 is corresponding to the material regime where there is no internal structure of the material. Phase 3 is corresponding to the material regime where there is a microstructure in translational motions of the particles. The relaxed energy thus reads

${\text{W}}^{rel}\text{=}\{\begin{array}{l}{\text{W}}_1{}^{rel}\text{if}{\Vert \kappa \Vert }^2\ge {\beta }^2{\Vert \text{dev}\epsilon \Vert }^2+\frac{\mu }{2\alpha {\beta }^2}\\ {\text{W}}_2{}^{rel}\text{if}\frac{-{\mu }_{\circ }}{2\alpha }\le {\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\le \frac{\mu }{2\alpha {\beta }^2}\\ {\text{W}}_3{}^{rel}\text{if}{\Vert \kappa \Vert }^2\le {\beta }^2{\Vert \text{dev}\epsilon \Vert }^2-\frac{{\mu }_{\circ }}{2\alpha }\end{array}$ (6)

where${\text{W}}_1{}^{rel},{\text{W}}_2{}^{rel}$and${\text{W}}_3{}^{rel}$are explicitly given as

${\text{W}}_1{}^{rel}=\{\begin{array}{l}\{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+{\mu }_c{\Vert asy\nabla \text{u}-E\text{.}\phi \Vert }^2-\frac{{\mu }^2}{4\alpha {\beta }^4}\\ +\frac{\overline{\lambda }}{2}{(}^{\text{tr}}+(\overline{\mu }-\overline{{\mu }_c}){\Vert sym\kappa \Vert }^2+({\mu }_{\circ }+\frac{\mu }{{\beta }^2}){\Vert \kappa \Vert }^2\end{array}\text{if}\overline{\mu }\ge \overline{{\mu }_c},\\ \{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+{\mu }_c{\Vert asy\nabla \text{u}-E\text{.}\phi \Vert }^2-\frac{{\mu }^2}{4\alpha {\beta }^4}\\ +\frac{\overline{\lambda }}{2}{(}^t+(\overline{\mu }-\overline{{\mu }_c}){\Vert \text{asy}\kappa \Vert }^2+({\mu }_{\circ }+\frac{\mu }{{\beta }^2}){\Vert \kappa \Vert }^2\end{array}\text{if}\overline{\mu }<\overline{{\mu }_c},\end{array}$ (7)

${}_{\text{W}}{}^{rel}=\{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+\mu {\Vert \text{dev}\epsilon \Vert }^2+{\mu }_c{\Vert asy\nabla \text{u}-E\text{.}\phi \Vert }^2+\frac{\overline{\lambda }}{2}{\left(tr\kappa \right)}^2\\ +\overline{\mu }{\Vert sym\kappa \Vert }^2+\overline{{\mu }_c}{\Vert asy\kappa \Vert }^2+\alpha {\left({\Vert sym\kappa \Vert }^2+{\Vert asy\kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)}^2\end{array}$ (8)

and

${\text{W}}_3{}^{rel}=\{\begin{array}{l}\{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+{\mu }_c{\Vert asy\nabla \text{u}-{\rm E}\text{.}\phi \Vert }^2+\frac{\overline{\lambda }}{2}{\left(tr\kappa \right)}^2\\ +(\overline{\mu }-\overline{{\mu }_c}){\Vert sym\kappa \Vert }^2+({\mu }_{\circ }{\beta }^2+\mu ){\Vert \text{dev}\epsilon \Vert }^2-\frac{{\mu }_{\circ }{}^2}{4\alpha }\end{array}\text{if}\overline{\mu }\ge \overline{{\mu }_c},\\ \{\begin{array}{l}\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)\text{(tr}\epsilon {\text{)}}^2+{\mu }_c{\Vert asy\nabla \text{u}-{\rm E}\text{.}\phi \Vert }^2+\frac{\overline{\lambda }}{2}{\left(tr\kappa \right)}^2\\ \text{-}(\overline{\mu }-\overline{{\mu }_c}){\Vert asy\kappa \Vert }^2+({\mu }_{\circ }{\beta }^2+\mu ){\Vert \text{dev}\epsilon \Vert }^2-\frac{{\mu }_{\circ }{}^2}{4\alpha }\end{array}\text{if}\overline{\mu }<\overline{{\mu }_c},\end{array}$ (9)

respectively.

The computation of this analytical expression for the relaxed energy corresponding to non-quasiconvex energy function in (5) enable to predict all microstructural features of the material which are carried safely from the microscopic to macroscopic computational scale. Hence it is possible to extract information regarding the development of microstructural regimes in the granular materials pertinent to observing its macro-mechanical behavior. For practical applications it is now more efficient and effective to reformulate the original non-quasiconvex problem in (2) to a relaxed energy minimization problem using this relaxed potential.

3 A TWO-FIELD VARIATIONAL FORMULATION

In the quasi-static case the balance equations for the motion of a granular particle in the framework of Cosserat continumm can be derived from the variational principle (2). Taking variations of equation (2) with respect to$u$and$\phi$yields the equation for balance of linear and angular momentum respectively. The local form of these equations are thus stated as

$\nabla .\sigma +b=\mathrm{0,}$ (10)

$\nabla .\mu +\sigma :E+m=\mathrm{0,}$ (11)

Where$E$is the third order permutation tensor,$b$is the body force and$m$is the body moment, subjected to certain Drichlet boundary conditions of the type$\sigma \cdot N_u=t_u$and$\mu \cdot N_{\phi }=t_{\phi }$. The Cosserat force-stress tensor$\sigma$and couple stress tensor$\mu$are computed as the derivatives of the strain energy function with respect to the Cosserat strain tensor$e$and curvature strain tensor$\kappa ,$respectively

$\sigma =\frac{\partial W^{rel}}{\partial e},\mu \text{=}\frac{\partial W^{rel}}{\partial \kappa }.$ (12)

To derive weak balance equations or variational formulation of (10) the arbitrary variations in displacement and micro-rotation fields are chosen as$\delta u$ $\delta \phi ,$ respectively. Further, let$t_u=\sigma \cdot N_u$be the applied traction on the surface${\Gamma }_u$of the continuum body$B$having a volume$\Omega ,$with$N_u$being the normal unit vector to${\Gamma }_u$. Let$t_u$and$t_{\phi }$be the traction and moment forces on the surface${\Gamma }_u$and${\Gamma }_{\phi },$respectively, of the body, where$n_u$and$n_{\phi }$are the normal unit vectors, then by using the divergence theorem and the properties$\nabla \cdot \left(\mu \cdot \delta \phi \right)=\nabla \cdot \left(\mu \right)\cdot \delta \phi +\mu :(\nabla \otimes \delta \phi )$and$E:\sigma =\sigma :E$the final set of variational equations for (10) and (11), that must be satisfied for the equilibrium configuration of a body by its each material point in a Cosserat continuum in an integrated sense, can be written, respectively as

${\displaystyle {\int }_{\Omega }\sigma :(\nabla \otimes \delta u)}d\Omega ={\displaystyle {\int }_{\Omega }b\cdot \delta ud\Omega }+{\displaystyle {\int }_{{\Gamma }_u}{t}_u\cdot \delta \text{u}d{\Gamma }_u,}$ (13)

and

${\displaystyle {\int }_{\Omega }(\mu :(\nabla \otimes \delta \varphi )}-\text{(}\sigma :E\text{)}\cdot \delta \varphi )d\Omega ={\displaystyle {\int }_{\Omega }\text{m}\cdot \delta \varphi d\Omega }+{\displaystyle {\int }_{{\Gamma }_{\phi }}{t}_{\phi }\cdot \delta \varphi d{\Gamma }_{\phi }.}$ (14)

4 NUMERICAL IMPLEMENTATION

In the absence of body force$b$and body couple$m$the system of linear and angular momentum weak-balance equations (13) and (14) are solved numerically using finite element descretization of the material domain$\Omega$as

$\Omega ={\displaystyle \underset{e=1}{\overset{n_e}{\cup }}{\Omega }^e,{\Omega }^i\cap }{\Omega }^j=\varnothing \text{for}i\ne j$ (15)

where$n_e$is the total number of elements in$\Omega .$ It is convenient to define a local coordinate system$\left({\xi }_1,{\xi }_2,{\xi }_3\right)$in each element${\Omega }^e$such that$-1\le {\xi }_1\le \mathrm{1,}-1\le {\xi }_2\le 1$and$-1\le {\xi }_3\le 1.$ The field variables in${\Omega }^e$can be interpolated using the nodal shape functions defined in the local coordinate system. In the context of Cosserat continuum these field variables are the displacement and micro-rotation which together contains a total number of nine degrees of freedom at a single material point in an element${\Omega }^e.$ The three displacement$\left\{u_i^e;\text{}i=\mathrm{1,2,3}\right\}$and three microrotation$\left\{{\phi }_i^e;\text{}i=\mathrm{1,2,3}\right\}$degrees of freedoms are approximated on nodal points in each element${\Omega }^e$by matrices of their corresponding nodal shape functions as

$u_i{}^e\approx {\displaystyle \sum _{m=1}^N{N}^m({\xi }_1},{\xi }_2,{\xi }_3){\overline{u}}_i{}^{e}and{\phi }_i{}^e\approx {\displaystyle \sum _{m=1}^N{N}^m({\xi }_1},{\xi }_2,{\xi }_3){\overline{\phi }}_i{}^e,;i=\mathrm{1,2,3.}$ (16)

here${\left\{N^m({\xi }_1,{\xi }_2,{\xi }_3)\right\}}_{m=1}^N$are the nodal shape functions. Associated to these displacements and micro-rotations are the strain measures called Cosserat strains and curvature strains, respectively. These strains are computed in each${\Omega }^e$as

$e_{ij}^e=u_{j,i}^e-E_{ijk}{\phi }_k^e,and{\kappa }_{ij}^e={\phi }_{j,i}^e$ (17)

herein$E_{ijk}$is the Cartesian permutations, takes a value of zero if any of the indices$ijk$is repeating, takes 1 if$ijk$is a cyclic permutation and -1 otherwise, wherein

$u_{j,i}^e\approx {\displaystyle \sum _m^N{\displaystyle \sum _n^3\frac{\partial N^m}{\partial {\xi }_n}}}\frac{\partial {\xi }_n}{\partial X_i}{\tilde{u}}_j^e,{\phi }_{j,i}^e\approx {\displaystyle \sum _m^N{\displaystyle \sum _n^3\frac{\partial N^m}{\partial {\xi }_n}}}\frac{\partial {\xi }_n}{\partial X_i}{\tilde{\phi }}_j^e$ (18)

with$X_i$being the physical coordinate in each direction$i.$ Conjugate to the Cosserat and curvature strains are the Cosserat force-stresses and the coupled stresses, respectively. In vector notation these quantities are expressed as

${\sigma }^e={[{\sigma }^e{}_{\mathrm{11,}}{\sigma }^e{}_{\mathrm{22,}}{\sigma }^e{}_{\mathrm{33,}}{\sigma }^e{}_{\mathrm{12,}}{\sigma }^e{}_{\mathrm{21,}}{\sigma }^e{}_{\mathrm{13,}}{\sigma }^e{}_{\mathrm{31,}}{\sigma }^e{}_{\mathrm{23,}}{\sigma }^e{}_{32}]}^T,$ (19)

${\mu }^e={[{\mu }^e{}_{\mathrm{11,}}{\mu }^e{}_{\mathrm{22,}}{\mu }^e{}_{\mathrm{33,}}{\mu }^e{}_{\mathrm{12,}}{\mu }^e{}_{\mathrm{21,}}{\mu }^e{}_{\mathrm{13,}}{\mu }^e{}_{\mathrm{31,}}{\mu }^e{}_{\mathrm{23,}}{\mu }^e{}_{32}]}^T.$ (20)

where${\sigma }_{ij}^e$and${\mu }_{ij}^e$are computed in each of the material phase as follows

In phase 1

${\sigma }_{ij}^e=2\left(\frac{\lambda }{2}+\frac{\mu }{d}\right){\epsilon }_{kk}^e{\delta }_{ij}+{\mu }_c\left(u_{j,i}^e-u_{i,j}^e\right)-2{\mu }_c{E}_{ijk}{\phi }_k^e$ (21)

${\mu }_{ij}^e=\{\begin{array}{l}\overline{\lambda }{\kappa }_{kk}^e{\delta }_{ij}+(\overline{\mu }-{\overline{\mu }}_c)\left({\kappa }_{ij}^e+{\kappa }_{ji}^e\right)+2\left({\mu }_{\circ }+\frac{\mu }{{\beta }^2}\right){\kappa }_{ij}^{e}if\overline{\mu }\ge {\overline{\mu }}_c,\\ \overline{\lambda }{\kappa }_{kk}^e{\delta }_{ij}-(\overline{\mu }-{\overline{\mu }}_c)\left({\kappa }_{ij}^e-{\kappa }_{ji}^e\right)+2\left({\mu }_{\circ }+\frac{\mu }{{\beta }^2}\right){\kappa }_{ij}^{e}if\overline{\mu }<{\overline{\mu }}_c.\end{array}$ (22)

In phase 2

${\sigma }_{ij}^e=\{\begin{array}{l}\overline{\lambda }{\epsilon }_{kk}^e{\delta }_{ij}+(\mu +{\mu }_c)u_{j,i}^e+(\mu -{\mu }_c)u_{i,j}^e-2{\mu }_c{E}_{ijk}{\phi }_k^e\\ -4\alpha {\beta }^2\left({({\kappa }_{ij}^e)}^2-{\beta }^2{({\epsilon }_{ij}^e-\frac{1}{d}{\epsilon }_{kk}^e{\delta }_{ij})}^2\right)\left({\epsilon }_{ij}^e-\frac{1}{d}{\epsilon }_{kk}^e{\delta }_{ij}\right)\end{array}$ (23)

${\mu }_{ij}^e=\overline{\lambda }{\kappa }_{kk}^e{\delta }_{ij}+(\overline{\mu }+{\overline{\mu }}_c){\kappa }_{ij}^e+(\overline{\mu }-{\overline{\mu }}_c){\kappa }_{ji}^e+4\alpha \left({({\kappa }_{ij}^e)}^2-{\beta }^2{({\epsilon }_{ij}^e-\frac{1}{d}{\epsilon }_{kk}^e{\delta }_{ij})}^2\right){\kappa }_{ij}^e$ (24)

In phase 3

${\sigma }_{ij}^e=\lambda {\epsilon }_{kk}^e{\delta }_{ij}+(\mu +{\mu }_c)u_{j,i}^e+(\mu -{\mu }_c)u_{i,j}^e-2{\mu }_c{E}_{ijk}{\phi }_k^e+2{\mu }_{\circ }{\beta }^2\left({\epsilon }_{ij}^e-\frac{1}{d}{\epsilon }_{kk}^e{\delta }_{ij}\right)$ (25)

${\mu }_{ij}^e=\{\begin{array}{l}\overline{\lambda }{\kappa }_{kk}^e{\delta }_{ij}+(\overline{\mu }-{\overline{\mu }}_c)\left({\kappa }_{ij}^e+{\kappa }_{ji}^e\right)\text{if}\overline{\mu }\ge {\overline{\mu }}_c,\\ \overline{\lambda }{\kappa }_{kk}^e{\delta }_{ij}-(\overline{\mu }-{\overline{\mu }}_c)\left({\kappa }_{ij}^e-{\kappa }_{ji}^e\right)\text{if}\overline{\mu }<{\overline{\mu }}_c.\end{array}$ (26)

Hence, using the representations of field vectors and tensor in matrix notation one can rewrite the weak-balance equations (13) and (14) in the following form

${\displaystyle {\int }_{\Omega }\delta {\tilde{u}}^{e,T}(B_u^T{\sigma }^e)}d\Omega -{\displaystyle {\int }_{\partial {\Omega }_u}\delta {\tilde{u}}^{e,T}(N_u^T{t}_u^e)d{\Gamma }_u=0.}$ (27)

${\displaystyle {\int }_{\Omega }\delta {\tilde{\phi }}^{e,T}(-H^T{\sigma }^e+B_{\phi }^T{\mu }^e)}d\Omega -{\displaystyle {\int }_{\partial {\Omega }_{\phi }}\delta {\tilde{\phi }}^{e,T}(N_{\phi }^T{t}_{\phi }^e)d{\Gamma }_{\phi }=0.}$ (28)

Above the matrices$B_u^{},H$and$B_{\phi }^{}$are computed according to the following

$B_u=D_u{N}_u,B_{\phi }=D_{\phi }{N}_{\phi }$ (29)

and the transformation matrix

$H\left({\xi }_1,{\xi }_2,{\xi }_3\right)=\left[H_1\left({\xi }_1,{\xi }_2,{\xi }_3\right)H_2\left({\xi }_1,{\xi }_2,{\xi }_3\right)H_3\left({\xi }_1,{\xi }_2,{\xi }_3\right)\mathrm{....}{H}_{NN}\left({\xi }_1,{\xi }_2,{\xi }_3\right)\right],$ (30)

where

$H_m^T\left({\xi }_1,{\xi }_2,{\xi }_3\right)=\left[\begin{array}{l}0000000N^m-N^m\\ 00000-N^m{N}^{m}00\\ 000N^m-N^{m}0000\end{array}\right]$ (31)

The matrices of nodal shape function$N_{\phi },N_{\phi }$and the differential operator matrices$D_u,D_{\phi }$are given by

$N_u\left({\xi }_1,{\xi }_2,{\xi }_3\right)=N_{\phi }\left({\xi }_1,{\xi }_2,{\xi }_3\right)=\left[N^{1}I{N}^{2}I{N}^{3}I\mathrm{....}{N}^{N}I\right],I={\delta }_{ij}{e}_i\otimes e_j;i,j=\mathrm{1,2,3.}$ (32)

and

$D_u=D_{\phi }={\left[\begin{array}{l}{\partial }_{1}000{\partial }_{2}0{\partial }_{3}00\\ 0{\partial }_{2}0{\partial }_{1}0000{\partial }_3\\ 00{\partial }_{3}00{\partial }_{1}0{\partial }_{2}0\end{array}\right]}^T;{\partial }_i=\frac{\partial }{\partial X_i},i=\left\{\mathrm{1,2,3}\right\}$ (33)

By using the principle of variational calculus and knowing the fact that$\delta {\tilde{u}}^{e,T}$and$\delta {\tilde{\phi }}^{e,T}$are arbitrarily chosen vectors one arrives at the following set of residual equations

${\Re }_{lin}^e={\displaystyle {\int }_{\Omega }{B}_u^T{\sigma }^e}d\Omega -{\displaystyle {\int }_{\partial {\Omega }_u}{N}_u^T{t}_u^{e}d{\Gamma }_u=\mathrm{0,}}$ (34)

${\Re }_{ang}^e={\displaystyle {\int }_{\Omega }(-H^T{\sigma }^e+B_{\phi }^T{\mu }^e)}d\Omega -{\displaystyle {\int }_{\partial {\Omega }_{\phi }}{N}_{\phi }^T{t}_{\phi }^{e}d{\Gamma }_{\phi }=0.}$ (35)

Thus the solution variables at the elemental level are obtained by solving these two equations. In order to do so the above set of nonlinear equations are linearized with respect to the elemental nodal vectors ${\tilde{u}}^e$and${\tilde{\phi }}^e,$i.e.

$\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{u}}^e}\delta {\tilde{u}}^e+\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{\phi }}^e}\delta {\tilde{\phi }}^e=\mathrm{0,}$ (36)

$\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{u}}^e}\delta {\tilde{u}}^e+\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{\phi }}^e}\delta {\tilde{\phi }}^e=0.$ (37)

The terms$\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{u}}^e},\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{\phi }}^e},\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{u}}^e}$and$\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{\phi }}^e}$here contributes to the tangent modulus operator as follows

$\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{u}}^e}=K_{uu}^e={\displaystyle {\int }_{\Omega }{B}_u^T{K}_{ee}{B}_{u}d\Omega ,}$ (38)

$\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{u}}^e}=K_{\phi \phi }^e={\displaystyle {\int }_{\Omega }\left(2{\mu }_c{H}^{T}H-H^T{K}_{ek}{B}_{\phi }+B_{\phi }^T{K}_{kk}{B}_{\phi }\right)d\Omega ,}$ (39)

$\frac{\partial {\Re }_{lin}^e}{\partial {\tilde{\phi }}^e}=K_{u\phi }^e={\displaystyle {\int }_{\Omega }{B}_u^T\left(-2{\mu }_{c}H+K_{ek}{B}_{\phi }\right)d\Omega ,}$ (40)

$\frac{\partial {\Re }_{ang}^e}{\partial {\tilde{\phi }}^e}=K_{\phi u}^e={\displaystyle {\int }_{\Omega }\left(-H^T{K}_{ee}+B_{\phi }^T{K}_{ke}\right)B_{u}d\Omega ,}$ (41)

The matrices$K_{ee},K_{ek},K_{ke}$and$K_{kk}$are to be computed from the fourth order tensors${\mathbb{K}}_{ee},{\mathbb{K}}_{ek},{\mathbb{K}}_{ke}$and${\mathbb{K}}_{kk},$respectively. These are defined according to the formulae

${\mathbb{K}}_{ee}=\frac{{\partial }^2{W}^{rel}}{\partial \nabla u\otimes \partial \nabla u},{\mathbb{K}}_{ek}=\frac{{\partial }^2{W}^{rel}}{\partial \nabla \phi \otimes \partial \nabla u},{\mathbb{K}}_{ke}=\frac{{\partial }^2{W}^{rel}}{\partial \nabla u\otimes \partial \nabla \phi },{\mathbb{K}}_{kk}=\frac{{\partial }^2{W}^{rel}}{\partial \nabla \phi \otimes \partial \nabla \phi }.$ (42)

The resulting expressions for these fourth order tensors are computed in Phase 1 as

${\mathbb{K}}_{ee}=2\left(\frac{\lambda }{2}+\frac{\mu }{d}\right)I\otimes I\text{+}{\mu }_c\mathbb{I}-{\mu }_c\overline{\mathbb{I}},{\mathbb{K}}_{ek}\text{=}{\mathbb{K}}_{ke}=\mathbb{O},$ (43)

${\mathbb{K}}_{kk}=\overline{\lambda }(I\otimes I\text{)+}\left(\overline{\mu }-{\overline{\mu }}_c+2\left({\mu }_{\circ }+\frac{\mu }{{\beta }^2}\right)\right)\mathbb{I}+\left(\overline{\mu }-{\overline{\mu }}_c\right)\overline{\mathbb{I}},$ (44)

in Phase 2 as

$\begin{array}{l}{\mathbb{K}}_{ee}=\{\begin{array}{l}\left(\lambda +2\alpha {\beta }^2\left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)\right)I\otimes I+\left(\mu +{\mu }_c-2\alpha {\beta }^2\left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)\right)\mathbb{I}\\ +\left(\mu -{\mu }_c-2\alpha {\beta }^2\left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)\right)\overline{\mathbb{I}}+8\alpha {\beta }^4\left(\text{dev}\epsilon \otimes \text{dev}\epsilon \right)\end{array},\\ {\mathbb{K}}_{ek}=-8\alpha {\beta }^2\left(\text{dev}\epsilon \otimes \kappa \right),{\mathbb{K}}_{ke}=-8\alpha {\beta }^2\left(\kappa \otimes \text{dev}\epsilon \right),\end{array}$ (45)

${\mathbb{K}}_{kk}=\overline{\lambda }I\otimes I+\left\{\mu +{\mu }_c+4\alpha \left({\Vert \kappa \Vert }^2-{\beta }^2{\Vert \text{dev}\epsilon \Vert }^2\right)\right\}\mathbb{I}+\left(\mu -{\mu }_c\right)\overline{\mathbb{I}}+8\alpha \left(\kappa \otimes \kappa \right),$ (46)

in Phase 3 as

${\mathbb{K}}_{ee}=\left(\lambda -{\mu }_{\circ }{\beta }^2\right)I\otimes I+\left(\mu +{\mu }_c+{\mu }_{\circ }{\beta }^2\right)\mathbb{I}+\left(\mu -{\mu }_c+{\mu }_{\circ }{\beta }^2\right)\overline{\mathbb{I}},$ (47)

${\mathbb{K}}_{ek}\text{=}{\mathbb{K}}_{ke}=\mathbb{O},{\mathbb{K}}_{kk}=\overline{\lambda }(I\otimes I\text{)+}\left(\overline{\mu }-{\overline{\mu }}_c\right)\mathbb{I}+\left(\overline{\mu }-{\overline{\mu }}_c\right)\overline{\mathbb{I}},$ (48)

where$\mathbb{O}$is the fourth order zero tensor,$I$is the second order identity tensor,$\mathbb{I}$ and$\overline{\mathbb{I}}$are the fourth order unit tensors computed as

$\mathbb{I}={\delta }_{ik}{\delta }_{ji}{e}_i\otimes e_j\otimes e_k\otimes e_l,\overline{\mathbb{I}}={\delta }_{il}{\delta }_{jk}{e}_i\otimes e_j\otimes e_k\otimes e_l;i,j,k,l=\mathrm{1,2,3.}$ (49)

The system of nonlinear equations to solve for each element is thus stated as

$\{\begin{array}{l}{K}_{uu}^e{u}^e+K_{u\phi }^e{\phi }^e-f_u^e\text{=}\text{0}\text{,}\\ K_{\phi u}^e{u}^e+K_{\phi \phi }^e{\phi }^e-f_{\phi }^e\text{=}\text{0}\text{.}\end{array}$ (50)

where$f_u^e$and$f_{\phi }^e$are respectively the vectors of external forces and couples acting on the body${\Omega }_{}^e$and are given by

$f_u^e={\displaystyle {\int }_{\partial {\Omega }_u}[N_u^T{t}_u^e]d{\Gamma }_u},f_{\phi }^e={\displaystyle {\int }_{\partial {\Omega }_{\phi }}[N_{\phi }^T{t}_{\phi }^e]d{\Gamma }_{\phi }}$ (51)

The global unknown vector for the displacement and micro-rotational degrees of freedoms are

$u={\left(u,\phi \right)}^T,\text{here}\text{u}\text{=}\underset{e=1}{\overset{n_e}{A}}{\text{u}}^e\phi \text{=}\underset{e=1}{\overset{n_e}{A}}{\phi }^e$ (52)

where$\underset{e=1}{\overset{n_e}{A}}$is the finite element assembly operator. The global system of nonlinear equations to solve thus can be written as

$\{\begin{array}{l}{\pi }_{lin}(u,\phi )=K_{uu}u+K_{u\phi }\phi -f_u=\mathrm{0,}\\ {\pi }_{ang}(u,\phi )=K_{\phi u}u+K_{\phi \phi }\phi -f_{\phi }=0.\end{array}$ (53)

where the global stiffness matrices

$K_{uu}\text{=}\underset{e=1}{\overset{n_e}{A}}{K}_{uu}^e\text{,}{K}_{u\phi }\text{=}\underset{e=1}{\overset{n_e}{A}}{K}_{u\phi }^e,K_{\phi u}\text{=}\underset{e=1}{\overset{n_e}{A}}{K}_{\phi u}^e,K_{\phi \phi }\text{=}\underset{e=1}{\overset{n_e}{A}}{K}_{\phi \phi }^e$ (54)

and the global external force and moment vectors are defined, respectively, as

$f_u\text{=}\underset{e=1}{\overset{n_e}{A}}{f}_u^e\text{,}{f}_{\phi }\text{=}\underset{e=1}{\overset{n_e}{A}}{f}_{\phi }^e\text{.}$ (55)

The computational algorithm for solving the set of non-linear system of equations in (53) is summarized in the Table 1.

Table 1
Solution Algorithm.

Step 1	Define material constants$E,\nu ,\lambda ,\mu ,\overline{\mu },\overline{\lambda },{\mu }_{\circ }=\min (\mu ,\overline{\mu }),{\mu }_c,{\overline{\mu }}_c,\alpha$and$\beta .$
Step 2	For each quasi-time increment: Divide the loading into tnts time steps as$t_{n+1}=t_n+\Delta t_{n+1}$where$n\in \left\{\mathrm{0,1,2,3,...,}tnts\right\}.$
Step 3	Initialize the field varialbles$u^i=\mathrm{0,}{\phi }^i=0.$
Step 4	Compute Cosserat strains$\nabla u^i$and coupled strain$\nabla {\phi }^i$
Step 5	Compute the term$\theta \left(\nabla u^i,\nabla {\phi }^i\right)={\Vert \nabla {\phi }^i\Vert }^2-{\beta }^2{\Vert \text{dev}\text{sym}\nabla u^i\Vert }^2$that characterize microstructural and non-microstructural material phase
Step 6	Calculate the material tensors${\mathbb{K}}_{{}_{ee}}^i,{\mathbb{K}}_{{}_{ek}}^i,{\mathbb{K}}_{{}_{ke}}^i$ and ${\mathbb{K}}_{{}_{kk}}^i$such that if $\theta \left(\nabla u^i,\nabla {\phi }^i\right)\ge \frac{\mu }{2\alpha {\beta }^2}$ then use equations (43) and (44), if$-\frac{{\mu }_{\circ }}{2\alpha }\le \theta \left(\nabla u^i,\nabla {\phi }^i\right)\le \frac{\mu }{2\alpha {\beta }^2}$then use equations (45) and (46), if$\theta \left(\nabla u^i,\nabla {\phi }^i\right)\le -\frac{{\mu }_{\circ }}{2\alpha }$then use equations (47) and (48).
Step 7	Determine the stiffness matrices$K_{uu},K_{u\phi },K_{\phi u}$and$K_{\phi \phi }$from equations (38), (40), (41) and (39) respectively.
Step 8	Update the solution ${\text{u}}_{n+1}^{i+1}{\text{=}\text{u}}_{n+1}^i-{\left(j_{n+1}^i\left({\text{u}}_{n+1}^i\right)\right)}^{-1}{\text{R}}_{n+1}^i\left({\text{u}}_{n+1}^i\right)$ (56) where ${\text{u}}_{n+1}^i\text{=}\left[\begin{array}{l}{\text{u}}_{n+1}^i\\ {\phi }_{n+1}^i\end{array}\right],{\text{R}}_{n+1}^i\left({\text{u}}_{n+1}^i\right)\text{=}\left[\begin{array}{l}{\pi }_{lin}^i\left({\text{u}}_{n+1}^i,{\phi }_{n+1}^i\right)\\ {\pi }_{ang}^i\left({\text{u}}_{n+1}^i,{\phi }_{n+1}^i\right)\end{array}\right]$ (57) and ${\text{J}}_{n+1}^i\left({\text{u}}_{n+1}^i\right)\text{=}{\left[\begin{array}{l}{K}_{uu}{K}_{u\phi }\\ K_{\phi u}{K}_{\phi \phi }\end{array}\right]}_{n+1}^i$ (58)
Step 9	Check convergence, if$\left|{\pi }_{lin}^{}\left({\text{u}}_{n+1}^{i+1}\right)\right|<tol$ane$\left|{\pi }_{ang}^{}\left({\text{u}}_{n+1}^{i+1}\right)\right|<tol$then Update the load increment to time$t_n\leftarrow t_{n+1}$ and go to Step 2 else Update the iteration counter$i\leftarrow i+1$and go to Step 3

5 RESULTS AND DISCUSSION

Here we present numerical experiments to demonstrate on the important features of the exact quasi-convex energy envelopes. The computations are performed to simulate granular material behavior in a Couette shear cell, under compression and in indentation. The geometry of the model is reduced to two dimensions thereby allowing computing three degrees of freedom at each point of the material domain. Two of them are the displacements$u_1$and$u_2$and third is the micro rotation${\phi }_3$. In the computations these degrees of freedom are approximated on each node of an element${\Omega }^e$using bi-quadratic interpolation functions. A plain strain assumption is used in all the three cases in consideration.

5.1 Extended microstructure in a Couette shear cell

Couette shear cell have been used to analyze the shear flows in granular materials in a number of numerical and experimental studies. For an overview on the comparison between numerical and experimental results obtained in a granular Couette shear the reader is referred to the paper by Latzel et al. (2003). These investigations provide clear evidence on the formation of different microstructural patterns in granular materials under shear deformation. These microstructural patterns appears as a result of localized deformations near the rotating cylinder in a Couette annular geometry Debregeas et al. (2001); Howell et al. (1999); Savage and Sayed (1984); Utter and Behringer (2009); Veje et al. (1999). These investigations have shown that under intense shearing different deformation patterns develops. The formation of these deformation patterns depends upon the interactions of the granular particles at micro scale. The kinematic of these particle rotations is shown in 1(c) where two possible kinds of particle rotations, namely, counter rotations and identical rotations can be seen. The micromotions of counter rotating particles contribute strong rotational effect to the proposed interaction energy potential in 5, whereas, in a situation of identical rotations of the particles it contributes much sliding effect to the interaction energy potential. Moreover, within a confined geometry and under certain boundary conditions the phenomenon of particle rotations is affected by particle size. Here, we are able to perform a numerical experiment to show that the microstructural patterns that appear in shearing granular materials can be classified and characterized as either microstructure due to translational motions of the particles or due to micromotions of the particles.

Figure 1
(a) Geometry of the two circular rotating cylinders, (b) Reduced Couette geometry and boundary conditions, (c) Kinematics of particle rotations: A schematic of rotating particle chain exhibiting two different rotational phenomenons.

Figure 2
(a) Phase with microstructure in micro-motions and no internal structure phase coexists, (b) All the three phases coexists, (c) Three phases coexists with more pronounced microstructure in micro-motions, (d) Coexistence of the three phases with almost vanishing microstructure in translational motions.

For this purpose a Couette annular geometry is taken into consideration to observe the formation of microstructure using the proposed theory. The granular material is confined between two concentric rigid circular cylinders as shown in Figure 1(a). The cylinders are subjected to rotations in opposite directions. Due to symmetry it was taken only first quadrant of an annular domain for the observation of microstructure formation in it. The annular domain is subjected to an in-plane shear deformation with the application of rotational motions at the outer boundaries. The width of the annular is taken to be 20 mm. The inner circular boundary is at a radius of 10 mm from the origin of the annulus. The circular boundaries are supposed to rotate in opposite direction and as a consequence two types of rotations among the particles inside the annular domain can be observed. These identical and counter rotations of the granular particles within the annular domain can be seen as in Figure 1(c). The boundary conditions for the numerical simulation using the proposed model uses fixed displacement along the circular boundaries, whereas, a small micromotion is prescribed at the boundaries.

The intention with this study is to observe the development of microstructural phases within the annular domain subjected to rotational deformation. The microstructure develops in both the translational and micro-rotational motions of the particles as shown in Figure 2. Where, in Figure 2(a) material exhibits a microstructure in micro-rotations of the continuum particles. In Figures 2(a), 2(b) and 2(c) it is shown that for a particular values as listed in Table 2 the material has microstructure in both the translational and micro-rotational motions of the particles. Thus allow to observe that all the material phases can coexist. Moreover, the deformed configurations with different values of$\beta$depict that decreasing the value of$\beta$causes the material to behave softly. Also, with the decrease in the particle size the microstructure in the micromotions of the particles is more pronounced.

Table 2
Material parameters for the shear test in a Couette geometry.

Figure	$E$	$\nu$	${\mu }_c$	$\alpha$	$\beta$	$\overline{\mu }$	${\overline{\mu }}_c$
-	(MPa)	-	(MPa)	(N.mm2)	mm-1	(N)	(N)
2(a)	2.0x102	0.3	2.0x100	2.0x105	5.8x10-1	8.0x101	5.0x101
2(b)	2.0x102	0.3	2.0x102	2.0x105	5.8x10-1	8.0x101	5.0x101
2(c)	2.0x102	0.3	2.0x102	2.0x105	4.0x10-1	8.0x101	5.0x101

5.2 Localized deformations in granular materials

The localization of deformations leads to possible material failure and such localized zones have been observed both numerically deBorst (1991); Trinch and Hackl (2014) and experimentally Kaus and Podladchikov (2006); Tejchman and Gudehus (1993) in a number of physical situations. The emphasis of this study is to show that it is possible to observe this phenomenon with the application of exact relaxed potentials. The exact relaxed potentials are enable enough to predict on the formation of microstructures within these localized zones in the material. Here, an investigation is presented on the formation of microstructure and thus provides a clue that how localized deformations develop that produce different possible bands. This may be important to know that under certain boundary conditions the material can form different localized deformation bands which depending upon the loading conditions afterward possibly leads to material failure. To illustrate on the formation of these localized deformation bands a tension-compression test performed on a rectangular specimen is presented. The formation of microstructural zones in the specimen clearly predicts the localized deformation mechanism observed by Kaus and Podladchikov (2006).

5.2.1 A rectangular specimen in compression

In this example, a rectangular specimen of a granular material is considered with a small imperfection in the form of a weak element at the center of the specimen as shown in Figure 3. The material parameters used for the simulation are given in Table 3. The geometry and boundary conditions are shown in Figure 3, where the vertical displacements on both top and bottom of the specimen are constrained. The material points are allowed to move horizontally at both the top and bottom boundary of the specimen except the point at the left lower corner of the specimen which is fixed in both the horizontal and vertical direction. Additionally a frictional boundary condition is used where the micro-rotation of the continuum points is allowed at both the top and bottom boundary of the specimen. A maximum displacement of 34.8 cm is applied on top boundary in vertically downward direction over a 1000 loading steps with a load step size of 4.35x10^-3.

Figure 3
(a) Geometry and boundary conditions of the rectangular specimen with weak element. (b) Selected lines along the width of the microstructural zone.

Figure 4
Rectangular specimen in compression. In first column: The deformed configuration of the specimen under compression with coarse mesh consisting of 765 elements. In second column: The deformed configuration of the specimen under compression with fine mesh consisting of 4214 elements.

The purpose of this analysis is to observe the development of localized zones related to material microstructure in the specimen under compression. Two different mesh sizes for discretization of the specimen into finite elements are used in the analysis. In the first analysis (see first column of Figure 4) the specimen is discretized into 765 finite elements whereas in the second (see second column of Figure 4) 4214 elements are used to discretize the domain of the specimen. The formation of microstructure in the material is triggered with an introduced inhomogeneity in the form of a weak element. It is observed that the microstructure in the material develops in zones where material failure may possibly occur. The developed microstructural zones in the material resembles to the localized deformation zones observed by Kaus and Podladchikov (2006). This development of the microstructure is gradually increasing with the increase in the loading. The two colors in Figure 4 depict the microstructural and non-microstructural zones of material. The red color zone is corresponding to the material phase where there is a microstructure in the material. This microstructure is due to the translational motions of the continuum particles. The onset of the localized zones is predicted by the formation of microstructure in the material as shown in Figure 4. The development of this microstructural zones leads to the information on the possible material failure. It is observed that the width of the microstructural band is not affected by the mesh size; this is highly due to the properties of the exact relaxed potentials. The microstructural bands are the regimes in material where dissipation occurs while significant deformation takes place and the maximum compressive stress exceeds the yield strength of the material.

Figure 5
Distribution of shear Cosserat rotation ${\phi }_3$(a) and (b), Cosserat strains (c), curvature strains (d), force-stress (e) and couple stress (f) components along the width of microstructural zones.

Figure 6
Geometry of the granular medium under indentation along with the prescribed boundary conditions.

Figure 5(a) and 5(b) clearly shows an increased value of the microrotation at the center across the width of the microstructural zones. The microrotations are more concentrated towards the center of the microstructural zones. This is why these microstructural zones can be seen as patterns on which the material can lead to possible failure, since inside the localized zones of deformations the particle rotations are found to be large. Since the material model used for these observations do not accommodate plastic deformations therefore this study cannot show but is able to predict on the possible material failure zones that are related to energy dissipation.

From Figure 5(c) and 5(d) it can be seen that the Cosserat strains and curvature strains are more concentrated towards the center of these microstructural zones. Figure 5(e) and 5(f) respectively shows that the couple stress and the Cosserat shear stress switch their direction at the center across the width of the microstructural zones. The physical significance of this phenomenon can be realized by considering the center line of the microstructural zones as a slip or shearing line where there is strong shearing effect which causes the flip of shear and couple stress direction. These results are in accordance to the observation of Alshibli et al. (2006) where they show the formation of strain localization in a rectangular specimen of granular material by applying Cosserat continuum theory.

Table 3
Material parameters for the specimen with introduced imperfection in compression.

-	$E$	$\nu$	${\mu }_c$	$\alpha$	$\beta$	$\overline{\mu }$	${\overline{\mu }}_c$
-	(MPa)	-	(MPa)	(N.mm2)	mm-1	(N)	(N)
Mesh	2.0x105	0.3	2.0x101	5.0x101	1.5	7.0	2.0x101
Weak element	2.0x103	0.3	2.0x101	1.0x103	1.5	3.0x102	4.0

5.2.2 Indentation test

The significance of observing the mechanical response of a granular medium subjected to indentation is evident from the load bearing capacity problems in geotechnical engineering. Particle rotations have always played an important role in load bearing capacity problems and therefore of keen interest to many researchers Aranson and Tsimring (2009); Bardet (1994); Tordesillas et al. (2005) where Cosserat continuum was used to describe the behavior of granular medium subjected to indentation.

Figure 7
Microstructure development beneath the indenter in a granular foundation. In first column: The deformed configuration of the specimen with coarse mesh consisting of 640 elements. In second column: The deformed configuration of the specimen with fine mesh consisting of 2560 elements.

In this example the microstructure formation is studied in a granular medium under indentation with the application of relaxation theory within the framework of Cosserat continuum. A plain strain assumption is used in this analysis. The granular medium of dimensions 200x100 cm² is subjected to indentation by a flat rigid indenter with a dimension of 50x5 cm² as shown in Figure 6. The geometry of the granular medium is discretized into 2560 finite elements whereas the geometry of the indenter is discretized into 250 finite elements. The indenter can only move in vertical direction and this constraint is applied by fixing the horizontal degrees of freedom of all the nodal points of indenter. Both the horizontal and vertical degrees of freedom on the right and left boundary of the granular medium are fixed. The continuum points can move only in the horizontal direction at the base of the granular medium which is ensured by fixing the vertical degrees of freedom at the base of the granular medium. The punching of the indenter is controlled by the applied vertical displacements, where a maximum displacement of 3.76 cm is applied at the top nodes of the indenter mesh in 1390 loading steps with a step size of 1.4x10^-3. The material parameters used for the indenter and the granular medium are shown in Table 4.

A large number of experiments have been performed on the granular foundations subjected indentation revealing similar bands of localized deformations as shown in this investigation. Also numerical simulation using finite element scheme for the Cosserat continuum by Walsh and Tordesillas (2006) has shown such kind of microstructure formation in a granular medium punched by rigid indenter. Here, a numerical solution is shown where the development of microstructure has been predicted in the localized zones around the indenter. The nucleation and the evolution of microstructural zone can be observed as the indenter moves downward into the material domain. The microstructural zone developed beneath the indentor is corresponding to the material regime where there are localized displacement gradients. These are the regions where the dissipation in the material takes place upon yielding in plasticity analysis. Results from the numerical simulations in Figure 7 are in accordance to the generalized Prandtl's solution of a rigid flat punch problem, where the dead material (there is no microstructure in this region) underneath the rigid indenter has a triangular shape as seen in the second column of Figure 7(h) and is in agreement with the experimental investigation Walsh and Tordesillas (2006).

Table 4
Material parameters for the indentation test on a granular medium.

-	$E$	$\nu$	${\mu }_c$	$\alpha$	$\beta$	$\overline{\mu }$	${\overline{\mu }}_c$
-	(N/cm2)	-	(N/cm2)	(N.cm2)	cm-1	(N)	(N)
Granular medium	2.0x104	0.3	2.0	5.0x104	0.5	7.0x103	2.0x102
Indenter	2.0x1012	0.3	2.0	5.0x103	0.5	7.0x103	2.0x102

The red color zones of the material are corresponding to the phase where there is a microstructure due to translational motions of the particles, whereas the purple colored zones are the regimes where there is no internal structure in the material.

6 CONCLUSION

In this article, a promising scope of the exact relaxed energy potentials is shown in modeling the granular materials and observing the extended and localized deformation microstructures. On the basis of an exact relaxed energy minimization principle a two field variational formulation is presented and a finite element method is developed to compute the deformation behavior of granular materials with microstructrues. To demonstrate on the possible applications of the proposed method different numerical examples are simulated. It is worth mentioning that these results show that the localized deformations and extended microstructures formation are not bifurcation phenomena, but arise as different aspects of the same variational model using exact quasiconvex energy envelope in the framework of generalized elasticity.

Acknowledgment

The author gratefully acknowledges financial support by the HEC-DAAD Pakistan Germany joint research cooperation scheme under the research grant number A/07/97693.

References

• Alsaleh M.I., Voyiadjis G.Z., Alshibli K.A., (2006) Modeling strain localization in granular materials using micropolar theory: Mathematical formulations. Int. J. Numer. Anal. Meth. Goemech. 30: 1501-1524.
• Alshibli K.A., Alsaleh M.I., Voyiadjis G.Z., (2006) Modelling strain localization in granular materials using micropolar theory: Numerical implementation and verification. Int. J. Numer. Anal. Meth. Goemech. 30: 1525-1544.
• Altenbach J., Altenbach H., Eremeyev V.A., (2012) On generalized Cosserat-type theories of plates and shells: a short review and bibliography. Arch. Apl. Mech. 80: 73-92.
• Aranson I., Tsimring L., (2009) Granular Patterns. Oxford University Press
• Ball J.M., (1976) Convexity conditions and existence theorems in nonlinear elasticity. Archive for rational mechanics and analysis. 63:337-403.
• Ball J.M., James R.D., (1987) Fine phase mixtures as minimizers of energy. Archive for Rational Mechanics and Analysis 100: 13-52.
• Bardet J.P., (1994) Observation on the effects of particle rotations on the failure of idealized granular materials. Mech. Mater. 8: 159-182.
• Carstensen C., Conti S., Orlando A., (2008) Mixed analytical-numerical relaxation in finite single-slip crystal plastictiy. Continuum Mech. Thermodyn 20: 275-301.
• Carstensen C., Hackl K., Mielke A., (2002) Nonconvex potentials and microstructures in finite-strain plasticity. Proc. R. Soc. Lond. A 458: 299-317.
• Conti S., Dolzmann G., Klust C., (2009) Relaxation of a class of variational models in crystal plasticity. Proc. R. Soc. A. 465: 1735-1742.
• Conti S., Hauret P., Ortiz M., (2007) Conurrent multiscale computing of deformation microstructure by relaxation and local enrichment with application to single-crystal plasticity. Multiscale Model. Simul. 6: 135-157.
• Conti S., Ortiz M., (2005) Dislocation microstructures and the effective behavior of single crystals. Arch. Rational Mech. Anal.176: 103-147.
• Conti S., Theil F., (2005) Single-slip elastoplastic microstructures. Arch. Rational Mech. Anal. 178: 125-148.
• deBorst R., (1991) Simulation of strain localization: a reappraisal of the Cosserat-continuum. Engng. comp. 8: 317-332.
• Debregeas G., Tabuteau H., di Meglio J.M., (2001) Deformation and flow of a two-dimensional foam under continuous shear. Phys. Rev. Lett. 87: 178-305.
• DeSimone A., Dolzmann G., (2002) Macroscopic Response of Nematic Elastomers via Relaxation of a Class of SO(3)-Invariant Energies. Arch. Rational mech. Anal. 161: 181-204.
• Dret H., Raoult A., (1995) The quasiconvex envelope of the Saint Venant-Kirchhoff stored energy function. Proc. Roy. Soc. Edinburgh 125A: 1179-1192.
• Forest S., Sievert R., (2006) Nonlinear microstrain theories. International Jounal of Solids and Structures 43: 7224-7245.
• Gajo A., Bigoni D., Wood D., (2004) Multiple shear band development and related instabilities in granular materials. Journal of the Mechanics and Physics of Solids 52: 2683-2724.
• Govindjee S., Hackl K., Heinen R., (2007) An upper bound to the free energy of mixing by twin-compatible lamination forn-variant martensitic phase transformations. Continuum Mech. Therm. 18: 443-453.
• Gudehus G., Nübel K., (2004) Evolution of shear bands in sand. G´eotechnique 54: 187-201.
• Gürses E., Miehe C., (2011) On evolving deformation microstructures in non-convex partially damaged solids. Journal of the Mechanics and Physics of Solids 59: 1268-1290.
• Hackl K., Heinen R., (2008) An upper bound to the free energy of n-variant polycrystalline shape memory alloys. J. Mech. Phys. Solids. 56: 2832-2843.
• Howell D., Behringer R.P., Veje C.T., (1999) Stress Fluctuations in a 2D Granular Couette Experiment: A Continuous Transition Phy. Rev. Lett., 82(26): 5241-5244.
• Kaus J.P.K, Podladchikov Y.Y., (2006) Initiation of localized shear zones in viscoelastic rocks. Journal of Geophysical Research. 111(B04412): 1-18.
• Khan M.S., Hackl K., (2012) Prediction of microstructure in a Cosserat continuum using relaxed energies, Proceedings in Applied Mathematics and Mechanics 12 (1): 265-266.
• Khan M.S., Hackl K., (2013) A constitutive model for granular materials with microstructures using the concept of energy relaxation, Proceedings in Applied Mathematics and Mechanics 13 (1): 187-188.
• Khan M.S., Hackl K., (2018) Modeling of Microstructures in a Cosserat Continuum Using Relaxed Energies. In: Rocca E., Stefanelli U., Truskinovsky L., Visintin A., (eds) Trends in Applications of Mathematics to Mechanics. Springer INdAM Series,vol 27. Springer, Cham.
• Kafadar C.B., Eringen A.C., (1976) Polar field theories. In: Eringen, A.C. (ed.) Continuum Physics, vol. IV, pp. 1-75. Academic Press, New York.
• Kohn R.V., (1991) The relaxation of a double-well energy. Cont. Mech. Thermodynam. 3: 193-236.
• Kohn R.V., Strang G., (1986a) Optimal design and relaxation of variational problems I. Comm. Pure Appl. Math. 39: 113-137.
• Kohn R.V., Strang G., (1986b) Optimal design and relaxation of variational problems II. Comm. Pure Appl. Math. 39: 139-182.
• Kohn R.V., Vogelius M., (1987) Relaxation of a variational method for impedance computed tomography. Comm. Pure Appl.Math. 40: 745-777.
• Lambrecht M., Miehe C., Dettmar J., (2003) Energy relaxation of non-convex incremental stress potentials in a strain-softening elastic-plastic bar. Internaltional Journal of Soids and Structures 40: 1369-1391.
• Latzel et al., (2003) Comparing simulation and experiment of a 2D granular Couette shear device. Eur. Phys. J. E 11: 325-333.
• Miehe C., Lambrecht M., Gürses E., (2004) Analysis of material instabilities in inelastic solids by incremental energy minimization and relxation methods: evolving deformation microstructures in finite plasticity. Journal of the Mechanics and Physics of Solids 52: 2725-2769.
• Oda M., Kazama H., (1998) Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils. G ´eotechnique 48: 465-481.
• Raoult A., (2010) Quasiconvex Envelopes in Nonlinear Elasticity, published in: Poly-,Quasi- and Rank-one Convexity in Applied Mechanics, pp. 17-51, J. Schröder, P. Neff (ed.), Springer Vienna.
• Savage S.B., Sayed M., (1984) Stresses developed by dry cohesionless granular materials sheared in an annular shear cell. J. Fluid Mech., 142: 391-430.
• Sawada K., Zhang F., Yashima A., (2006) Rotation of granular material in laboratory tests and its numerical simulation using TIJ-Cosserat continuum theory. Computational Methods, 1701-1706.
• Schaefer H., (1967) Das CosseratKontinuum. ZAMM 47: 485-498.
• Steinmann P., (1994) A micropolar theory of finite deformation and finite rotations multiplicative elasto-plasticity. International Journal of Solids and Structures 31: 1063-1084.
• Sulem J., Derrolaza M., (2002) Finite element analysis of the indentation test on rocks with microstructure. Computers and Geotechnics 29: 95-117.
• Tejchman J., Bauer E., (1996) Numerical simulation of shear band formation with a polar hypoplastic model. Comp. Geotech. 19: 221-244.
• Tejchman T., Gudehus G., (1993) Silo-music and sio-quake experiments and a numerical Cosserat approach. Powder Technol. 76: 201-212.
• Tordesillas A., Peters J.F., Gardiner B., (2004) Shear band evolution and accumulated microstructural development in Cosserat media. Int. J. Numer. Anal. Methods Geotech. Eng. 29: 981-1010.
• Tordesillas A., Peters J.F., Muthuswamy M., (2005) Role of particle rotations and rolling resistance in a semi-infinite particulate solid indented by a rigid flat punch. ANZIAM Journal 46: C260-275.
• Trinch T.B., Hackl K., (2014) Modelling of shear localization in solids by means of energy relaxation. Asia Pacific Journal onComputational Engineering 1-9.
• Utter B., Behringer R.P., (2009) Multiscale Motion in the Shear Band of Granular Couette Flow, AIP Conf. Proc. 1145: 339-342.
• Veje C.T., Daniel W., Howell W., (1999) Kinematics of a two-dimensional granular Couette experiment at the transition to shearing. Phys. Rev. E 59: 739-745.
• Walsh S.D.C., Tordesillas A., (2006) Finite element methods for micropolar models of granular materials. Applied Mathematical Modelling, 30(10): 1043-1055.

• Available online: March 25, 2019.

Publication Dates

• Publication in this collection
  25 Apr 2019
• Date of issue
  2019

History

• Received
  06 Dec 2018
• Reviewed
  30 Jan 2019
• Accepted
  13 Mar 2019