A 3D fourth order fabric tensor approach of anisotropy in granular media

Computational Materials Science 46 (2009) 869–880

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A 3D fourth order fabric tensor approach of anisotropy in granular media Jamila Rahmoun a, Djimédo Kondo b,*, Olivier Millet c a b c

Laboratory of Industrial and Human Automatic, Mechanic and Informatics, UMR CNRS 8530, University of Valenciennes, Le Mont Houy, 59313 Valenciennes, France Laboratoire de Mécanique de Lille, Université de Lille 1, Boulevard Paul Langevin, 59655 Villeneuve d’Ascq, France Laboratoire d’Etude des Phénomènes de Transfert Appliqués au Bâtiment, Université de La Rochelle, Avenue Michel Crépeau, 17000 La Rochelle, France

a r t i c l e

i n f o

Article history: Received 14 March 2009 Received in revised form 12 April 2009 Accepted 21 April 2009 Available online 3 June 2009 Keywords: Granular media Micromechanics Fabric tensor Anisotropy Homogenization

a b s t r a c t We propose a micromechanical approach for granular media, with a particular account of the textureinduced anisotropy and of the strain localization rule. The approach is mainly based on the consideration of a fourth order fabric tensor able to capture general anisotropy which can be induced by complex distribution of contacts. Incorporation of this fourth order fabric tensor in a suitable homogenization scheme allows to determine the corresponding macroscopic elastic properties of the granular material. For this purpose, in addition to the classical Voigt upper bound, a new kinematics-based localization rule is proposed. It generalizes the one formulated by Cambou et al. [B. Cambou, Ph. Dubujet, F. Emeriault, F. Sidoroff, Eur. J. Mech. A/Solids 14 (1995) 225–276] in the case of an isotropic contact distribution. The results of the complete model compare well to numerical simulations results when available [C.S. Chang, C.L. Liao, Appl. Mech. Rev. 47 (1 Part 2) (1994) 197–207] (case of isotropic distribution of contacts). Finally, the interest of the fourth order fabric tensor based approach combined with the proposed localization rule is shown for different distributions of contacts by comparing its predictions to those given by a second order fabric tensor approach. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Granular media are involved in many engineering applications including soils mechanics, powder technology, etc. It is well known that the macroscopic properties and behavior of such class of materials, which can be considered as a grains assembly, is the result of particle interactions mainly due to interparticles contacts. The relevant microscopic variables are strongly related to the contact forces between grains and to the relative displacement of grains. Whatever is the mechanical nature of these contacts, their geometrical network leads to a particular structure which constitutes a fundamental component of the material characteristics. It has been experimentally shown that the contact distribution can be very complex depending on the compaction mode or on the loading path applied to the material [3,4]. The internal structure of granular media corresponds mainly to their texture which can be influenced by different factors [5–10]. The first factor is the particles shape. Indeed, for non spherical grains but with prolate or oblate shape, one has a preferential orientations of grains in the assembly and consequently a directional character of the contact distribution. The methodology of preparation of a granular material sample (deposit of grains by a source

* Corresponding author. E-mail addresses: [email protected] (J. Rahmoun), djimedo. [email protected] (D. Kondo), [email protected] (O. Millet). 0927-0256/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2009.04.030

point, pluviation, etc.) and the deformation history may also strongly affect the distribution of the contacts [11–14]. Since there is an obvious link between the contact directions and the anisotropy of the granular assembly referring to its contact network,1 it is particulary interesting to investigate the possibility to predict the macroscopic anisotropic properties from the microstructure. This is precisely the main purpose of different modeling procedures including homogenization techniques (see [1,15–18]). In most of these works the anisotropy is described by means of fabric tensors representing the contact distribution. However, the existing analyses are limited either by their 2D formulation or by the restricted type of anisotropy (orthotropy) that they considered due to the use of a second order fabric tensor. The main objective of the present paper is to determine the elastic anisotropic properties of granular media by using not only a second order fabric tensor, but also a fourth order one [19]. In order to simplify the problem, we consider that at any contact between grains, the local relation between forces and relative displacements are linear. We then investigate a new kinematical rule which provides the link between microscopic strain field and macroscopic strain. Comparisons between the obtained closed form expression of the elastic stiffness and numerical results available in literature [2,20] show the interest of the proposed 1 It must be emphasized that only anisotropy due to contact distribution is addressed in the present work and that no attempt is done to model load-induced anisotropy which can evolves.

870

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

homogenization procedure based on the fabric tensors. The present paper generalizes approach recently developed by Millet et al. in a 2D context [21]. The following tensorial notations are adopted:

ða  bÞijkl ¼ aij bkl ;

1.5 1 0.5

1 ðabÞ ¼ ðab þ abÞ 2

-1

-0.5

0.5

1

-0.5

2. Description of the texture the granular material -1

In is well known that the contact probability PðnÞ in the normal direction n is the most considered characteristics to define fabric tensors of granular media [22]. Classically, the second order tensorial representation of the fabric tensor reads in the 3D context (see for instance [7,23,24]):

D¼

1 4p

Z

PðnÞðn  nÞds

ð1Þ

s2

in which s2 is the boundary of the unit sphere. PðnÞ satisfies the following normalization condition:

1 4p

Z

PðnÞds ¼ 1

ð2Þ

s2

for an anisotropic distribution of contacts, while for an isotropic one, PðnÞ ¼ 1. Following the work of Lubarda and Krajcinovic [23] (see also [24]) in the context of cracks density parameter distribution, we approximate the distribution of the contact probability PðnÞ in the normal direction n by a second order tensor d as PðnÞ ¼ d : ðn  nÞ. It can be show that PðnÞ can be expressed with respect to the macroscopic variable D. Indeed, in 3D, we have:

D¼

1 4p

Z

PðnÞðn  nÞds

and trD ¼ 1

ð3Þ

s2

from which one gets:

d¼

    15 1 15 1 D  ðtrDÞd ¼ D d 2 5 2 5

ð4Þ

Indeed:

  15 1 PðnÞ ¼ D : ðn  nÞ  2 5

ð5Þ

Z

PðnÞðn  n  n  nÞds and trðDÞ ¼ D

ð6Þ

s

Indeed, according to Lubarda and Krajcinovic [23]:

Q¼

  15 21 1 D  7D  d þ d  d 4 2 2

ð7Þ

It follows that:

PðnÞ ¼

Fig. 1. Comparsion of the different theoretical contact distributions with experimental data [4].

where trðDÞ ¼ D, that is Dijkk ¼ Dij ; similar property is valid for any other contraction of D on two indexes. It must be emphasized that the description of the contact distribution by means of a second order fabric tensor D is sufficient when the materials exhibits at most an orthotropic behavior. In order to illustrate the degree of induced anisotropy, we compare (see Fig. 1) the contact distributions obtained from a second order and fourth order fabric tensors approach with experimental data. The experiments description can be found in [4]: deposit of 722 discs of radius 13, 18 and 28 mm under the effect of gravity. From Fig. 1, it is noted that contact probability PðnÞ suggests the anisotropy induced by the vertical pile of the grains. We also represent on the same Fig. 1 the contact distributions given by (5) and (8), calculated from the definitions (1) and (6), respectively. It is noticed that both probabilities are significantly different. The fourth order fabric tensor provide results in agreement with experimental data [4], whereas the second order fabric tensor leads to a rough approximation of the induced anisotropy.

Before presenting the results provided by the Voigt assumption, let us first recalling the principle and methodology of kinematic approaches.

Let us consider now a fourth order fabric tensor in the 3D context. To this end, it is convenient to represent PðnÞ with the help of a fourth order tensor Q such as PðnÞ ¼ Q<ðn  n  n  nÞ. Similar to the above reasoning, we aim, in a first time, at expressing PðnÞ with respect to the macroscopic variable D. In 3D, one has:

