Micromechanical modelling of the anisotropic creep behaviour of granular medium as a fourth-order fabric tensor

The main objective of the present work is to develop a micromechanics approach to predict the macroscopic anisotropic creep behaviour of granular media. To this end, the linear viscoelastic behaviour of the inter-particle interaction at contact is adopted, and the contact distribution is characterized by a fourth-order fabric tensor in the local scale. Then, fourth-order tensor fabric-based micromechanical approaches based on Voigt and Reuss localization assumptions are applied to granular media in the Laplace–Carson space. With help of the inverse Laplace–Carson transformation of these obtained models, the macroscopic anisotropic creep behaviour of granular media submitted to a constant external loading is examined. Finally, the obtained results by specializing the Burgers model into the obtained models are compared with the numerical simulations in the particle flow code (PFC2D) to illustrate the validation and the accuracy of the analytical models for the macroscopic anisotropic creep behaviour of granular media.


Introduction
Granular media are commonly used in engineering structures, such as rockfill dams, foundations, highway embankments, etc. They play significant roles in durability and service life of these relevant practical engineering feats. 1 Thus, investigating the mechanical behaviour of granular media has attracted a lot of research attention in many fields. 2,3 Granular media are typical composites 4 and their overall mechanical behaviours are mainly dependent on their internal characteristics. 5 At the grain scale, granular media are defined as an assembly of small discrete particles whose interactions are transmitted from one particle to another during contact. 6 Generally, the particles in the granular assembly are assumed to be rigid or semi-rigid, so that the inter-particle interaction can be considered as point forces acting at contacts. 7 Additionally, the inter-particle contact distribution in granular assembly exhibits a distinct anisotropic feature that is generally expressed by a fabric tensor. 8 The particle interaction behaviour and the anisotropic contact distribution greatly influence the macroscopic characteristics of granular media. 9 Therefore, to predict the macroscopic properties of granular media, the 1 effects of the contact behaviour and the contact contribution at the local scale must be taken into account. 10 A great amount of work has been dedicated to predicting the mechanical behaviour of granular media. 11,12 For example, in the framework of linear elasticity, Walton 13 predicted the elastic modulus of granular media without considering the local anisotropic contact distribution effect. Millet et al. 14 proposed a micromechanical approach for the behaviour of granular media in which the anisotropic contact distribution was taken into account by means of a fourth-order fabric tensor. In the framework of viscoelasticity, some approaches have been proposed to simulate the creep behaviour of granular media. For example, Misra and Singh 12 applied the average method to obtain the integral form of the macroscopic creep property of granular media and took the numerical method to illustrate its macroscopic creep behaviour. In these works, the local anisotropic contact distribution effect was taken into account. 15 However, as mentioned, the local contact anisotropic effect plays an important role on the macroscale behaviour of granular media. Thus, in the present work, we propose a micromechanical approach to predict the creep behaviour of granular media in which a fourth-order fabric tensor is introduced to simulate the local contact anisotropic distribution.
The proposed approach has two significant features. The first considers the closed-form creep behaviour of granular media, which is obtained by combining the granular micromechanical approach based on Voigt and Reuss localization assumptions and classical Laplace-Carson (LC) transformation. The second considers local anisotropic contact distribution by introducing a fourth-order fabric tensor. Additionally, the comparison between the simulation results of the particle-flow code (PFC2D) and the results of obtained models used as examples of granular media submitted to constant external loading are illustrated via the validation and accuracy of the analytical models. The present work can be viewed as an expansion of Millet et al. 14 regarding the setting of linear viscoelasticity.
The paper is organized as follows. In the next section, the local and global viscoelastic constitutive equations associated to inter-particle interactions and granular media, the Laplace-Carson transformation and its inverse are presented. In Section 3, the micromechanical derivation of the viscoelastic behaviour of granular media is presented in detail. In Sections 4 and 5, the specialization of obtained models is taken with the Burgers model, and a numerical example is provided to verify the obtained models. Finally, some conclusions are drawn in Section 6.

