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

Abstract

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. ¹ 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 ⁴ and their overall mechanical behaviours are mainly dependent on their internal characteristics. ⁵ 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. ⁶ 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. ⁷ Additionally, the inter-particle contact distribution in granular assembly exhibits a distinct anisotropic feature that is generally expressed by a fabric tensor. ⁸ The particle interaction behaviour and the anisotropic contact distribution greatly influence the macroscopic characteristics of granular media. ⁹ Therefore, to predict the macroscopic properties of granular media, the effects of the contact behaviour and the contact contribution at the local scale must be taken into account. ¹⁰

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 ¹³ predicted the elastic modulus of granular media without considering the local anisotropic contact distribution effect. Millet et al. ¹⁴ 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 ¹² 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. ¹⁵ 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. ¹⁴ 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.

Basic equations for linear viscoelasticity

LC transformation and its inverse

In the following, the LC transformation and its inverse formulation are introduced. ¹⁶ For a continuous function, $f(t)$ , defined at time $t$ , the LC transform, $L[f(t)]$ or $\tilde{f}(s)$ has been defined by ¹⁷ :

$${\tilde{f}}_c(s)=L[f(t)]=s{\int }_0^{\infty }{e}^{-\mathrm{st}}f(t)\mathrm{dt}$$

where $s$ is the variable in the LC space and $L[]$ is the LC transform operator. The inverse LC transform of $\tilde{f}(s)$ , expressed by

$$f(t)=L^{-1}[{\tilde{f}}_c(s)]=\frac{1}{2\pi i}{\int }_{\alpha -i\infty }^{\alpha +i\infty }\frac{1}{s}{e}^{\mathrm{st}}{\tilde{f}}_c(s)\mathrm{ds}$$

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

$$f(t)={\int }_0^{t}g(t-\tau )\frac{\mathrm{dh}\left(\tau \right)}{d\tau }d\tau =g⊛h$$

(1)

and its LC, $L[f(t)]$ , can be written by

$$L[f(t)]=g\left(s\right)h\left(s\right)$$

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

$${\mathit{f}}^c(t)={\int }_0^t\mathit{K}(t-\tau ).\frac{d{\mathit{u}}^c\left(\tau \right)}{d\tau }d\tau =\mathit{K}(t)⊛{\mathit{u}}^c(t)$$

(2)

where $\mathit{K}$ is the second-order tensor stiffness, ${\mathit{f}}^c(t)$ is the force exerted on the contact surface and ${\mathit{u}}^c(t)$ is the relative displacements between the particles. Applying LC to the above equation, we have

$$L[{\mathit{f}}^c(t)]={\tilde{\mathit{f}}}^c=\tilde{\mathit{K}}\left(s\right).\widetilde{{\mathit{u}}^c}\left(s\right)$$

(3)

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:

$$\mathit{K}(t)=k_n(t){\mathit{n}}^c\otimes {\mathit{n}}^c+k_t(t){\mathit{t}}^c\otimes {\mathit{t}}^c$$