1 4p

fourth_order_fabric_tensor second_order_fabric_tensor isotropic_distribution experimental_results

Z

We have then trðDÞ ¼ 13 trðdÞ ¼ 1, and trðdÞ ¼ 3. From these results, d h i D  15 d . It follows that: is given by d ¼ 15 2

D¼

Legend

3. Homogenization based Voigt kinematic assumption

Z 1 PðnÞðn  nÞds ¼ d : n  n  n  nds 4p s2 s2 h i 1 1 ¼d: ðd  d þ 2ddÞ ¼ trðdÞd þ 2d 15 15

1 D¼ 4p

-1.5

  15 21 1 D<ðn  n  n  nÞ  7D : ðn  nÞ þ 4 2 2

ð8Þ

3.1. Principle and methodology We aim at deriving by means of homogenization techniques, the macroscopic elastic behavior of the granular medium with account of anisotropy due to the contact distribution. For this purpose, a linear elastic contact law is considered. Denoting F c the contact forces between grains, a classical definition of the quasi static average Cauchy stress tensor is [25–28]:

R¼

N 1X F c  Lc V c¼1

ð9Þ

in which V denotes the volume of the representative elementary volume (rev) considered, and where the summation is performed on the N interior contacts. From (9), it appears that the macroscopic stress tensor depends on the distribution of the local contact forces. However, R as defined by (9) is not symmetric; therefore, an another definition established by a direct use of the principle of the virtual work [29] is generally adopted. Indeed:

R : DE ¼

N 1X ðF c :Duc Þ V c¼1

ð10Þ

871

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

t

in which use has been done of the fact that nc  F ct ¼ 0. Reporting this result in (15), we obtain:

c nc

R¼

c

g1

Lc

g2

N h i 2r X F cn nc  nc þ T cT  F ct V c¼1

ð16Þ

3.2. The Voigt assumption Fig. 2. Representative scheme of the contact between two grains.

where Duc denotes the virtual displacements of grains associated to a virtual macroscopic deformation E. Again, the summation is considered only on the active contacts c ¼ 1; . . . N. Considering now Eq. (10), we get: N 1X R : DE ¼ ðF c  Lc Þ : DE V c¼1

ð11Þ

P since, DE is symmetric, only the symmetric part of V1 Nc¼1 ðF c  Lc Þ contributes in (11). Consequently, a suitable definition of the average stress tensor is:

" #sym N 1 X R¼ F c  Lc V c¼1 sym

1 ða 2

ð12Þ

c

c

c

L ¼ kL kn ¼ 2rn

ð13Þ

Let us now introduce the following decomposition of the displacement uc in the local basis ðt; nÞ at the contact:

uc ¼ ucn nc þ uct uc and nc denote, respectively, the displacements vector at the contact and the unit vector normal to the contact. As already indicated a linear contact law is considered. Then, the contact forces writes:

F c ¼ F cn nc þ F t

ð14Þ

We aim now at deriving an expression of the homogenized elastic stiffness tensor Chom of the granular medium. Consideration of (12)–(14) leads to:

R¼

" N  2r X V

uci ¼ Eij Lcj

ð17Þ

where Eij are the components of the macroscopic strain tensor and Lci is the branch vector which joins the centers of grains at contact point c. Then, Eq. (13) leads to:

uci ¼ Eij Lcj ¼ 2rEij ncj ¼ 2rðE  nc Þi

ð18Þ

From (18), it is readily seen that:

ucn ¼ uc  nc ¼ 2rnc  E  nc ¼ 2rðnc  nc Þ : E uct ¼ uc  ucn nc ¼ 2rE  nc  ucn nc ¼ 2rT c : E

T

where ðaÞ ¼ þ a Þ. It must be emphasized that in the case of a great number of grains, the antisymmetrical part of the expression (9) is negligible [30], and (9) and (12) lead to the same results. For simplicity, we assume in the present study that the granular assembly is constituted of spherical grains having the same radius r. Therefore, we have (see Fig. 2): c

The kinematic assumption of Voigt, which assumes that the deformations are uniform inside the rev, provides a simple relation between the relative displacements at grains contacts and the macroscopic deformation [31,6,32]. It is recalled that the Voigt assumption reads [33,34]:

 F cn nc þ F ct  nc

#sym ð15Þ

c¼1

where the summation are carried out on the active contacts c ¼ 1; . . . ; N. On the other hand, it can be easily observed that:

ðF ct  nc Þsym ¼ T cT  F ct where T c ¼ nc  I  ðnc  nc  nc Þ. Indeed:

h i ðT cT  F ct Þij ¼ ðnc  I  nc  nc  nc ÞT  F ct ij   ¼ ðI  nc  nc  nc  nc Þ  F ct ij   ¼ ðI  nc Þ  F ct ij ¼ Iijkl ncl ðF t Þck 1 ¼ ðdik djl þ dil djk Þncl ðF t Þck 2 i 1h ¼ ncj ðF t Þci þ nci ðF t Þcj 2

c

c

c

c

ð19Þ

c

with T ¼ n  I  ðn  n  n Þ in which I denotes the symmetric fourth order unit tensor:

Iijkl ¼

1 ðdik djl þ dil djk Þ 2

Considering the contact law (14) and the kinematic assumption (18), we have:

F cn ¼ K n ucn ¼ 2rK n ðnc  nc Þ : E F ct ¼ K t uct ¼ 2rK t T c : E

ð20Þ

where K n and K t represent, respectively, the normal and the tangential rigidity at the contact. By reporting this result in (16), we get:

R¼

N h i 4r 2 X K n ðnc  nc  nc  nc Þ þ K t ðT cT :T c Þ : E V c¼1

ð21Þ

It follows that the macroscopic elastic constitutive relation R ¼ Chom : E is given by the following expression of the homogenized elastic stiffness tensor:

Chom ¼

N h i 4r 2 X K n ðnc  nc  nc  nc Þ þ K t ðT cT  T c Þ V c¼1

ð22Þ

It is now possible to illustrate the anisotropy of the elastic macroscopic behavior by expressing the macroscopic stiffness Chom in term of PðnÞ as:

Chom ¼

 Z 4Nr2 1 PðnÞðn  n  n  nÞds Kn 4 p s2 V  Z 1 þ Kt PðnÞðT T :TÞds 4 p s2

ð23Þ

This paths the way to the use of the fabric tensors of various orders introduced in Section 2 as the appropriate mathematical descriptors of the anisotropy. 3.3. Predictions based on the second order fabric tensor We take now advantage of the fact that PðnÞ can be defined by means of the second fabric tensor D in order to express Chom . To do this, we first replace the following integrals:

872

1 4p

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

Z

PðnÞðn  n  n  nÞds and

s2

1 4p

Z

4. The proposed kinematic localization rule

PðnÞðT T  TÞds

s2

by their respective expressions (B.2) and (B.3) derived in Appendix B. It comes:

Chom ¼

  4Nr 2 1 1 K n  ðd  d þ 2ddÞ þ ðD  d þ d  DÞ 7 5 V i 1  1 þ 2ðDd þ dDÞ  K t  ðd  d þ 2ddÞ 7 5   3 þ D  d þ d  DÞ  ðDd þ dDÞ 2

ð24Þ

which satisfies the usual symmetries of the elastic stiffness tensor. In the general case, Chom has the standard symmetries (21 independent coefficients). Its components, whose detailed expressions are given in Appendix D, depend on the components of D and on the contact rigidities K n and K t .