LC transformation and its inverse
In the following, the LC transformation and its inverse formulation are introduced. 16 For a continuous function, f (t), defined at time t, the LC transform, L½ f (t) or e f (s) has been defined by 17 : where s is the variable in the LC space and L½ is the LC transform operator. The inverse LC transform of e f (s), expressed by From the definition of the LC transform and its inverse transform, we observe that the LC transform is just its Laplace transform multiplied by s and the inverse LC transform is just its Laplace transform divided by s in the integral formula. The domain of integration of the LC transform is the same as the Laplace transform. Thus, it is obvious that the LC transform and its inverse transform have the same property as the Laplace transformation and its inverse transform. With the help of LC, the linear viscoelastic problem described in real space can be transformed into the LC space, which is equal to the linear elastic problem in its form. Next, we consider two functions: g(t) and h(t). Their corresponding convolution function, f (t), is expressed by and its LC, L½f (t), can be written by where~is the convolution symbol.

Linear viscoelastic interface model
In the problem under consideration, the interface at the contact point, c between two contact particles at the microscale is supposed to be linear viscoelastic and is expressed by where K is the second-order tensor stiffness, f c (t) is the force exerted on the contact surface and u c (t) is the relative displacements between the particles. Applying LC to the above equation, we have It is obvious that the above equation in LC space is equal to the well-known elastic spring-type interface model in form. Additionally, the second-order creep tensor function is decomposed into normal and tangent parts: where n c and t c denote the unit normal and tangent vectors at a contact, and k n (t) and k t (t) are the normal and tangent scalar creep functions, respectively. Thus, equation (3) is rewritten by where f u c n and f u c t are the normal and tangent displacement components, and f k c n and f k c t are the corresponding terms of the normal and tangent scalar-creep functions in LC space, respectively.

Constitutive equation for linear viscoelasticity
At the macroscopic scale, the constitutive stress-strain relation for a viscoelastic material is given classically by a Stieljes integral as with 18 : where S hom is the fourth-order tensor-creep function.
With help of the LC transform, we have The above equation can also be written as with g C hom = g S hom À1 . We observe that equations (6) and (7) in the LC space are mathematically identical to Hooke's law in form.

Micromechanical derivation of the granular viscoelastic behaviour of
In the previous section, the constitutive equations defined at the micro and macro scales were presented. In the following, we derive the closed-form of the twodimensional viscoelastic behaviour of granular by the homogenization procedure based on Voigt and Reuss localization assumptions, 19,20 respectively. The main feature of our derivation shows that the anisotropy effects caused by the contact distribution on the macroscopic viscoelastic behaviour can be handled by introducing a fourth-order fabric tensor.

Basic equations
Considering a representative volume element (RVE), V , of granular materials (see Figure 1), its average stress, S, on time t is defined by where sym represents the symmetry of the corresponding tensor, and the summation, c, is on the contacts with c = 1, Á Á Á , N. A is the total area occupied by RVE. The branch vector l c joins the centres of two contacting particles. f c is the force exerted between two particles. Furthermore, it should be noted that the contact force, f c , is linear viscoelastic and expressed by equation (2). Additionally, to simplify the problem, we assume that the granular assembly is comprised of circular grains, and the branch vector can be expressed as where r denotes the average radius of the grains, and n c denotes the unit normal vector at a contact.
In the case of a great number of contacts, it is preferable to introduce the probability of contact, P(n), in the direction of the vector, n. Thus, equation (8) reduces to where P(n) is the probability of contact in the direction of the normal, n, and s, the unit circle. With help of LC, the above equation can be transformed into a formulation in the LC space: where e f c is expressed by equation (4).