where ${\mathit{n}}^c$ and ${\mathit{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

$$\widetilde{{\mathit{f}}^c}=\widetilde{k_n}\widetilde{u_n^c}{\mathit{n}}^c+\widetilde{k_t}\widetilde{u_t^c}{\mathit{t}}^c$$

(4)

where $\widetilde{u_n^c}$ and $\widetilde{u_t^c}$ are the normal and tangent displacement components, and $\widetilde{k_n^c}$ and $\widetilde{k_t^c}$ 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 ¹⁸ :

$$\mathit{E}(t)={\int }_0^t{\mathbb{S}}^{\mathrm{hom}}(t-\tau ):\frac{d\mathrm{\Sigma }\left(\tau \right)}{d\tau }d\tau ={\mathbb{S}}^{\mathrm{hom}}⊛\mathrm{\Sigma }$$

(5)

where ${\mathbb{S}}^{\mathrm{hom}}$ is the fourth-order tensor-creep function. With help of the LC transform, we have

$$\tilde{\mathit{E}}\left(s\right)=\widetilde{{\mathbb{S}}^{\mathrm{hom}}}\left(s\right):\tilde{\mathrm{\Sigma }}\left(s\right)$$

(6)

The above equation can also be written as

$$\tilde{\mathrm{\Sigma }}=\widetilde{{\mathbb{C}}^{\mathrm{hom}}}:\tilde{\mathit{E}}$$

(7)

with $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}={\widetilde{{\mathbb{S}}^{\mathrm{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 two-dimensional 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, $\mathrm{\Sigma }$ , on time $t$ is defined by

$$\mathrm{\Sigma }(t)=\frac{1}{A}\sum _{c=1}^N\frac{1}{2}({\mathit{f}}^c\otimes {\mathit{l}}^c+{\mathit{l}}^c\otimes {\mathit{f}}^c)=\frac{1}{A}\sum _{c=1}^N{({\mathit{f}}^c\otimes {\mathit{l}}^c)}^{\mathrm{sym}}$$

(8)

where $\mathrm{sym}$ represents the symmetry of the corresponding tensor, and the summation, $c$ , is on the contacts with $c=1,\cdots ,N$ . $A$ is the total area occupied by RVE. The branch vector ${\mathit{l}}^c$ joins the centres of two contacting particles. ${\mathit{f}}^c$ is the force exerted between two particles. Furthermore, it should be noted that the contact force, ${\mathit{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

$${\mathit{l}}^c=2r{\mathit{n}}^c$$

(9)

where $r$ denotes the average radius of the grains, and ${\mathit{n}}^c$ denotes the unit normal vector at a contact.

Figure 1.

Contact geometry of particles $p$ and $q$ .

In the case of a great number of contacts, it is preferable to introduce the probability of contact, $P(\mathit{n})$ , in the direction of the vector, $\mathit{n}$ . Thus, equation (8) reduces to

$$\mathrm{\Sigma }(t)=\frac{\mathrm{Nr}}{\pi A}[{\int }_{S}P(\mathit{n}){({\mathit{f}}^c\otimes {\mathit{n}}^c)}^{\mathrm{sym}}\mathrm{dS}]$$

(10)

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:

$$\tilde{\mathrm{\Sigma }}\left(s\right)=\frac{\mathrm{Nr}}{\pi A}[{\int }_{S}P(\mathit{n}){(\widetilde{{\mathit{f}}^c}\otimes {\mathit{n}}^c)}^{\mathrm{sym}}\mathrm{dS}]$$

(11)

where $\widetilde{{\mathit{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, ²¹ a fourth-order fabric tensor $\mathbb{R}$ is introduced and defined by

$$\mathbb{R}=〈\mathit{n}\otimes \mathit{n}\otimes \mathit{n}\otimes \mathit{n}〉=\frac{1}{N}\sum _{k=1}^N{\mathit{n}}^k\otimes {\mathit{n}}^k\otimes {\mathit{n}}^k\otimes {\mathit{n}}^k$$

(12)

where 〈〉 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(\mathit{n})$ , the fourth-order fabric tensor, $\mathbb{R}$ , can be rewritten as

$$\mathbb{R}=\frac{1}{2\pi }{\int }_{S}P\left(\mathit{n}\right)(\mathit{n}\otimes \mathit{n}\otimes \mathit{n}\otimes \mathit{n})\mathrm{dS}$$

(13)

Referring to, ¹⁴ the probability of contact, $P(\mathit{n})$ , can be expressed by

$$P(\mathit{n})=16[\mathbb{R}::(\mathit{n}\otimes \mathit{n}\otimes \mathit{n}\otimes \mathit{n})-\frac{3}{4}\mathit{D}:(\mathit{n}\otimes \mathit{n})+\frac{1}{16}]$$

(14)

where $\mathit{D}$ is the second-order tensor, which is defined by

$$\mathit{D}=〈\mathit{n}\otimes \mathit{n}〉=\frac{1}{2\pi }{\int }_{S}P(\mathit{n})(\mathit{n}\otimes \mathit{n})\mathrm{dS}$$

(15)

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, ${\mathit{u}}^c$ , and the macroscopic strain, $\mathit{E}$ , as

$${\mathit{u}}^c=\mathit{E}.{\mathit{l}}^c=2r\mathit{E}.{\mathit{n}}^c$$

(16)

where $\mathit{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

$$\mathrm{\Sigma }(t)=\frac{2N{r}^2}{\pi A}[{\int }_{S}P(\mathit{n}){(\mathit{K}(t)⊛\mathit{E}.{\mathit{n}}^c\otimes {\mathit{n}}^c)}^{\mathrm{sym}}\mathrm{dS}]$$

(17)

Thus, equation (11) can be given by

$$\begin{array}{c}\tilde{\mathrm{\Sigma }}=\frac{2N{r}^2}{\pi A}{\int }_{s}P\left(\mathit{n}\right)[{\tilde{k}}_n({\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c)\hfill \\ +{\tilde{k}}_t({\mathit{T}}^{\mathrm{cT}}·{\mathit{T}}^c)]\mathrm{dS}:\tilde{\mathit{E}}\hfill \end{array}$$

(18)

where the third-order tensor, ${\mathit{T}}^c$ , defined by ${\mathit{T}}^c=(\mathit{I}-{\mathit{n}}^c\otimes {\mathit{n}}^c)\otimes {\mathit{n}}^c$ , which is introduced to simplify our notations. ${\mathit{T}}^{\mathrm{cT}}$ denotes the transposed ${\mathit{T}}^c$ .

By identifying equations (7) and (18), we derive

$$\begin{array}{c}\widetilde{{\mathbb{C}}^{\mathrm{hom}}}=\frac{2N{r}^2}{\pi A}{\int }_{s}P\left(\mathit{n}\right)[\widetilde{k_n}({\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c)\hfill \\ +\widetilde{k_t}({\mathit{T}}^{\mathrm{cT}}·{\mathit{T}}^c)]\mathrm{ds}\hfill \end{array}$$

(19)

where $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}$ denotes the macroscopic stiffness tensor in LC space.

Considering the definition of the fourth-order fabric tensor, $\mathbb{R}$ , and the second-order fabric tensor, $\mathit{D}$ , we can obtain^24,25

$$\frac{1}{2\pi }{\int }_{s}P\left(\mathit{n}\right)({\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c)\mathrm{ds}=\mathbb{R}$$

$$\frac{1}{2\pi }{\int }_{s}P\left(\mathit{n}\right)({\mathit{T}}^{\mathrm{cT}}·{\mathit{T}}^c)\mathrm{ds}=[{\left(\mathit{D}\overline{\underset{¯}{\otimes }}\mathbf{\delta }\right)}^{\mathrm{sym}}-\mathbb{R}]$$

where $\overline{\underset{¯}{\otimes }}$ denotes ${(\mathit{A}\otimes \mathit{B})}_{\mathrm{ijkl}}=\frac{1}{2}(A_{\mathrm{ij}}{B}_{\mathrm{kl}}+A_{\mathrm{ik}}{B}_{\mathrm{jl}})$ for arbitrary symmetry second-order tensors $\mathit{A}$ and $\mathit{B}$ . to simplify the problem, we make the following assumption: $\widetilde{k_t}=\alpha \widetilde{k_n}$ . Thus, introducing the two above equations into the equation (19), we get the expression, $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}$

$$\widetilde{{\mathbb{C}}^{\mathrm{hom}}}=\frac{4N{r}^2}{A}\widetilde{k_n}\{\mathbb{R}+\alpha [{\left(\mathit{D}\overline{\underset{¯}{\otimes }}\mathbf{\delta }\right)}^{\mathrm{sym}}-\mathbb{R}]\}$$

(20)

It is obvious that the anisotropy property is involved in the macroscopic stiffness tensor, $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}$ , in LC because of the fabric tensors, $\mathit{D}$ and $\mathbb{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:

$${\mathit{f}}^c=\mathrm{\Sigma }.\mathit{B}.{\mathit{n}}^c$$

(21)

where the symmetric two-order tensor, $\mathit{B}$ , is a localization operator. The substitution of equation (11) into equation (10) gives us

$$\mathrm{\Sigma }=\frac{2\mathrm{rN}}{A}{[(\mathrm{\Sigma }.\mathit{B}).\frac{1}{2\pi }{\int }_{S}P(\mathit{n})({\mathit{n}}^c\otimes {\mathit{n}}^c)\mathrm{ds}]}^{\mathrm{sym}}$$

With the definition of the second-order fabric tensor, $\mathit{D}$ , of equation (15), the above equation can be rewritten as

$$\mathrm{\Sigma }=\frac{2\mathrm{Nr}}{A}{(\mathrm{\Sigma }.\mathit{B}.\mathit{D})}^{\mathrm{sym}}$$

(22)

According to equation (22), the localization operator, $\mathit{B}$ , can analytically be expressed by

$$\mathit{B}=\frac{A}{2\mathrm{Nr}}{\mathit{D}}^{-1}$$

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, ${\mathit{u}}^c$ , and the localization operator, $\mathit{B}$ , in LC space:

$$\tilde{\mathrm{\Sigma }}:\tilde{\mathit{E}}=\frac{1}{A}\sum _{c=1}^N\widetilde{{\mathit{f}}^c}.\widetilde{{\mathit{u}}^c}=\frac{1}{A}\sum _{c=1}^N\tilde{\mathrm{\Sigma }}:(\mathit{B}.{\mathit{n}}^c\otimes \widetilde{{\mathit{u}}^c})$$

With the probability of contact, $P(\mathit{n})$ , the above equation can be rewritten in integral form:

$$\tilde{\mathrm{\Sigma }}:\tilde{\mathit{E}}=\tilde{\mathrm{\Sigma }}:\frac{N}{2\pi A}{\int }_{s}P\left(\mathit{n}\right)(\mathit{B}.{\mathit{n}}^c)\otimes \widetilde{{\mathit{u}}^c}\mathrm{ds}$$

which leads to

$$\tilde{\mathit{E}}=\frac{N}{2\pi A}{[{\int }_{s}P\left(n\right)(\mathit{B}.{\mathit{n}}^c)\otimes \widetilde{{\mathit{u}}^c}\mathrm{ds}]}^{\mathrm{sym}}$$

(23)

Similarly, we introduce the inverse formulation of equation (4):

$$\widetilde{{\mathit{u}}^c}=\widetilde{h_n^c}\widetilde{f_n^c}{\mathit{n}}^c+\widetilde{h_t^c}\widetilde{{\mathit{f}}_t^c}$$

(24)

where $\widetilde{f_n^c}$ is the normal contact force component, and $\widetilde{{\mathit{f}}_t^c}$ is the tangent contact force vector. $\widetilde{h_n^c}$ and $\widetilde{h_t^c}$ are the corresponding terms of the normal and tangent scalar creep functions in LC space, respectively.

Introducing equation (24) into equation (23), we have

$$\tilde{\mathit{E}}=\frac{N}{2\pi A}{[{\int }_{s}P\left(n\right)\mathit{B}.{\mathit{n}}^c\otimes (\widetilde{h_n^c}\widetilde{f_n^c}{\mathit{n}}^c+\widetilde{h_t^c}\widetilde{{\mathit{f}}_t^c})\mathrm{ds}]}^{\mathrm{sym}}$$

(25)

According to equation (21), we get

$$\begin{array}{c}\widetilde{f_n^c}=\widetilde{f^c}.{\mathit{n}}^c=(\tilde{\mathit{\Sigma }}.B):\left({\mathit{n}}^c\otimes {\mathit{n}}^c\right)\end{array}$$

$$\begin{array}{c}\widetilde{{\mathit{f}}_t^c}=\widetilde{{\mathit{f}}^c}-\widetilde{f_n^c}{\mathit{n}}^c=(\tilde{\mathit{\Sigma }}.\mathit{B}):\left(\widetilde{{\mathit{n}}^c}\otimes \mathit{\delta }-{\mathit{n}}^c\otimes {\mathit{n}}^c\otimes {\mathit{n}}^c\right)\end{array}$$

(26)

Introducing equation (26) into equation (25) leads to

$$\begin{array}{c}{\tilde{E}}_{\mathrm{ij}}=\frac{1}{2}\frac{N}{2\pi A}{B}_{\mathrm{po}}\{B_{\mathrm{ik}}{\int }_{S}P(\mathit{n})[(\widetilde{h_n^c}-{\tilde{h}}_t^c)n_o^c{n}_q^c{n}_j^c{n}_k^c\hfill \\ +{\tilde{h}}_t^c{n}_o^c{\delta }_{\mathrm{qj}}{n}_k^c]\mathrm{ds}+B_{\mathrm{jk}}{\int }_{S}P(\mathit{n})[({\tilde{h}}_n^c-{\tilde{h}}_t^c)\hfill \\ n_o^c{n}_q^c{n}_i^c{n}_k^c+{\tilde{h}}_t^c{n}_o^c{\delta }_{\mathrm{qi}}{n}_k^c]\mathrm{ds}\}{\tilde{\mathrm{\Sigma }}}_{\mathrm{pq}}\hfill \end{array}$$

By identification with equation (6) and considering the definitions of $\mathit{B}$ , $\mathit{D}$ and $\mathbb{R}$ , we get the explicit expression of the homogenized compliance tensor in LC:

$$\begin{array}{c}\widetilde{{\mathbb{S}}^{\mathrm{hom}}}=\frac{A}{4N{r}^2}[(\widetilde{h_n^c}-\widetilde{h_t^c}){\mathit{D}}^{-1}.\mathbb{R}.{\mathit{D}}^{-1}\hfill \\ +\widetilde{h_t^c}({\mathit{D}}^{-1}\underset{¯}{\otimes }\mathit{I}+\mathit{I}\underset{¯}{\otimes }{\mathit{D}}^{-1})]\hfill \end{array}$$

Assuming $\widetilde{h_t}=\frac{1}{\alpha }\widetilde{h_n}$ , we obtain

$$\begin{array}{c}\widetilde{{\mathbb{S}}^{\mathrm{hom}}}=\frac{A}{4N{r}^2}\widetilde{h_n^c}[(1-\frac{1}{\alpha }){\mathit{D}}^{-1}.\mathbb{R}.{\mathit{D}}^{-1}\hfill \\ +\frac{1}{\alpha }({\mathit{D}}^{-1}\underset{¯}{\otimes }\mathit{I}+\mathit{I}\underset{¯}{\otimes }{\mathit{D}}^{-1})]\hfill \end{array}$$

(27)

with the definition of the symbol, $\underset{¯}{\otimes }$ by ${(\mathit{A}\otimes \mathit{B})}_{\mathrm{ijkl}}=A_{\mathrm{ik}}{B}_{\mathrm{jl}}$ .

It is clear that the final homogenized stiffness and compliance tensors, $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}$ and $\widetilde{{\mathbb{S}}^{\mathrm{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, $\widetilde{{\mathbb{C}}^{\mathrm{hom}}}$ and $\widetilde{{\mathbb{S}}^{\mathrm{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_n^c(t)$ and $h_t^c\left(T\right)$ , are described by the Burger’s model and expressed by ²⁶

$$\begin{array}{c}{k}_n={\beta }_{n1}^{-1}(1-\exp (-\frac{{\beta }_{n1}}{{\eta }_{n1}}T))+{\beta }_{n2}^{-1}+{\eta }_{n2}^{-1}T\hfill \\ \multicolumn{1}{c}{k_t={\beta }_{t1}^{-1}(1-\exp (-\frac{{\beta }_{t1}}{{\eta }_{t1}}T))+{\beta }_{t2}^{-1}+{\eta }_{t2}^{-1}T}\end{array}$$

(28)

where independent material constants ${\beta }_{n1}$ , ${\beta }_{n2}$ and ${\beta }_{t1}$ , ${\beta }_{t2}$ are related to two elasticites, whereas independent material constants ${\eta }_{n1}$ , ${\eta }_{t1}$ and ${\eta }_{n2}$ , ${\eta }_{t2}$ are related to two viscosities for the normal and tangent creep functions, respectively.

Figure 2.

Burger’s model.

With the help of LC, according to equation (28), we can express $\widetilde{k_n}$ and $\widetilde{k_t}$ as

$$\begin{array}{c}\multicolumn{1}{c}{\widetilde{k_n}=1/(\frac{1}{{\beta }_{n2}}+\frac{1}{s{\eta }_{n2}}+\frac{1}{{\beta }_{n1}+s{\eta }_{n1}})}\\ \multicolumn{1}{c}{\widetilde{k_t}=1/(\frac{1}{{\beta }_{t2}}+\frac{1}{s{\eta }_{t2}}+\frac{1}{{\beta }_{t1}+s{\eta }_{t1}})}\end{array}$$

(29)

The substitution of equation (29) into equation (20) leads to

$$\begin{array}{c}\widetilde{{\mathbb{C}}^{\mathrm{hom}}}=\frac{4N{r}^2}{A(\frac{1}{{\beta }_{n2}}+\frac{1}{s{\eta }_{n2}}+\frac{1}{{\beta }_{n1}+s{\eta }_{n1}})}\hfill \\ \times \{\mathbb{R}+\alpha [{\left(D\overline{\underset{¯}{\otimes }}\delta \right)}^{\mathrm{sym}}-\mathbb{R}]\}\hfill \end{array}$$

(30)

Similarly, $\widetilde{S_{\mathrm{ijpq}}^{\mathrm{hom}}}$ associated to the Reuss assumption can expressed by

$$\widetilde{S_{\mathrm{ijpq}}^{\mathrm{hom}}}=\frac{A}{4N{r}^2}(\frac{1}{{\beta }_{n2}}+\frac{1}{s{\eta }_{n2}}+\frac{1}{{\beta }_{n1}+s{\eta }_{n1}})$$

(31)

$$\begin{array}{c}[(1-\frac{1}{\alpha }){\mathit{D}}^{-1}.\mathbb{R}.{\mathit{D}}^{-1}\hfill \\ +\frac{1}{\alpha }({\mathit{D}}^{-1}\underset{¯}{\otimes }\mathit{I}+\mathit{I}\underset{¯}{\otimes }{\mathit{D}}^{-1})]\hfill \end{array}$$

(32)

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 × 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, ${\Sigma }_{\mathrm{yy}}^0=0.1\mathrm{kPa}$ , 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 {\theta }_i,\sin {\theta }_i)$ from $0$ to $2\pi$ for the numerical granular sample is numerically expressed by

$$\frac{1}{2\pi }P(n_i)\mathrm{\Delta }\theta =\frac{\mathrm{\Delta }{N}_i}{N}$$

where $\mathrm{\Delta }{N}_i$ denotes the number of contacts in the direction, $n_i$ contained in the discretization interval $[{\theta }_i-\frac{\mathrm{\Delta }\theta }{2},{\theta }_i+\frac{\mathrm{\Delta }\theta }{2}]$ of $[0,2\pi ]$ and $N$ denotes the total number of contacts. In this study, we take $\mathrm{\Delta }\theta ={10}^0$ or $\pi /18$ , so that the corresponding $\mathrm{\Delta }{N}_i$ in all directions, $n_i$ , is given as shown in the following Table 1.

Figure 3.

Numerical sample prior to testing.

Table 1.

Contact distribution.

$\theta$	[−5, 5]	[5, 15]	[15, 25]	[25, 35]	[35, 45]	[45, 55]
$\mathrm{\Delta }N$	18	26	16	27	30	42
$\theta$	[55, 65]	[65, 75]	[75, 85]	[85, 95]	[95, 105]	[105, 115]
$\mathrm{\Delta }N$	73	83	77	65	74	75
$\theta$	[115, 125]	[125, 135]	[135, 145]	[145, 155]	[155, 165]	[165, 175]
$\mathrm{\Delta }N$	75	50	33	22	17	12

The obtained results clearly indicate the anisotropy of the granular material as plotted in Figure 4.

Figure 4.

The distribution of contact normal in different planes, prior to testing.

The corresponding fabric tensors are obtained by

$$\left[\mathit{D}\right]=\left[\begin{array}{cc}0.316\hfill & -0.031\hfill \\ \multicolumn{1}{c}{-0.031}& 0.684\hfill \end{array}\right]$$

$$\left[\mathbb{R}\right]=\left[\begin{array}{ccc}0.194\hfill & 0.123\hfill & -0.014\hfill \\ \multicolumn{1}{c}{0.123}& 0.561\hfill & -0.017\hfill \\ \multicolumn{1}{c}{-0.014}& -0.017\hfill & 0.123\hfill \end{array}\right]$$

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

$$\begin{array}{c}P(n)=1+16R_{1111}\stackrel{4}{\cos }\theta +64R_{1112}\stackrel{3}{\cos }\theta \sin \theta \hfill \\ +96R_{1122}\stackrel{2}{\sin }\theta \stackrel{2}{\cos }\theta +64R_{1222}\stackrel{3}{\sin }\theta \cos \theta \hfill \\ +16R_{2222}\stackrel{4}{\sin }\theta -12D_{11}\stackrel{2}{\cos }\theta \hfill \\ -24D_{12}\sin \theta \cos \theta -12D_{22}\stackrel{2}{\sin }\theta \hfill \end{array}$$

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.

Table 2.

Materials properties.

$\alpha$	${\beta }_{n1}$ (MN/m)	${\beta }_{n2}$ (MN/m)	${\eta }_{n1}$ (MN · h/m)	${\eta }_{n2}$ (MN · h/m)
0.5	10	10	10	10,000

With these numerical values, under the constant uniaxial compressive loading of ( ${\Sigma }_{\mathrm{yy}}=0.1\mathrm{MPa},{\Sigma }_{\mathrm{xx}}={\Sigma }_{\mathrm{xy}}=0$ ), the analytical uniaxial compressive strain, $E_{\mathrm{yy}}$ , based on Voigt and Reuss assumptions are given by, respectively,

$$\begin{array}{c}{E}_{\mathrm{yy}}^{\mathrm{Voigt}}=0.0765[\frac{1}{10}+\frac{t}{10000}+\frac{1}{10}(1-e^{-t})]\hfill \end{array}$$

$$\begin{array}{c}{E}_{\mathrm{yy}}^{\mathrm{Reuss}}=0.2165[\frac{1}{10}+\frac{t}{10000}+\frac{1}{10}(1-e^{-t})]\hfill \end{array}$$

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.

Figure 5.

Comparison of analytical and numerical solutions for the specimen.

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 example to illustrate the features of our obtained model. Consequently, conclusions are as follows:

• 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.