The interest of a kinematical localization rule more elaborated than the Voigt assumption for granular media, has been already demonstrated in several studies in literature (see for instance [1]). However, due to the limitation of these studies to the isotropic case, we aim here at proposing a general kinematic assumption appropriate to anisotropic granular media. As usually, this kinematic rule must link the relative displacements at the contact uc to the uniform deformation E applied at the boundary of the rev. For this purpose, we will take advantage of representation theorems as described for instance in [37]. Let us before that introduce the vector U C defined by:

Uc ¼

1 PðnÞuc 2r

ð28Þ

4.1. The case of second order method 3.4. The case of the fourth order fabric tensor By using representation theorems described in [37,38], we get: Now, let us derive the general expression (23) of Chom by using the description of PðnÞ by means of the fourth order fabric tensor D. The computation involves the following integrals:

1 4p

Z

PðnÞðn  n  n  nÞds and

s2

1 4p

Z

h b 00 þ A c0 trðDÞ þ A c1 trðEÞ þ A c2 trðE  DÞ þ A c3 E : ðn  nÞ Uc ¼ A c4 ðE  DÞ : ðn  nÞ þ A c5 trðEÞD : ðn  nÞ þA

PðnÞðT T  TÞds

c6 E : ðn  nÞD : ðn  nÞ þ A c7 trðEÞtrðDÞ þA i c8 trðDÞðn  E  nÞ þ A c9 D : ðn  nÞ n þA h i c1 D : ðn  nÞ þ B c2 trðDÞ þ B c3 þE  n B h i þ D  n C^1 E : ðn  nÞ þ C^2 trE þ C^3

s2

which are expressed in Appendix C (Eqs. (C.2) and (C.3)). It follows that:

Chom ¼

  4Nr 2 4Nr 2 1 KnD þ K t ðDd þ dDÞ  D 2 V V

ð25Þ

Obviously, the presence of D in this expression of Chom appears to be more relevant for the description of anisotropy than (24), the latter being able to describe only orthotropic macroscopic behavior.

c1 ðE  D  n  D  E  nÞ þD

ð29Þ

Since trðDÞ ¼ 1, it follows that:

h U c ¼ A0 þ A1 trðEÞ þ A2 trðE  DÞ þ A3 E : ðn  nÞ

3.5. Isotropic distribution of contacts As a first illustration, we consider an isotropic distribution of the contacts. In this case, as indicated before, PðnÞ ¼ 1 and then:

1 1 n  nds ¼ d 4p s2 3 Z 1 1 D¼ n  n  n  nds ¼ ðd  d þ 2ddÞ 4p s 2 15

D¼

þ A4 ðE  DÞ : ðn  nÞ þ A5 trðEÞD : ðn  nÞ i þ A6 E : ðn  nÞD : ðn  nÞ þ A7 D : ðn  n n h i þ E  n B1 D : ðn  nÞ þ B2 h i þ D  n C 1 E : ðn  nÞ þ C 2 trE þ C 3

Z

ð26Þ

Expression (24) reduced then to:

Chom ¼



o 4Nr 2 n

K n d  d þ 2dd  K t d  d þ 3dd 15V hom

from which the macroscopic bulk modulus k lhom read:

Chom ¼ k

lhom ¼

hom

þ D1 ðE  D  n  D  E  nÞ

and shear modulus

ðd  dÞ þ 2lhom ðddÞ

2 Nr2 K n 4 Nr 2 K n hom ¼ ð2 þ 3aÞ and k 15 V 9 V

ð27Þ

where a ¼ KKnt . This result coincides with the one provided by Ref. [2] precisely in the isotropic case (see also [31,35,36]). This verification confirms only the coherence of the present approach developed in a more general context, e.g. anisotropic. However, a quantitative comparison of this result to the numerical values obtained from discrete element simulations2 indicates that the Voigt assumption is too simple and not sufficient to predict the overall properties of the granular material. A more relevant kinematical localization rule is due.

2

See forthcoming comparison in Section 5.

ð30Þ

c0 þ A b 00 , A1 ¼ A c1 þ A c7 , A3 ¼ A c3 þ A c8 and B2 ¼ B c2 þ B c3 . with A0 ¼ A Di . For the other terms, we have Ai ¼ Abi ; Bi ¼ Bbi ; C i ¼ C^i and Di ¼ c Considering a linear elastic material without initial stresses and strains (in the natural equilibrium state), in the absence of external loading without possible rigid solid motion, we obtained a linear relation between the displacement uc and E. Therefore, for the general form (30) of U C one gets A0 ¼ A7 ¼ C 3 ¼ 0. Now, in order to establish the general kinematic assumption (Eq. (30)), let us consider the average strain, which utilize equally the macroscopic strain, written in the following form:

E¼

N 4p V

Z

sym PðnÞA  nc  uc ds

ð31Þ

s2

in which enters a second order tensor A which is introduced in Appendix E. Componentwise, one has:

" #sym  sym Z N 1 X N 1 c c Eij ¼ ðui Ajk nk Þ ¼ PðnÞðuci Ajk nck Þds V c¼1 V 4p s 2

ð32Þ

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

Taking into account (30), it is readily seen that:

Eij ¼



2rN 1 V 4p

Z s2

sym U ci Ajk nck ds

4.2. Determination of the macroscopic elastic stiffness tensor

ð33Þ

Then, by reporting expression (30) of U c in (33) and following the procedure detailed in Appendix E, it is shown that:

i 3 h 5n  E  n  trðEÞ n 2 4 h i 15 ðn  D  nÞE  n þ b n  5D : ðn  nÞn þ 2D  n þ 2  1 2 þ k1  trðEÞn  trðE  DÞn þ E : ðn  nÞD : ðn  nÞn 35 35  2 2  D : ðn  nÞE  n  E : ðn  nÞD  n 7 7  2 4 þ k2  trðEÞn  trðE  DÞn þ E : ðn  nÞn 5 5

U c ¼ gE  n 



g

þ

 1 5 þ k3  trðEÞn  trðE  DÞn þ E : ðn  nÞn 2 2

ð38Þ



 2r 3g 15 g 3 trðEÞ E : ðn  nÞ þ  þ þ PðnÞ 4 2 2 4 15 þ D : ðn  nÞE : ðn  nÞ 2   2r 15 c c ut ¼ ðu  ucn nÞ ¼ g þ D : ðn  nÞ T : E PðnÞ 2 ucn ¼ uc  n ¼

þ 5ðE  DÞ : ðn  nÞn  5D : ðn  nÞE  n þ E  D  n   1 5  D  E  n þ k4  trðEÞn  2trðE  DÞn þ E : ðn  nÞn 2 2

where T ¼ I  n  n  n  n and T : E ¼ E  n  E : ðn  n  nÞ due to

ð34Þ

It should be mentioned that in the particular isotropic case where PðnÞ ¼ 1 and D ¼ 13 d, one has:



 1 1 1 1 1 uc ¼ 2rU c ¼ 2r A0 þ A7 þ C 3 þ A1 þ A2 þ A5 þ C 2 trðEÞ 3 3 3 3 3     1 1 1 1 þ A3 þ A4 þ A6 þ C 1 ÞE : ðn  nÞ n þ B1 þ B2 E  n 3 3 3 3

ð35Þ which leads to: 820 1 > >

<6B C Bb  5 b þ 2 bC þ 3 þ 1 g  1 k1  1 k2  5 k3  5 k4 uc ¼ 2r 6 4@ > 3ffl{zfflfflfflfflfflfflfflfflffl 3 ffl}A 4 2 21 3 6 6 > |fflfflfflfflfflfflfflfflffl : 0 

15 5 5 5 25 25  g k 1  k2  k3  k4 E : ðn  nÞ n  trðEÞ þ  4 2 21 3 6 6   5 2 5 5 5 k 1 þ k 2  k3  k 4 E  n ð36Þ þ gþ  2 21 3 3 3 The searched kinematical assumption is then obtained in the following form:



5 2 2 5 5 uc ¼ 2r g þ  k1  k2  k3  k4 E  n 2 21 3 3 3

h i 3 g 1 1 5 5 þ  k1  k2  k3  k4 5n  E  n  trðEÞ n  4 2 21 3 6 6 ð37Þ One can note that the expression of uc as function of E is not affected by b. On the other hand, it is also interesting to notice that in this iso2 k1  23 k2  53 k3  53 k4 ¼  32, tropic case, for the particular value g  21 we recover the kinematic Voigt assumption (see Eq. (17)):

U ¼ 2rE  n



 i 2r g 3 h 5n  E  n  trðEÞ n gE  n  þ PðnÞ 2 4 15 D : ðn  nÞE  n þ 2

We aim now at applying the homogenization scheme based on this kinematic assumption (38). In the same way as the assumption of Voigt, we consider together the expression of the average stress tensor (12), the contact law (14) and the kinematical rule (38) in order to determine the macroscopic properties of the granular medium. For this end, let us first decompose the displacement in the local basis associated to a given contact. We have:

 þ 4ðE  DÞ : ðn  nÞn þ trðEÞD : ðn  nÞn

c

The anisotropic kinematical rule proposed here is obviously more general than the Voigt assumption. It depends on five adjustable parameters k1 ; k2 ; k3 ; k4 ; b as well as of D. For sake of simplicity, we limit here the study to the particular case of one adjustable parameter3 g by considering ki ¼ 0; i ¼ 1 to 4 and b ¼ 0. This choice allows to reduce the number of the free parameters without loss of generality regarding the description of the material anisotropy:

uc ¼

 2D : ðn  nÞE  n  2E : ðn  nÞD  n

þ 10ðE  DÞ : ðn  nÞn  5D : ðn  nÞE  n   5E : ðn  nÞD  n þ trE D  n

873

the symmetry of E. On the other hand, with the help of the linear elastic contact law, we have:



  2r 3g 15 g 3 trðEÞ E : ðn  nÞ þ Kn  þ þ PðnÞ 4 2 2 4 15 D : ðn  nÞE : ðn  nÞ þ 2

F cn ¼ K n ucn ¼

and

F ct ¼ K t uct ¼

2r PðnÞ



gþ

  15 D : ðn  nÞ K t T : E 2

By reporting these two last expressions in (12), and by using the relation ðF ct  nÞsym ¼ T T  F ct , we get:

R¼

Z   4Nr 2 1 3g 15 Kn  E : ðn  nÞ þ 4 V 4p s 2 2

 15 g 3 trðEÞ ðn  nÞ D : ðn  n  n  nÞ : E þ þ þ 2 2 4   15 T þ gþ D : ðn  nÞ K t T :T : E ds 2

The above expression of the macroscopic stress R reads also:  

Z 4r 2 N 3g 15 1 R¼ n  n  n  nds þ Kn  2 V 4 4p s2

Z g 3 1 n  n  dds þ þ 2 4 4p s2  Z 15 1 D: n  n  n  n  n  nds þ 2 4p s2   Z Z 1 15 1 T T  Tds þ D: ðn  nÞ  T T  Tds :E þ Kt g 4p s2 2 4p s2 ð39Þ

3 The detailed study of the kinematic assumption (34) has to be done in order to provide the physical significance of the parameters ki . However, this does not constitute the main purpose of the present study.

874

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

Identification with R ¼ Chom : E leads to the following expression of the homogenized elastic stiffness tensor (cf. (B.5)–(B.8) and (B.9) in Appendix B):

Chom



  4r 2 N g 5 1 2 g 1 ¼ dd Kn  ðd  dÞ þ ðddÞ þ þ þ V 5 2 4 5 6 4   1 1 ðd  d þ 2ddÞ þ ðD  d þ d  D þ 2Dd þ 2dDÞ þ 14 7     4r 2 N g 1 2 1 þ Kt  ðd  dÞ þ ðddÞ  ðd  d  5ddÞ V 5 14 3 5  1 3 3 g ð40Þ þ ðD  d þ d  D  Dd  dDÞ þ dd 7 2 2 3

Note that for g ¼  32, the macroscopic stiffness tensor provided by the Voigt assumption is recovered. 4.3. Anisotropic modeling by means of the fourth order fabric tensor The starting point is here, the relation between PðnÞ and the fourth order fabric tensor D (Eq. 8). By analogy with (38) given in the case of a second order fabric tensor, the kinematic assumption reads:

i g 3 h U c ¼ gE  n  þ 5n  E  n  trðEÞ n þ aD<ðn  n  n  nÞE  n 2 4

and

F ct

1 1 ¼ 2r PðnÞ

R¼

U c ¼ gE  n 

g 2

þ

3 h 4

i 1 5n  E  n  trðEÞ n þ aE  n 5

Now, if we compare with Eq. (38) for the isotropic case in which results must coincide for a second or a fourth order tensor. We deduce:

a¼

25 2

The kinematic assumption obtained from the fourth order fabric tensor is then:

uc ¼



 i 2r g 3 h 5n  E  n  trðEÞ n gE  n  þ PðnÞ 2 4 25 D<ðn  n  n  nÞE  n þ 2

ð42Þ

The normal and tangential components of the displacement become:

ucn ¼



 2r 3g 15 g 3 trðEÞ E : ðn  nÞ þ  þ þ PðnÞ 4 2 2 4 25 D<ðn  n  n  n  n  nÞ : E þ 2

and

uct ¼

2r PðnÞ



gþ

 25 D<ðn  n  n  nÞ T : E 2

The normal and tangential contact forces are obtained from the local contact law as:



  2r 3g 15 g 3 trðEÞ F cn ¼ E : ðn  nÞ þ Kn  þ þ PðnÞ 4 2 2 4  25 D<ðn  n  n  n  n  nÞ : E þ 2

ð43Þ

and then:

Chom ¼

1 ðd  d þ 2ddÞ 15

 25 D<ðn  n  n  nÞ K t T : E 2

 

Z 4r 2 N 3g 15 1 n  n  n  nds Kn  þ V 4 4p s2 2

Z g 3 1 þ n  nds  d þ 2 4 4 p s2  Z 25 1 sn  n  n  n  n  n  n  nds D< þ 2 2p s  Z  1 25 g þ D<ðn  n  n  nÞ T T  Tds : E þ Kt 2p s 2

In order to identify the constant a, we are going to consider the isotropic case, for which:

Eq. (30) becomes:

gþ

In the same way, from (12) we get:

ð41Þ

D¼



  Z 4r 2 N 3g 15 1 Kn  þ n  n  n  nds V 4 4p s 2 2

Z g 3 1 n  nds  d þ þ 2 4 4p s 2  Z 25 1 n  n  n  n  n  n  n  nds þ D< 2 2p s  sym Z  1 25 þ Kt g þ D<ðn  n  n  nÞ T T  Tds 4p s2 2

By referring to Appendix C, in which can be found the analytical expressions of the integrals (cf. (C.2), (C.4), (C.5) and (C.6)), we get: 



4r2 N g 1 g 1 Chom ¼ dd Kn  ðd  d þ 2ddÞ þ þ þ V 10 4 6 4 5 þ ðd  d þ 2dd þ 8D þ 4D  d þ 4d  D þ 8Dd þ 8dDÞ 126  4r 2 N g g 5  ðd  d þ 2ddÞ þ dd  ðd  d þ 2dd þ Kt V 126 15 3 þ 8D þ 4D  d þ 4d  D þ 8Dd þ 8dDÞ 5 ðdd þ 2Dd þ 2dDÞ þ 14

ð44Þ