Description of the fourth-order fabric tensor
For considering the anisotropic effect caused by the contact distribution at the micro scale, 21 a fourth-order fabric tensor R is introduced and defined by where hi denotes the average over all the orientations (over a unit circle). It should be noted that the contact distribution of granular material at the micro scale is usually termed as its Fabric and is introduced for describing its microstructural feature. With help of the probability of contact, P(n), the fourth-order fabric tensor, R, can be rewritten as Referring to, 14 the probability of contact, P(n), can be expressed by P(n)=16 R :: (n n n n) À 3 4 D : (n n) where D is the second-order tensor, which is defined by In this subsection, the second-and fourth-order fabric tensors are defined by the distribution of contact orientations. With these fabric tensors, the anisotropy effect caused by the contact distribution at the micro scale is introduced into the macro-scale creep behaviour of granular materials, which is presented in detailed in the following subsection.

Voigt kinematic assumption
Taking the Voigt kinematic assumption 22,23 where the strain field in the representative volume element (RVE) is uniform, we write the relationship between the local relative displacement, u c , and the macroscopic strain, E, as where E is considered to be uniform and imposed on the boundary of RVE. With the help of equations (16), (2) and (10) is expressed by Thus, equation (11) can be given by where the third-order tensor, T c , defined by T c = (I À n c n c ) n c , which is introduced to simplify our notations. T cT denotes the transposed T c . By identifying equations (7) and (18), we derive where g C hom denotes the macroscopic stiffness tensor in LC space.
Considering the definition of the fourth-order fabric tensor, R, and the second-order fabric tensor, D, we can obtain 24 for arbitrary symmetry second-order tensors A and B. to simplify the problem, we make the following assumption: e k t = a e k n . Thus, introducing the two above equations into the equation (19), we get the expression, g C hom It is obvious that the anisotropy property is involved in the macroscopic stiffness tensor, g C hom , in LC because of the fabric tensors, D and R. To this point, the LC inverse can be applied to obtain the closed form of the creep behaviour of granular materials where the anisotropic property is taken into account in the framework of the Voigt kinematic assumption.

Reuss static localization assumption
When it comes to the Reuss static localization assumption where the stress field in the representative volume element (RVE) is uniform, we take the following expression: where the symmetric two-order tensor, B, is a localization operator. The substitution of equation (11) into equation (10) gives us With the definition of the second-order fabric tensor, D, of equation (15), the above equation can be rewritten as According to equation (22), the localization operator, B, can analytically be expressed by Next, capitalizing on the virtual work principle, we can derive a consistent definition of the average strain tensor involving relative displacement at the point of contact, u c , and the localization operator, B, in LC space: With the probability of contact, P(n), the above equation can be rewritten in integral form: Similarly, we introduce the inverse formulation of equation (4): where f f c n is the normal contact force component, and e f c t is the tangent contact force vector. f h c n and f h c t are the corresponding terms of the normal and tangent scalar creep functions in LC space, respectively.
Introducing equation (24) into equation (23), we have According to equation (21), we get Introducing equation (26) into equation (25) leads to By identification with equation (6) and considering the definitions of B, D and R, we get the explicit expression of the homogenized compliance tensor in LC: Assuming e h t = 1 a e h n , we obtain with the definition of the symbol, by A B ð Þ ijkl = A ik B jl . It is clear that the final homogenized stiffness and compliance tensors, g C hom and g S hom , are general and can be adopted to specific viscoelastic interface models, such as Maxwell, Kelvin and Burgers. In addition, it should be noted that with help of Minimum potential energy principle and complementary energy principle, the Voight and Reuss localization assumptions lead to the Voight upper and Reuss lower bounds for the effective elastic moduli of composites, respectively.

Specialization of the viscoelastic behaviour to the Burger's model
In this section, we particularize the obtained homogenized tensors, g C hom and g S hom , into the Burgers viscoelastic interface model (see Figure 2). The normal and tangent scalar interface creep functions, k n (t) and k t (t) or h c n (t) and h c t T ð Þ, are described by the Burger's model and expressed by 26 where independent material constants b n1 , b n2 and b t1 , b t2 are related to two elasticites, whereas independent material constants h n1 ,h t1 and h n2 , h t2 are related to two viscosities for the normal and tangent creep functions, respectively. With the help of LC, according to equation (28), we can express e k n and e k t as The substitution of equation (29) into equation (20) leads to Similarly, g S hom ijpq associated to the Reuss assumption can expressed by 1 À 1 a D À1 :R:D À1 Note that we adapt the inverse LC transformation to obtain the closed-form of the creep behaviour of the granular media for equations (30) and (31).

Numerical discussions
In this section, we consider a numerical sample of granular material subject to a constant unaxial loading to verify the validity of the analytical model by comparing the analytical results and the PFC2D numerical results.
The sample of granular material under consideration is taken as a two-dimensional rectangle having a height and width of 0.6 m 3 0.3 m established by PFC2D. The radius of generated circle granular particles of this sample ranges from 0.005 to 0.02 m. For considering numerical full contact, a prior uniaxial pressure, S 0 yy = 0:1kPa, is firstly applied to this granular sample. Finally, the numerical granular sample is given, which contains 333 circle particles and 815 contacts (see Figure 3). The associated contact distribution, or, equivalently, the probability of contact in the direction of n i = ( cos u i , sin u i ) from 0 to 2p for the numerical granular sample is numerically expressed by where DN i denotes the number of contacts in the direction, n i contained in the discretization interval ½u i À Du 2 , u i + Du 2 of 0, 2p ½ and N denotes the total number of contacts. In this study, we take Du = 10 0 or p=18, so that the corresponding DN i in all directions, n i , is given as shown in the following Table 1.
The obtained results clearly indicate the anisotropy of the granular material as plotted in Figure 4.
The corresponding fabric tensors are obtained by  On the other hand, the distribution of contact orientations can be calculated by the fourth-order fabric tensors. 27,28 Note that equation (14) can be further deduced as P(n) = 1 + 16R 1111 cos 4 u + 64R 1112 cos 3 u sin u + 96R 1122 sin 2 u cos 2 u + 64R 1222 sin 3 u cos u + 16R 2222 sin 4 u À 12D 11 cos 2 u À 24D 12 sin u cos u À 12D 22 sin 2 u The distribution of contact orientations from the fourth-order fabric tensor for this sample are shown in the same Figure 4. It is observed that the use of a fourth-order fabric tensor seems to qualitatively and quantitatively represent the anisotropic contact distribution of the numerical sample.
In the following, the obtained sample of granular material was subjected to a constant uniaxial compressive loading for numerically illustrating our results. The corresponding material properties are listed in Table 2.
With these numerical values, under the constant uniaxial compressive loading of (S yy = 0:1MPa, S xx = S xy = 0), the analytical uniaxial compressive strain, E yy , based on Voigt and Reuss assumptions are given by, respectively, The numerical and analytical results from PFC2D are plotted in Figure 5. The curves based on Voigt and Reuss formulations are the upper and the lower bounds, and, as expected, the numerical curve is situated between these two derived bounds.

Conclusion
The objective of this paper was to derive the closed form of the anisotropic creep behaviour of the granular media in the framework of a two-dimensional problem. For linear viscoelastic granular media, the stress-strain relationship was nonlinear in real space. Thus, Laplace transform was used to cause the nonlinear relationship in real space to become a linear problem in Laplace space, which is formally equivalent to the problem of elasticity. Then, the creep constitutive equation was derived through a series of linear operations in Laplace space. Furthermore, a fourth-order fabric tensor was used to reflect the contact distribution in the local scale. Finally, we took the Burgers model as a special   It is clear that the fourth order fabric tensor can well describe a general anisotropy viscoelastic granular material. According to the analytical solutions of equations (20) and (27), we know the interfacial properties of grain were very important to the creep behaviour. In other words, the evolution trend of creep is dependent on the interfacial properties. We got two analytical solutions about creep behaviour of granular media based on different assumption respectively, and then we verified the accuracy of analytical solutions based on Burgers model.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.