4.4. The particular case of an isotropic distribution of contacts In the isotropic case, it is convenient to compare the results given by the kinematical rule with to the one provided the Voigt assumption as well as to existing numerical results reported by Ref. [2]. Indeed in isotropic case, expression (40) reads4:

Chom ¼

  



4Nr2 g 1 g 1 dd dd Kn þ þ V 15 6 5 6  



g 1 g 1 þ Kt  ddþ dd þ þ 15 6 5 2

from which it follows that:

8 < lhom ¼ E ¼ 1 2ð1þmÞ 15

Nr2 K n V

: khom ¼ E ¼ 4 3ð12mÞ 9

Nr2 K n V

ðð5 þ 6gÞ þ 3að5 þ 2gÞÞ

ð45Þ

4 It is recalled that in the isotropic case, the use of a second or fourth order fabric tensor leads to the same result.

875

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

5. Analysis of the results

.

5.1. The isotropic case

.

The results obtained from the two localization rules (Voigt and kinematical) are represented in Table 1 in the isotropic case. We note, from these results (see Table 1), that the macroscopic hom does not depend on the considered localization bulk modulus k rule and is a linear function of the normal rigidity in the contact K n . We also observe (cf. Fig. 3) that the variation of lhom obtained from the kinematical rule is linear with respect to a. This is in agreement with the numerical results of [2], which corresponds to the discrete points added in Fig. 3. For g ¼ 1; 2, the variations of lhom (linear in a) seems to be in agreement with the numerical results. The same result is found by considering the variations of hom and the macroscopic Young modulus Ehom (obtained from k hom l by standard relationship) with respect to a (see Fig. 4). We can then conclude that the general kinematic assumption proposed in the present paper allows to approach accurately the numerical simulations results.

. . . . . . .

.

.

.

.

.

5.2. Illustration of the results for an anisotropic distribution of contacts For this purpose of the analysis, it now convenient to introduce a spherical coordinate system. Then an arbitrary unit vector n can be represented in this coordinate system as:

n ¼ ½sinðhÞ cosð/Þ; sinðhÞ sinð/Þ; cosðhÞT ;

h½0; p;

/½0; 2p

Consider a random packing of spherical particles with r ¼ 0:7 mm radius, kn ¼ 2 and kt ¼ 1 rigidities and flexibilities and N ¼ 1956

2

k

lhom

Proposed kinematical rule 4 Nr 2 K n 9 V 1 Nr 2 K n ðð5 þ 6gÞ þ 3að5 þ 2gÞÞ 15 V

4 Nr K n 9 V 2 Nr 2 K n ð2 þ 3aÞ 15 V

hom

total number of contacts. The packing volume is V ¼ 4000 mm3 . The value g ¼ 1:2 is chosen for the adjustable parameter. The Young modulus related to the direction of unit vector m reads:

1 EðmÞ ¼ ðm  mÞ : Shom ðm  mÞ

Table 1 Comparaison of the considered localization rules. Voigt assumption

3V Ehom with respect to a by using the assumption of Voigt, 2Nr 2 kn kinematic 1 (g ¼ 1; 4), kinematic 2 (g ¼ 1; 35) and kinematic 3 (g ¼ 1; 2Þ.

Fig. 4. Variation of

ð46Þ

5.2.1. Trigonometrical orientation distribution of contacts Consider the following trigonometrical distribution of grains’s contact orientation:

PðnÞ ¼

3 2 sin ðarccosðn  e3 ÞÞ 2

ð47Þ

T

where e3 ¼ ½0; 0; 1 : The representation of this distributions is given in Fig. 5. The second and fourth order fabric tensor determined for this contact distributions, are given by:

. .

2

.

6 ½D ¼ 4

.

2 .

0:4000

0

0

0

0:4000

0

3 7 5 and

0 0 0:2000 0:2571 0:0857 0:0571

0

0

0

3

6 0:0857 0:2571 0:0571 0 0 0 7 7 6 7 6 6 0:0571 0:0571 0:0857 0 0 0 7 7 ½D ¼ 6 6 0 0 0 0:0571 0 0 7 7 6 7 6 4 0 0 0 0 0:0571 0 5 0 0 0 0 0 0:0857

. . . .

.

.

.

.

.

.

. .

. . 3V lhom with respect to a by using the assumption of Voigt, 2Nr 2 K n kinematic 1 (g ¼ 1; 4), kinematic 2 (g ¼ 1; 35) and kinematic 3 (g ¼ 1; 2Þ.

Fig. 3. Variation of

.

Fig. 5. Tensorial representation of the trigonometrical distribution of contact orientation.

876

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

.

.

.

.

.

.

. .

.

.

.

.

. .

Fig. 6. Young modulus giving by the Voigt and kinematical assumption and by means of second order fabric tensor of general distribution (in the left the Voigt assumption and in the right the kinematical assumption).

It must be emphasized that in this case of trigonometric distribution, the second order fabric tensor as well as the fourth order tensor produce the exact description of the contact distribution (shown in Fig. 5). The homogenized elastic stiffness tensors giving by the Voigt assumption (Cv ) and kinematical assumption (Ck ) in the case of a second order fabric tensor take the form:

2

0:6298 0:0821 0:0547

0

0

6 0:0821 0:6298 0:0547 0 0 6 6 6 0:0547 0:0547 0:2738 0 0 6 ½Cv  ¼ 6 6 0 0 0 0:1985 0 6 6 0 0 0 0:1985 4 0 0

0

0

0 0 0 0 0

0

0

0:2738

0:5914 0:1013 0:0739

0

0

0

6 0:1013 0:5914 0:0739 0 0 6 6 6 0:0739 0:0739 0:2355 0 0 6 ½Ck  ¼ 6 6 0 0 0 0:1697 0 6 6 4 0 0 0 0 0:1697 0

0

0

0

0

0 0 0 0

Pðh; /Þ ¼

63 189 sinðhÞ cosð/Þ cosðhÞ þ sinðhÞcos2 ð/Þ sinð/Þ cosðhÞ 2 4 þ

3; 224; 991 37; 799; 937 cos2 ð/Þ þ cos4 ð/Þ cos2 ðhÞ 1; 000; 000 4; 000; 000

3 7 7 7 7 7 7 7 7 7 5



37; 799; 937 cos4 ð/Þ cos4 ðhÞ 8; 000; 000

þ

656; 187 cos2 ðhÞ cos2 ð/Þ 2; 500; 000



17; 437; 329 37; 799; 937 cos4 ðhÞ cos2 ð/Þ  cos4 ð/Þ 5; 000; 000 8; 000; 000

3



189 sinðhÞ cos2 ð/Þ sinð/Þ cos3 ðhÞ 4

þ

4; 218; 757; 003 3; 600; 009 cos4 ðhÞ  cos2 ðhÞ 500; 000; 000 800; 000



63 63 sinðhÞ cosðhÞ cos3 ð/Þ  cos3 ðhÞ sinðhÞ cosð/Þ 2 2

þ

63 3374991 sinðhÞ cos3 ðhÞ cos3 ð/Þ þ 2 4; 000; 000

and 2

5.2.2. General orientation distribution of contacts Consider the following contact distribution:

7 7 7 7 7 7 7 7 7 5

0:2450

Although the are based on the same representation of contact distribution, these predictions present some differences on various modulus. The relative difference attains 40% Young modulus (Fig. 6).

which corresponds to a distribution exactly described by the fourth order fabric tensor whose components are:

. . . . .

ð48Þ

.

.

. .

. .

. . .

Fig. 7. Second and fourth tensorial representation of the general distribution of grain orientations.

.

877

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

.

. . .

.

.

.

.

Fig. 8. Young modulus giving by the Voigt assumption and by means of second and fourth order fabric of general distribution.

3 0:0357 0:0671 0:0471 0 0 0 6 0:0671 0:2071 0:0857 0 0 0 7 7 6 7 6 6 0:0471 0:0857 0:4071 0 0 0 7 7 6 ½D ¼ 6 7 7 6 0 0 0 0:0571 0 0:1000 7 6 7 6 4 0 0 0 0 0:0857 0:1500 5 2

0

0

0

0:1000

With the fourth order fabric tensor we get:

2

0:1500 0:0571

The approximation of this contact distribution by means of a second order fabric tensor is given by:

2 6 ½D ¼ 4

3

0:1500

0

0

0

0:3500

0

0

0

0:5000

Representations of the exact (fourth order) contact distribution and the approximate distribution (49) are shown in Fig. 7. The homogenized elastic stiffness tensor given by the Voigt assumption (Cc ) and by the kinematical assumption (Ck ) combined with the second order fabric tensor takes the form:

0:1848

0:0410 0:0616

0

0

0:1779

3

7 6 6 0:0410 0:5408 0:0889 0 0 0:1779 7 7 6 7 6 7 6 0:0616 0:0889 0:8078 0 0 0:0273 7 6 hom ½Cv  ¼ 6 7 7 6 0 0 0 0:2926 0 0 7 6 7 6 6 0 0 0 0 0:2173 0 7 5 4 0:1779 0:1779 0:0273

0

0

0:1608

0:1465

0:0807

0

0

0

0:5024 0:1081

0

0

0

0:1081 0:7694

0

0

0

and 2

6 6 0:0602 6 6 6 0:0807 6 ½Chom  ¼ 6 k 6 0 6 6 6 0 4

0:0602

0

0

0:2639

0

0

0

0

0

0:1886

0

0

0

0:1321

0:1779 0:1779 0:0273

0:0451

0

0

3

0

7 6 6 0:0643 0:5339 0:0821 0 0 0 7 7 6 7 6 7 6 0:0451 0:08215 0:8694 0 0 0 7 6 ½Chom  ¼ 7 6 v 7 6 0 0 0 0:2584 0 0:0958 7 6 7 6 7 6 0 0 0 0 0:2378 0:1437 5 4 0 2

9 21 15 3 sinðhÞ2 cosð/Þ2 þ sinðhÞ2 sinð/Þ2 þ cosðhÞ2  8 8 4 2 ð49Þ

2

0:0643

0

0

0:0958 0:1437 0:1745

and

7 5

which provides:

Pðh; /Þ ¼

0:1779

3 7 7 7 7 7 7 7 7 7 7 7 5

0:0044 0:0578 0:0745

0

0

0

3

7 6 6 0:0578 0:5130 0:1167 0 0 0 7 7 6 7 6 7 6 0:0745 0:1167 0:9161 0 0 0 7 6 ½Chom  ¼ 7 6 k 7 6 0 0 0 0:2826 0 0:0304 7 6 7 6 7 6 0 0 0 0 0:1924 0:0456 5 4 0

0

0

0:0304

0:0456 0:1095

These results shows strong difference between the two localisations assumptions. The values of the derived Young moduli (Fig. 8) confirm this observation. 6. Conclusion In the present paper, a multiscale approach of the anisotropic elastic properties of the granular materials is presented. This approach takes advantage of the representation of the grains contact distribution by means of second order and fourth order fabric tensor. The consideration of these fabric tensors with various homogenization schemes allows to quantify the overall anisotropy of the material. In particular, we obtain new results based on a proposed kinematical localization rule. In the case of an isotropic material, the predictions of elastic moduli are in a very good agreement with the available numerical results of Chang and Liao [2]. This agreement provides a first validation of the developed approach. For a more advanced validation, we compare the results obtained here by applying the models to some theoretical contact distribution. The last illustration allows to demonstrate the interest of the fourth order fabric tensor. It will be particulary interesting to proceed to a more complete validation, for instance by performing numerical discrete computations of elastic moduli [39].

878

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

In the same way, we compute the integral:

Appendix A. Tensorial identities Let us first recall the following identities (see for instance [23,40]) which will required in various computations:

1 4p 1 4p 1 4p 1 4p 1 4p

Z

1 4p

2

Z Z

s2

ni nj ds ¼

1 dij 3

1 ni nj nk nl ds ¼ J ijkl 5 s2

PðnÞn  d  nds i R h15 ¼ D : ðn  n  32Þ n  d  nds s2 2 R ¼ 15 D : 41p s2 n  n  n  d  nds 2 R  32 41p s2 n  d  nds ¼ Dd

ðB:4Þ

The following identities are also obtained by using the results (A.1):

ðA:1Þ

1 J 7 ijklmn s2 Z 1 ni nj nk nl nm nn np nq ds ¼ J ijklmnpq 9 s2 ni nj nk nl nm nn ds ¼

Z

1 4p

n  n  n  nds ¼

s2

1 ðd  d þ 2ddÞ 15

ðB:5Þ

and

Z

1 4p

with

n  nds  d ¼

s2

1 dd 3

ðB:6Þ

In the same way, by using (A.1), we get:

1 ðdij dkl þ dik djl þ dil djk Þ 3 1 J ijklmn ¼ ðdij J klmn þ dik J jlmn þ dil J jkmn þ dim J jkln þ din J jklm Þ 5 1 J ijklmnpq ¼ ðdij J klmnpq þ dik J jlmnpq þ dil J jkmnpq þ dim J jklnpq 7 þ din J jklmpq þ dip J jklmnq þ diq J jklmnp Þ J ijkl ¼

1 4p

D:

Z

n  n  n  n  n  nds

s2

1 2 ðd  d þ 2ddÞ þ ðD  d þ d  D þ 2Dd þ 2dDÞ 105 105 ðB:7Þ

¼

and

Appendix B. Computation of some required integrals in the case of the second order fabric tensor

Z

1 4p

T T  Tds Z  h i 1 1 ðn  nÞd þ dðn  nÞ  n  n  n  n ds ¼ 4 p s2 2 1 1 ¼ dd  ðd  d þ 2ddÞ 3 15

Following Lubarda and Krajcinovic [2], we use the expression (5) of PðnÞ and the results (A.1) (see [40,23]) to compute the following integral:

Z 1 PðnÞðn  nÞds 4p s2   Z 1 15 1 ¼ D : ðn  nÞ  ðn  nÞds 4p s2 2 5   Z Z 1 15 3 ¼ D: n  n  n  nds  ðn  nÞds 4p 2 2 s2 s2

s2

1 4p

ds ¼ 1

Zs

R

s2

ðB:8Þ

Finally, we have:

D¼

D: ðB:1Þ

In the same way, from (5) and (A.1), the following identity can established: Z 1 PðnÞðn  n  n  nÞds 4p s2   Z 1 15 1 D : ðn  nÞ  ðn  n  n  nÞds ¼ 4p s 2 2 5   Z Z 1 15 3 D: ðn  n  n  n  n  nÞds  ðn  n  n  nÞds ¼ 4p 2 2 s2 s2   1 1 2  ðd  dÞ  ðddÞ þ ðD  d þ d  DÞ þ 2ðDd þ dDÞ ¼ 7 5 5

1 4p

Z

ðn  nÞ  T T  Tds Z  h i 1 1 ¼D: ðn  n  n  nÞd þ n  n  dðn  nÞ 4p s2 2 s2

 n  n  n  n  n  ngds 1 ¼ ðdd þ Dd þ dDÞ 15  1 2 ðd  d þ 2ddÞ þ ðD  d þ d  D þ 2Dd þ 2dDÞ  105 105 ðB:9Þ where the use has been done of the following definition:

ðABÞijklmn ¼

1 ðAijkn Blm þ Aijknm Bln Þ 2

ðB:2Þ Knowing that T ¼ n  I  n  n  n and T T ¼ I  n  n  n  n, we have:

1 4p

Z

PðnÞðT T :TÞds Z h i 15 1 ¼ ðn  nÞ  dðn  nÞ þ ðn  nÞd ds D: 2 4p s 2 Z h i 3 1  dðn  nÞ þ ðn  nÞd ds 4 4 p s2 Z 1  PðnÞ  ðn  n  n  nÞds 4p s2   1 1 2 3 d  d þ ðddÞ  ðD  d þ d  DÞ þ ðDd þ dDÞ ðB:3Þ ¼ 7 5 5 2 s2

Appendix C. Computation of some required integrals in the case of the fourth order fabric tensor We are going to use the expression (6) of PðnÞ obtained when considering the fourth order fabric tensor, and the expressions (A.1). We get:

Z 1 PðnÞðn  nÞds 4p s 2  Z 15 1 ¼ ðn  n  n  n  n  nÞds 21D< 8 4p s 2  Z Z 1 1 ðn  n  n  nÞds þ ðn  nÞds  14D : 4p s2 4p s 2

D¼

ðC:1Þ

879

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

We have then:

1 D¼ 4p ¼

Z 

C 11

PðnÞðn  n  n  nÞds s2

C 13

Z

15 1 ðn  n  n  n  n  n  n  nÞds 21D< 8 4p s2 Z 1  14D : ðn  n  n  n  n  nÞds 4p s 2  Z 1 ðn  n  n  nÞds þ 4 p s2

C 15 C 22

ðC:2Þ

C 26

and

1 4p

C 24

Z

C 34

PðnÞðT T :TÞds

s2

¼

¼

     3 3 1 ¼ 6D11  K n þ D11 þ K t C 12 ¼ D11 þ D22  ðK n  K t Þ 5 5 5   1 ¼ D11 þ D33  ðK n  K t Þ C 14 ¼ D23 ðK n  K t Þ 5



1 1 ¼ D13 3K n þ K t C 16 ¼ D12 3K n þ K t 2 2       3 3 1 ¼ 6D22  K n þ D22 þ K t C 23 ¼ D22 þ D33  ðK n  K t Þ 5 5 5

1 ¼ D23 3K n þ K t C 25 ¼ D31 ðK n  K t Þ 2    

1 3 3 ¼ D12 3K n þ K t C 33 ¼ 6D33  K n þ D33 þ K t 2 5 5



1 1 ¼ D23 3K n þ K t C 35 ¼ D31 3K n þ K t 2 2     1 3 3 1 ¼ D22 þ D33  K n þ D33 þ D22 þ K t C 36 ¼ D12 ðK n  K t Þ 5 4 4 5



3 3 ¼ D12 K n þ K t C 46 ¼ D13 K n þ K t 4 4    

1 3 3 1 3 ¼ D11 þ D33  K n þ D11 þ D33 þ K t C 56 ¼ D32 K n þ K t 5 4 4 5 4     1 3 3 1 ¼ D11 þ D22  K n þ D22 þ D11 þ K t 5 4 4 5

 Z h i 15 21 1 ðn  n  n  nÞ  ðn  nÞd þ dðn  nÞ ds D< 8 2 4p s2 Z h i 1  7D : ðn  nÞ  ðn  nÞd þ dðn  nÞ ds 4p s2 Z h i 1 þ n  nÞd þ dðn  nÞ ds 8p s2  Z 1  PðnÞðn  n  n  nÞds 4p s2

C 44

1 ðDd þ dDÞ  D 2

Appendix E. Elements of the kinematical localization rule

C 45 C 55 C 66

ðC:3Þ We first have:

As well as

1 4p



Z

PðnÞn  d  nds ¼ s2

1 4p

Eij ¼

Z

15 ½21D<ðn  n  n  nÞ 8 i  14D : ðn  nÞ þ 1 n  d  ndsd s2

ðC:4Þ

Now, by using the results (A.1), we obtain:

Z 25 1 n  n  n  n  n  n  n  nds D< 2 4p s2  5

d  d þ 2dd þ 8D þ 4D  d þ 4d  D þ 8Dd þ 8dD ¼ 126 ðC:5Þ and finally: Z 25 1 ðn  n  n  nÞT T  Tds D< 2 4p s2 Z  h 25 1 1 ¼ D< ðn  n  n  n  n  nÞd 2 4p s2 2 i þ n  n  n  n  dðn  nÞ

  Z Z 2Nr 1 1 Akj A0 ni nk ds þ A1 trðEÞ ni nk ds V 4p s 2 4p s 2 Z Z 1 1 ni nk ds þ A3 Epq np nq ni nk ds þ A2 trðE  DÞ 4p s2 4 p s2 Z Z 1 1 þ A4 ðE  DÞpq np nq ni nk ds þ A5 trðEÞDpq np nq ni nk ds 4p s2 4p s 2 Z 1 þ A6 Epq np nq ni nk nm nn dsDnm 4p s 2 Z Z 1 1 þ A7 Dpq np nq ni nk ds þ B1 Ein nn nq np nk dsDpq 4p s 2 4p s2 Z Z 1 1 þ B2 Eip np nk ds þ C 1 Din nn nq np nk dsEpq 4p s2 4p s2 Z Z 1 1 þ C 2 trðEÞDip np nk ds þ C 3 Dip np nk ds 4p s 2 4p s 2  sym Z np nk ds ðE:1Þ þ D1 ðE  D  D  EÞip s2

 n  n  n  n  n  n  n  ngds  5

¼ dd þ 2Dd þ 2dD 14  5

d  d þ 2dd þ 8D þ 4D  d þ 4d  D þ 8Dd þ 8dD  126 ðC:6Þ

Appendix D. Components of the homogenized stiffness, in Voigt notation By adopting for Chom the classic notation of Voigt which allows to represent Chom as a symmetrical matrix 6  6, we get from (24):

V At this step, it is convenient to take A ¼ 2Nr D1 . Computing of the integrals of (E.1) (cf. Appendix B2 ), we get:  i 2Nr 1 h Eij ¼ A0 Aij þ A1 trðEÞAij þ A2 trðE  DÞAij þ A3 trðEÞAij þ 2Eik Akj 3V 5   1 V þ A4 trðE  DÞAij þ Eij þ Ekl Dli Akj 5 2Nr   1 V þ A5 trðEÞAij þ trðEÞdij 5 Nr  1 V þ A6 trðEÞAij þ trðEÞdij þ 2trðE  DÞAij þ 4Elk Akj Dil 35 Nr      V 1 V 1 V þ 2 Eij þ 2Eik Akj þ A7 Aij þ dij þ B1 Eij þ Eik Akj Nr 5 Nr 5 Nr   1 V V þ B2 Eik Akj þ C 1 2Ekl Dli Akj þ trðEÞdij þ C 2 trðEÞdij 5 2Nr 2Nr   sym V V þ C3 ðE:2Þ dij þ D1 Eij  Elk Akj Dil 2Nr 2Nr

880

J. Rahmoun et al. / Computational Materials Science 46 (2009) 869–880

Since this last relation must be satisfied for any value of D and E, we get by identification the following system of eight equations with 14 unknowns A0 to A7 ; B1 ; B2 ; C 1 ; C 2 ; C 3 and D1 :

8 A0 þ 15 A7 ¼ 0; 15 A7 þ 12 C 3 ¼ 0 > > > > 1 > A6 ¼ 0 A1 þ 15 A3 þ 15 A5 þ 35 > > > > 1 2 > þ A þ A ¼ 0 A > < 2 5 4 35 6 1 1 1 A þ 35 A6 þ 10 B1 þ 12 B2 ¼ 0 5 3 > > 1 1 1 > A þ 35 A6 þ 10 C 1 þ 12 C 2 ¼ 0 > > 5 5 > > 1 2 > > 10 A4 þ 35 A6 þ 15 C 1  12 D1 ¼ 0 > > :1 4 2 A þ 105 A6 þ 15 B1 þ 13 D1 ¼ 1 15 4 To solve this system, let us fix 6 of the 14 unknowns, for example A6 ¼ k1 ; A5 ¼ k2 ; D1 ¼ k3 ; C 2 ¼ k4 ; A0 ¼ b; B2 ¼ g and let us express the other unknowns with respect to ki ; b and g. We get:

8 A0 > > > > > A1 > > > > > A2 > > > > > > < A3

¼ b;

B2 ¼ g

1 ¼ 34 þ g2  35 k1  25 k2  12 k3  12 k4 2  ¼  35 k1 þ 45 k2 þ k3 þ 2k4

¼  15  52 g þ k2 þ 52 k3 þ 52 k4 4

A4 ¼ 4k2 þ 5k3 þ 10k4 > > > A5 ¼ k2 ; A6 ¼ k1 ; A7 ¼ 5b > > > > 15 2 > > B 1 ¼ 2  7 k1  2k2  5k3  5k4 > > > > 2 > C 1 ¼  7 k1  2k2  5k4 > > : C 2 ¼ k4 ; C 3 ¼ 2b; D1 ¼ k3 k1 ; k2 ; k3 ; k4 ; b and g are constants, considered as the characteristic of the granular medium. Then, by introducing the values of Ai ; Bi ; C i and Di in (30) of U c , we get the general kinematic assumption given by (34). References [1] B. Cambou, Ph. Dubujet, F. Emeriault, F. Sidoroff, Eur. J. Mech. A/Solids 14 (1995) 225–276. [2] C.S. Chang, C.L. Liao, Appl. Mech. Rev. 47 (1 Part 2) (1994) 197–207. [3] M. Oda, S. Nemat-Nasser, J. Konisjh, Soils Found. 25 (1985) 85–97. [4] J. Lanier, M. Jean, Powder Technol. 109 (1) (2000) 206221. [5] R.J. Bathurst, L. Rothenburg, Mech. Mater. 9 (1990) 65–80. [6] C.S. Chang, Micromechanical Modelling of Constitutive Equation for Granular Materials. Micromechanics of Granular Materials, Elsevier Science Publishers, 1988. pp. 271–278. [7] K. Kanatani, Int. J. Eng. Sci. 22 (1984) 149–164. [8] D. Krajcinovic, Damage Mechanics, North-Holland, The Netherlands, 1996. [9] M. Madadi, O. Tsoungui, M. Lätzel, S. Luding, Int. J. Solids Struct. 41 (2004) 2563–2580.

[10] F. Sidoroff, B. Cambou, A. Mahboubi, Mech. Mater. 16 (1993) 83–89. [11] L. Louis, C. David, V. Metz, P. Robion, B. Menendz, C. Kissel, Int. J. Rock Mech. Min. Sci. 42 (2005) 911923. [12] L. Louis, P. Robion, C. David, D. Frizon de Lamotte, J. Struct. Geol. 28 (2006) 549–560. [13] C. Nougier, Simulation des Interactions Outil-sol: Application aux Outils de Traitement des Sols, Ph.D. Thesis, Département de Modélisation Physique et Numérique de l’Institut de Physique du Globe de Paris, 1999. [14] C. Nouguier-Lehon, P. Dubujet, B. Cambou, Analysis of granular material behaviour from two kinds of numerical modelling, in: 15th ASCE Engineering Mechanics Conference, Columbia University, New York, NY, June 2002. [15] B. Cambou, M. Chaze, F. Dedecker, Eur. J. Mech. A/Solids 19 (2000) 999–1014. [16] B. Cambou, Ph. Dubujet, C. Nouguier-lehon, Mech. Mater. (2004) 1–10. [17] C.S. Chang, S.J. Chao, Y. Chang, Int. J. Solids Struct. 32 (14) (1995) 1989– 2008. [18] E. Masad, D. Little, R. Lytton, Int. J. Geomech. 4 (2004) 254–263. [19] J. Rahmoun, Modélisation du Comportement des Matériaux Granulaires par des Approches Discrètes et Continues, Ph.D. Thesis, University of Lille, 2006. [20] C.S. Chang, J. Gao, Acta Mech. 115 (1996) 213–229. [21] O. Millet, S. Gu, D. Kondo, Comp. Rend. Mécanique 335 (2007) 231–237. [22] M. Satake, Constitution of mechanics of granular materials through graph theory, in: S.C. Cowin, M. Satake (Eds.), Proceedings of the US–Japan Seminar Continuum Mechanical and Statistical Approaches in the Mechanics of Granular Materials, Elsevier, Amsterdam, The Netherlands, 1978, pp. 47–62. [23] V.A. Lubarda, D. Krajcinovic, Int. J. Solids Struct. 30 (20) (1993) 2859–2877. [24] V. Pensée, Contribution de la Micromécanique à la Modélisation Tridimensionnelle de l’endommagement par Mésofissuration, Ph.D. Thesis, University of Lille, 2002. [25] J. Christoffersen, M.M. Mehrabadi, S. Nemat-Nasser, J. Appl. Mech., ASME 48 (2) (1981) 339–344. [26] A.E.H. Love, A Treatise of Mathematical Theory of Elasticity, University Press, Cambridge, 1927. [27] L. Rothenburg, A.P.S. Selvadurai, Micromechanical Definition of the Cauchy Stress Tensor for Particulate Media, Mechanics of Structured Media, Elsevier, Amsterdam, The Netherlands, 1981. pp. 469–486. [28] J. Weber, Recherches Concernant les Contraintes Intergranulaires dans les Milieux Pulvirulents. bul. Liaison P. and Ch., Number Juil.-août, 1966. [29] B. Cambou, M. Jean, Micromécanique des Matériaux Granulaires, Hermes Science, 2001. [30] R.B. Chapuis, De la Structure des Milieux Granulaires en Relation avec leur Comportement Mécanique, Ph.D. Thesis, D.E., Montréal, Canada, 1976. [31] K. Walton, J. Mech. Phys. Solids 35 (3) (1987) 213–226. [32] R.J. Bathurst, L. Rothenberg, J. Appl. Mech., ASME 55 (1988) 17–23. [33] F. Emeriault, B. Cambou, A. Mahboubi, Mech. Cohes.-Frict. Mater. 1 (1996) 199–218. [34] F. Emeriault, C.S. Chang, Comput. Geotech. 20 (3/4) (1997) 223–244. [35] C.S. Chang, C.L. Liao, Int. J. Solids Struct. 26 (4) (1990) 437–453. [36] J.T. Jenkins, Inelastic behavior of random arrays of identical spheres, in: N.A. Fleck, A.C.E. Cocks (Eds.), IUTAM Symposium on Mechanics of Granular and Porous Materials, Kluwer Academic Publishers, 1997, pp. 11–22. [37] J.P. Boehler, Application of Tensors Functions in Solids Mechanics, CISM Courses and Lectures, vol. 292, Springer-Verlag, Wien, New York, 1987. [38] A.J.M. Spencer, Isotropic polynomial invariants and tensor functions, in: J.P. Beohler (Ed.), Applications of Tensor Functions in Solid Mechanics, CISM Courses and Lectures, vol. 292, Springer, Berlin, 1987, pp. 141–169. [39] J. Fortin, O. Millet, G. de Saxcé, Int. J. Numer. Meth. Eng. (62) (2005) 639–663. [40] Q.-C. He, A. Curnier, Int. J. Solids Struct. 32 (10) (1995) 1433–1457.