Enhanced micropolar model for wave propagation in ordered granular materials

(1)

Contents lists available at ScienceDirect

International

Journal

of

Solids

and

Structures

journal homepage: www.elsevier.com/locate/ijsolstr

Enhanced

micropolar

model

for

wave

propagation

in

ordered

granular

materials

A. Merkel

a , ∗

_,

_S.

_Luding

a

a Multi Scale Mechanics, TS, CTW, UTwente, P.O. Box 217, 7500 AE Enschede, Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 2 December 2015 Revised 31 August 2016

Available online 29 November 2016 Keywords:

Wave propagation Granular Micropolar

a

b

s

t

r

a

c

t

Thevibrationalproperties ofaface-centeredcubic granularcrystalofmonodisperse particlesare pre-dictedusingadiscretemodelaswellastwomicropolarmodels,firsttheclassicalCosseratandsecond anenhancedCosserat-typemodel,thatproperlytakesintoaccountalldegreesoffreedomatthecontacts betweentheparticles.The continuummodelsarederived fromthediscretemodel viaamicro–macro transitionofthediscreterelativedisplacementsandparticlerotationstotherespectivecontinuumfield variables.Next,onlythelongwavelengthapproximationsofthemodelsarecomparedand,considering the discretemodel as reference, theCosserat model showsinconsistent predictions ofthe bulkwave dispersionrelations.Thiscanbeexplainedby aninsufficientmodelingofsliding modeofparticle in-teractionsinthe Cosseratmodel.Anenhancedmicropolar model isproposedincludingonlyone new elastictensorfromthemorecompletesecond–ordergradientmicropolartheory.Thisenhanced micropo-larmodeltheninvolvestheminimumnumberofelasticconstantstoconsistentlypredictthedispersion relationsinthelongwavelengthlimit.

1. Introduction

The classical theory of elasticity consists of a macroscopic material description. The material is not described at the micro-level by considering the displacement of the different particles in interaction, but is described as a continuum by considering macroscopic quantities as stress and strain. The classical elasticity theory can be viewed as a ﬁrst gradient of the displacement ﬁeld approximation of solid state theory ( Ashcroft and Mermin, 1976 ) and is thus valid in the long wavelength limit only. Granular media, due to their micro-inhomogeneous character, are not well described by the standard continuum theory of elasticity. In contrast to classical continua, where the sizes of the vibrating atoms can be assumed to be negligible compared to the macroscale, the sizes of the particles in a granular assembly are comparable to it ( Schwartz et al., 1984 ). In addition, considering the sliding, twisting and rolling resistances at the level of the contacts between the particles, a consistent description of the elasticity of a granular medium needs to take into account all the rotational degrees of freedom of each in- dividual particle and thus all the relative degrees of each pair.

By including the rotational degrees of freedom into the analysis, the evaluation of the different elastic constants in the quasistatic

∗ _{Corresponding author at: Multi Scale Mechanics, TS, CTW, UTwente.}

E-mail addresses: aurelien.merkel.etu@univ-lemans.fr (A. Merkel), s.luding@utwente.nl (S. Luding).

behavior of granular assemblies becomes complex ( Jenkins, 1990; Jenkins and Ragione, 20 01; 20 03; Ragione and Jenkins, 2009; Ra- gione and Magnanimo, 2012 ). The elastic behaviors of crystalline structures of monodisperse beads can be eﬃciently described by a discrete model, where the displacement and rotation of each indi- vidual bead are taken into account. One of the major differences with classical elasticity is the existence of optical-type rotational- related modes of wave propagation. Especially, the dispersion relations given by the discrete model of the bulk waves propagating in crystalline structures of contacting monodisperse beads, without solid bridges between them, is consistent with experimental results ( Merkel et al., 2010; 2011 ) in a hexagonal close-packed structure, and with numeric simulation results in a face-centered cubic structure ( Mouraille et al., 20 06; 20 09; Mouraille, 20 08 ). Neverthe- less, the discrete model can be solved analytically only for well- known regular crystalline structures, the case of a random assembly of beads can be done ( Kruyt, 2010; 2012 ) but is too complex for large systems. Moreover, due to diffraction scattering and attenua- tion effects, even for a small level of randomness, only long wavelength waves will propagate in a granular assembly ( Mouraille and Luding, 2008; Dazel and Tournat, 2010 ). More speciﬁcally, only the long wavelength waves can propagate as a coherent ballistic wave ( Jia et al., 1999; Jia, 2004; Langlois and Jia, 2014 ). Considering that only the long wavelength waves will propagate in random assemblies, which differs from the ideal crystalline case, and that their discrete modeling is too complex, a continuum formulation seems http://dx.doi.org/10.1016/j.ijsolstr.2016.11.029

(2)

more suitable. A continuum model is also relevant when considering the wave propagation in a granular assembly in contact with an elastic solid ( Wallen et al., 2015 ).

The generalization of the classical elasticity theory accounting for the rotational degrees of freedom of bodies is known as the Cosserat theory ( Cosserat and Cosserat, 1909; Eringen, 1999 ). For instance, a Cosserat model can explain the strain localization in a sheared fault gauge ( Pasternak et al., 20 03; 20 04; 20 06 ). But, from a direct comparison between the dispersion relations predicted with a discrete model and with the ones of a Cosserat model, the latter does not predict correctly the dispersion relations of the modes of propagation related to the rotational motion of the beads ( Merkel et al., 2011 ). Similarly, the simulation results of a Couette shear cell differ from the predictions of a Cosserat model ( Mohan et al., 1999; Lätzel et al., 2001; Mohan et al., 2002; Latzël, 2003 ). Despite many theoretical effort s, the comparison between experimental results and predictions from the Cosserat theory remains inconclusive; a consistent continuum description of granular assemblies is still challenging, see Maugin and Metrikine (2010) and Goddard (2014) and the references therein. As it will be shown below, one of the interactions, or relative degrees of freedom, between the beads involving the rotational degrees of freedom is not modeled properly leading to inconsistent results. A reconstruc- tion of the Cosserat moduli ( Pasternak and Dyskin, 2014 ) can lead to misleading results. The drawback of the Cosserat theory can be overcome using a second-order gradient theory ( Mühlhaus and Oka, 1996; Suiker et al., 20 01a; 20 01b; Suiker and de Borst, 2005 ), where the homogenization is performed by differential expansions ( Pasternak and Mühlhaus, 2002b; Pasternak and Mühlhaus., 2005 ). Nevertheless, the second-order gradient micropolar theory intro- duces three new elastic tensors, which involve too many new elastic constants to represent a feasible alternative to the discrete modeling. A homogenization by integral transformation resulting into a non-local Cosserat continuum theory gives the same results as the discrete model, but this approach does not provide any simpliﬁcation compared to the discrete model, see Pasternak and Mühlhaus (2002a) and Pasternak and Mühlhaus. (2005) and the references therein, and it is not considered here.

In this work, after an identiﬁcation of the drawbacks of the Cosserat model, we propose a model based on a continuum formulation that correctly describes the wave propagation in granular media in the long wavelength limit. In Section 2 and in order to get a reference for the comparison of the continuum models, the general theoretical evaluation of the bulk wave propagation in a granular assembly using a discrete model, which follows the derivation in Merkel et al. (2010) , is presented. The different interactions between the beads due to contact forces and torques are discussed. In Section 3 , the dispersion relations of the bulk eigen- modes propagating in a Face-Centered Cubic (FCC) structure along the x-axis are derived using the discrete model. The long wavelength approximations of the dispersion relations, which can be directly compared with the predictions of the continuum models, are then derived. Two cases of contacts between the beads are considered. In the ﬁrst one, the contacts are considered without solid bridges between the beads and the surface roughness is negligible; this case is called the frictionalcase corresponding to normal and sliding resistant contacts. In the second one, the contacts between the beads are considered with solid bridges and this case is called the rolling andtwistingresistant case. In Section 4 , the dispersion relations of the discrete model are compared to those obtained trough numerical simulations of wave propagation in a FCC structure. In Section 5 , the macroscopic continuum models are derived from the microscopic relations of the discrete model following the homogenization techniques proposed in Suiker et al. (2001a ; 2001b ). In Section 6 , the dispersion relations of the bulk eigen- modes in a FCC structure are derived with the Cosserat model. The

problems and drawbacks of this model are then discussed. Finally in Section 7 , the enhanced micropolar model is presented and it is shown that its approximations for small wavenumber are exactly equal to those of the discrete model in both frictional, rolling and twisting resistant cases.

2. Descriptionoftheproblem,startingfromthediscretetheory

An assembly of monodisperse beads is considered, all of them being composed by the same material. The diameter of the beads is a, the mass density of the material constituting the beads is

ρ

b,

its Poisson’s ratio is

ν

. The mass of one bead with homogeneous density is mb =

πρ

ba3/6 , its moment of inertia is Ib =mba2/10 .

The problem is considered in Cartesian coordinates with unit vectors

(

xˆ _,yˆ _,ˆ z

)

. The position of a bead

α

is deﬁned by the vector Rα. A local coordinate system ( n, s, t) at the level of the surface of contact between two beads is deﬁned: The unit vector n, normal to the surface of contact between two beads

α

and

β

, is deﬁned as Merkel et al. (2010) and Chang and Gao (1995a ); 1995b ) n=

(

Rβ− Rα

)

/

|

Rβ− Rα

|

=cos

φ

xˆ+sin

φ

cos

θ

yˆ+sin

φ

sin

θ

zˆ

(

Rβ− Rα

)

/a, (1)

where it is assumed that the static and dynamic overlaps between the particles are negligible compared to their diameter,

φ

₌ arccos

(

n· ˆx

)

,

θ

=arccos

(

n· ˆy/sin

φ

)

if

φ

= 0 and

θ

=

φ

if

φ

=0 or

π

( Chang and Gao, 1995a ). The two unit vectors s and t, which are in the contact plane, are deﬁned as

s=

∂

n/

∂

φ

=− sin

φ

xˆ+cos

φ

cos

θ

yˆ+cos

φ

sin

θ

zˆ,

t=n× s=− sin

θ

yˆ+cos

θ

zˆ. (2)

The inﬁnitesimal displacement of bead

α

is uα, its inﬁnitesimal angular rotation is wα. The dispersion relations are deduced below from the equations of motion for translation

m b

∂

2_uα

∂

t 2 = β Fβα, (3) and rotation I b

∂

2_wα

∂

t 2 = β Mβα+1 2 β Dβα× Fβα, (4)

where the summation is over all the beads

β

in contact with the bead

α

and the branch vector is

Dβα=Rβ− Rα. (5)

The direct use of the branch vector in Eq. (5) in the equation of motion for rotation in Eq. (4) is only valid in the case of monodisperse beads, i.e.,

|

Rβ

|

=

|

Rα

|

, see Suiker et al. (2001a ), Chang and Gao (1995a ), Chang and Gao (1995b ) and Luding (2008) for more general formulations. From the linearization of the Hertz- Mindlin contact model between two beads, the contact interactions can be modeled by using a normal and a shear stiffness KN and KS, respectively ( Duffy and Mindlin, 1956; Johnson, 1985;

Thornton and Yin, 1991; Gilles and Coste, 2003 ). It should be noticed that this excludes all the nonlinear effects from the analysis ( Nesterenko, 2001; Tournat and Gusev, 2010 ). In the case where all the monodisperse beads are composed by the same material, the ratio of shear to normal stiffness is given by

K =KS/KN =

2

(

1 −

ν

)

/

(

2 −

ν

)

. Since the overlap of the beads in contact is small, also the diameter of the surface of contact d is assumed to be small compared the diameter of the beads, e.g. d  a. Considering the projections on the x,y, z axis of the force applied by bead

β

on bead

α

as Merkel et al. (2010) and Suiker et al. (2001a )

(3)

Fig. 1. Sketch of the interactions (relative degrees of freedom) between two beads: a) due to normal resistance, b) due to sliding resistance, c) due to rolling resistance and d) due to twisting (or torsion) resistance.

×

u β_j − uαj+ 1 2

ε

jklD βα k

(

w βl +w αl

)

=

K Nn in j+K S

(

s is j+t it j

)

u β_j − uαj+ 1 2

ε

jklD βα k

(

w βl +w αl

)

= K _{i j}βα

δ

βα_j , (6) where K_{i j}βα=KNninj+KS

(

sisj +titj

)

and

δ

βαj =uβj − uαj+

ε

jklDβαk

(

wβl +wαl

)

/2 are the stiffness matrix and the relative

displacement, respectively, and

ε

_ijkis the permutation symbol. While the unit of the siffness constants KN and KS is the New-

ton per meter (N/m), the bending stiffness GBand the torsion stiff-

ness GT of the contact are considered here as linear stiffness con-

stants and have the unit Newton meter per radians (Nm/rad) and they include a length squared, i.e., GT =KTa2and GB =KBa2. A con-

tact model including rolling (or bending) and twisting (or torsion) resistances in agreement with experimental results is beyond the scope of this paper, see Luding (2008) , Brendel et al. (2011) and Fuchs et al. (2014) for more details. From the equivalent bending stiffness GB and torsion stiffness GTof the contact, the torque ap-

plied by bead

β

on bead

α

is

M βα_i = G Tn in j

(

w βj − wαj

)

+G B

(

s is j+t it j

)(

w βj − wαj

)

, =

G Tn in j+G B

(

s is j+t it j

)

w β_j − wαj

=G βα_{i j}

θ

_jβα, (7)

where Gβα_{i j} =GTninj +GB

(

sisj +titj

)

and

θ

βαj =wβj − wαj are the

rotational stiffness matrix and the relative rotation, respectively. From Eqs. (6) and (7) , it is possible to sketch the degrees of freedom involved in the different interactions as displayed in Fig. 1 . The terms depending on uβ_j _{− u}α_j are the classical translational normal and shear interactions, which come from the normal and sliding resistance, respectively. The terms depending on wβ_j − wαj

are the moment interactions, which are related to the torsion (or twisting) resistance parallel to the vector n and the rolling resistance in the plane composed by the vectors s and t. The term depending on wβ_l ₊wα_l is the rotational contribution to the translational surface shear displacement, which activates the sliding resistance. The bead rotation is here transformed into a relative displacement in the plane of the surfaces of the spheres through a vectorial cross product. Considering a frictional assembly without solid bridges between the beads, with a negligible roughness of the contact surface, and due to the small size of the contact surface between the beads, the equivalent bending and twisting stiffness are assumed to be negligible, e.g. Gb =Gt 0 , leading to M_iβα 0 .

The bending and the twisting stiffness become non negligible in a rolling and twisting resistant assembly.

3. DispersionrelationsintheFCCstructureinthex=x1 direction

This section is dedicated to the discrete model and its long wavelength limit.

The FCC structure is a vertical packing of hexagonal contacting layers A, B and C, which are in the closest position related to each other in an ABCABC... sequence. It can also be seen as a vertical packing of cubic contacting layers A and B, which are in the closest position related to each other, in an ABAB... sequence ( Ashcroft and Mermin, 1976 ). Without loss of generality, the x=x1axis is chosen

to be perpendicular to the cubic layers, the y=x2and z=x3 axes

are in the plane and along the cubic layers. The position of the beads can be found using integer indices

α

1,

α

2and

α

3with

Rα1,α2,α3=_a_[

α

1xˆ/ √

2+

(

α

2+

α

1/ 2

)

yˆ+

(

α

3+

α

1/ 2

)

zˆ]. (8)

The FCC structure is transversely isotropic in the plane perpendicular to the x direction. With Eqs. (3) –(7) and after a plane wave substitution u₌Uei(ωt−k·x)_and_w₌_W_ei(ωt−k·x)_,_{the dynamical matrix}

is built. Its eigenvalues give the dispersion relations of the bulk modes that propagates within the crystal with their corresponding eigenvectors X₌

(

Ux,Uy,Uz,Wx,Wy,Wz

)

.

3.1. Frictionalcase

Considering only the stiffnesses KN and KS and the wave vector k in the x direction, a pure longitudinal L mode propagates, whose eigenvector is X=

(

Ux,0 ,0 ,0 ,0 ,0

)

and has the dispersion relation

ω

2 L= 8 m b

(

K N+K S

)

sin2

κ

, (9) where

κ

=kxa/

(

2 √

2

)

is the normalized wavenumber,

κ

=

π

/2 is the minimal wavelength propagating. A pure rotational mode R exists, whose eigenvector is X=

(

0 ,0 ,0 ,Wx,0 ,0

)

and has the disper-

sion relation

ω

2 R= 20K S m b [1+cos2

_κ

_]_. ₍₁₀₎

Two degenerate Transverse–Rotational (TR) and two degenerate Rotational–Transverse (RT) modes propagate, whose eigenvectors are X₌

(

0 _,Uy,0 ,0 ,0 ,Wz

)

and X=

(

0 ,0 ,Uz,0 ,Wy,0

)

and have the

dispersion relation

ω

2 TR,RT= 1 m b

2K Nsin2

κ

+9K Scos2

κ

+11K S ±

4K 2 Nsin 4

_κ

₊_cos2

_κ

₍

₃₉₉_K2 S − 4K NK S

)

− sin2

₍

₂

_κ

₎₍

₁₂₁_K2 S/ 4+21K NK S

)

+K S2+4K NK S

1/2

, (11)

where the minus and plus signs correspond to the TR and RT modes, respectively. The small wavelength cut-off frequencies of all the modes for

κ

=

π

/2 are

ω

c L=

8

(

K N+K S

)

m b =

42

(

K S+K N

)

πρ

ba 3 , (12)

ω

c R=

20K S m b =

120K S

πρ

ba 3 , (13)

(4)

Fig. 2. Dispersion curves in the x direction of the frictional FCC structure from the discrete model (black curves) and their approximations (red dashed curves) with ν= 0 . 3 , i.e. K ࣃ 0.82. The cyclic frequencies are normalized with the cutoff fre-

quency of the longitudinal mode ωc

L . (For interpretation of the references to color

in this ﬁgure legend, the reader is referred to the web version of this article.)

ω

c RT=2

K N+3K S m b = 2

6K N+18K S

πρ

ba 3 , (14)

ω

c TR =

10K S m b =

60K S

πρ

ba 3. (15)

The cut-off frequencies of the modes R and RT for

κ

= 0 are

ω

RT,R

(

0

)

=

40K S m b =

240K S

πρ

ba 3 . (16)

It should be noticed that

ω

RT

(

0

)

=

ω

RT

(

κ

=

π

/2

)

if KS =KN/7 , but

ω

RT is not constant between these two points; the mode RT still

propagates in this case.

For long wavelength (small wavenumber

κ

=kxa/

(

2

√ 2

)

π

/2 ), approximations of the dispersion relations can be found from a Taylor expansion

ω

L

K N+K S m b k xa =

6

(

K N+K S

)

πρ

ba k x, (17)

ω

TR

K N+K S 2m b k xa =

3

(

K N+K S

)

πρ

ba k x, (18)

ω

RT

40K S m b

1− 11 320k 2 xa 2

=

240K S

πρ

ba 3

1− 11 320k 2 xa 2

, (19)

ω

R

40K S m b

1− 10 320k 2 xa 2

=

240K S

πρ

ba 3

1− 10 320k 2 xa 2

. (20)

The scaled dispersion relations and their approximations are displayed in Fig. 2 along the segment

−X of the ﬁrst Brillouin zone. The coordinates of the points

and X are

=

(

0 ,0 ,0

)

and X =

(

π

√ 2 /a,0 ,0

)

in the ( kx, ky, kz)-space. Modeling the contacts

with a relation from the Hertz-Mindlin theory, the curves are determined only by the ratio of shear to longitudinal stiffness

K,

which only depends on the Poisson’s ratio of the material constituting the beads ( Merkel et al., 2010 ).

3.2.Rollingandtwistingresistantcase

Considering all the stiffnesses KN,KS,GT,GB and the wave vec-

tor k in the x direction, a longitudinal mode L propagates, whose

eigenvector is X=

(

Ux,0 ,0 ,0 ,0 ,0

)

ω

2 L =

8

m b

(

K N+K S

)

sin2

κ

, (21)

which is identical to the dispersion relation in Eq. (9) as expected. A pure rotational mode R exists, whose eigenvector is X=

(

0 ,0 ,0 ,Wx,0 ,0

)

ω

2 R= 20K S m b

1+cos2

κ

₊80

(

G T+G B

)

m ba 2 sin2

κ

. (22)

Two degenerate transverse–rotational TR and two degenerate rotational–transverse RT modes propagate, whose eigenvectors are

X=

(

0 ,Uy,0 ,0 ,0 ,Wz

)

and X=

(

0 ,0 ,Uz,0 ,Wy,0

)

and have the dis-

persion relation

ω

2 TR,RT = 1 m b

60G B+20G T a 2 +2K N

sin2

κ

+9K Scos2

κ

+11K S ±

4

900G 2 B+100G 2T a 4 +K 2 N− 60K NG B a 2 +600G TG B a 4 − 20K NG T a 2

sin4

κ

+K S

630G B+210G T a 2 − 21K N

sin2

(

2

κ

)

+4K S

K N− 30G B+10G T a 2

sin2

κ

+K 2 S 4

(

121cos 2

₍

₂

_κ

₎

₊₇₉₈_cos

₍

₂

_κ

₎

₊₆₈₁

₎

1/2

, (23)

where the minus and plus signs correspond to the TR and RT mode, respectively. The cutoff frequencies of all the modes for

κ

=

π

/2 can then be predicted as

ω

c L =

8K N+K S m b =

42K S+K N

πρ

ba 3 , (24)

ω

c R=

20K S+4

(

G T+G B

)

/a 2 m b =

120K S+4

(

G T+G B

)

/a 2

πρ

ba 3 , (25)

ω

c TR,RT =

10K Sa 2+12G B+4G T a 2_m b =

60K S+

(

12G B+4G T

)

/a 2

πρ

ba 3 , (26)

ω

c TR,RT =2

K N+3K S m b =2

6K N+3K S

πρ

ba 3 . (27)

If

κ

=

π

/2 ,Eq. (23) reduces to Eqs. (26) or (27) without distinction whether the sign considered is plus or minus. To discriminate between the two cutoff frequencies, the repulsion property between the mode TR and RT may be used. If Eq. (26) > Eq. (27) , then Eq. (26) is the cutoff frequency of the mode RT and vice versa. The cutoff frequencies of the modes R and RT for

κ

= 0 can then be predicted as

ω

RT,R

(

0

)

=

40K S m b =

240K S

πρ

ba 3 . (28)

For long wavelength (small wavenumber

κ

=kxa/

(

2

√

2

)

π

/2 ), the approximations of the dispersion relations can be found from a Taylor expansion

ω

L

K N+K S m b k xa =

6

(

K N+K S

)

πρ

ba k x, (29)

(5)

Fig. 3. Dispersion curves in the x direction of the rolling and twisting resistant FCC structure from the discrete model (black continuous curves), their approximations (red dotted curves) and the dispersion curves of the frictional FCC structure from the discrete model (blue dashed curves) with K = K S /K N 0 . 82 and a = 2 mm. (a)

GB = K N a 2 / 10 and G T = K S a 2 / 10 . (b) G B = K N a 2 and G T = K S a 2 . The cyclic frequencies are normalized with the cut-off frequency of the longitudinal mode ωcL . The rotational

stiffnesses G B and G T are inspired by Brendel et al. (2011) . (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this

article.)

ω

TR

K N+K S 2m b k xa =

3

(

K N+K S

)

πρ

ba k x, (30)

ω

RT

40K S m b

1+3G B+G T 16K S k 2 x− 11a 2 320k 2 x

=

240K S

πρ

ba 3

1+3G B+G T 16K S k 2 x− 11a 2 320k 2 x

, (31)

ω

R

40K S m b

1+G B+G T 8K S k 2 x− 10a 2 320k 2 x

=

240K S

πρ

ba 3

1+G B+G T 8K S k 2 x− 10a 2 320k 2 x

. (32)

In the rolling and twisting resistant case, the dispersion curves depend on the 4 stiffness KN, KS, GT, GBand the diameter of the beads

a. In Fig. 3 , the dispersion relations and their approximations are displayed in two different cases. In the ﬁrst one, the cutoff frequencies of the modes R and RT at

κ

=

π

/2 are smaller than the ones at

κ

= 0 . In the second case, the situation is opposite.

In conclusion of this section, the long wavelength approximation of the discrete model is derived in the frictional (sliding resistant) case in Eqs. (17) –(20) and in the rolling and twisting resistant case in Eqs. (29) –(32) . These equations are considered as references in the comparison with the continuum models that are derived in the following sections.

4. Numericalwavepropagation

In this section, in order to validate the theoretical predictions of the discrete model, the wave propagation is numerically simulated ( Mouraille et al., 2006 ) with a packing arranged in a FCC structure of equal sized beads. The crystal is formed by square cubic layers of 4 _{× 4} particles in the y_{− z plane,} which are stacked in the closest position to each other in the x direction in an ABAB... sequence, where the B layers are shifted in the y− z plane to be able to ﬁll the deepest holes in the A layers. The crystal is composed of 200 layers, giving a total of 3200 beads, while periodic bound- ary conditions are considered in all directions. The beads have a diameter a=2 mm and are composed of glass with a mass density

ρ

b= 20 0 0 kg m −3, a Young modulus E = 69 GPa, and a Poisson’s

ratio

ν

=0 .3 . The normal static force applied on each contact in

the initial conﬁguration is F0

N=2 .10 −3 N, which, according to the

Hertz law ( Merkel et al., 2010; Duffy and Mindlin, 1956; Johnson, 1985; Thornton and Yin, 1991; Gilles and Coste, 2003 ) where K N=

(

3aFN0/ 8

)

1 3E 23

(

1−

ν

2

)

−23, (33) and

δ

0=2

3 4

(

1−

ν

2

₎

_F0 N E

a/ 2

2/3 , (34)

gives a normal stifness KN =2 .051 .10 5 N m −1and an initial inter-

particle relative displacement

δ

0 = 1 .463 .10 −8 m. The wave is ex-

cited on a single layer (here, the second layer) by imposing, depending on the mode to excite, a velocity, a rotation angle and a rotation/angular velocity. To excite the L mode (planar compres- sive P-wave) of propagation in x-direction, all the particles in the second layer get a (small) velocity

v

0

x. In contrast, to excite the TR

and RT modes (shear-, or S-waves) of propagation in x-direction, all the particles in the second layer get a (small) velocity

v

0

y or

v

0z.

Finally, to excite the R mode, all the particles of the second layer get a rotation/angular velocity w0

x and a rotation angle q0x, both

parallel to the x direction of propagation. The dispersion relations are found by computing a two-dimensional Fast Fourier Transform (FFT) in time and space using the velocities and angles per layer, sampled in space at every layer and in time at every

t = 10 −6 _s

with 2048 points.

The dispersion relations for the ”frictional” case with tangential stiffness 1 are plotted in Fig. 4 . Depending on the different types of excitation, the dispersion relations reﬂect the respective modes and the numerical simulations correspond exactly to those of all the modes predicted by the discrete model.

The dispersion relations of the TR and RT modes in the rolling and twisting resistant case are plotted in Fig. 5 and also show and conﬁrm the good agreement between the numerical simulation results and the predictions of the discrete model. It should be noticed that the complete dispersion curves are found from one single numerical simulation, since the agitation is a discontinuous change in velocity and thus activates all frequencies, which leads to a very clean signal for the perfect granular crystal structure used here. The introduction of, even weak, polydispersity in the diam- eters of the beads affects the wave propagation and, as a result,

1 The tangential elasticity K

S is active and the coeﬃcient of friction is set to a

very large value in order to avoid slip at the contacts and such as to ensure the propagation of an elastic wave.

(6)

Fig. 4. Dispersion relations from numerical wave propagation in the FCC structure with normal stiffness K N = 2 . 051 . 10 5 N m −1 and tangential stiffness K S = K K N =

1 . 689 . 10 5 N m −1 , i.e. _K_{≈ 0.82, where the inset shows the magnitude of the two-dimensional FFT in a logscale with arbitrary reference. The dispersion relations from}

the discrete mode (dashed green curves) correspond to the solid lines in Fig. 2 , with activated (a) L mode, (b) TR and RT modes, and (c) R mode. Note that in contrast to Fig. 2 , the vertical axis is not scaled. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 5. Dispersion relations from numerical wave propagation in the rolling and twisting resistant FCC structure with K N = 2 . 051 . 10 5 N m −1 ,

KS = K K N = 1 . 689 . 10 5 N m −1 , G B = K N a 2 / 10 = 8 . 20 . 10 −2 N m/rad and

GT = K S a 2 / 10 = 6 . 76 . 10 −2 N m/rad, with the magnitude of the two-dimensional FFT

in logscale with arbitrary reference. Dispersion relations (green curves) from the discrete model correspond to the solid lines in Fig. 3 . Only the TR and RT modes are shown. The dispersion relation of the L mode is identical to the one of the frictional case. The R mode is very slowly propagating (almost ﬂat dispersion curve in this case) and requires a much longer time of simulation. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

only the small wave number part of the dispersion relation (long wavelengths) can be retrieved ( Mouraille and Luding, 2008 ).

In summary, these simulations conﬁrm/validate the discrete model, on which we continue by simplifying and trying to iden- tify a simpler model for large wavelengths.

5. Fromthediscretedescriptiontothecontinuumdescription

In this section, the macroscopic stress and couple stress tensors are derived through a micro–macro transition of the applied force and torque between two beads of the discrete model in Eqs. (6) and (7) . It follows the derivations in Suiker et al. (2001a ). In Section 5.1 , the discrete displacements and rotations are transformed into continuum ﬁelds. The order of the Taylor expansion is chosen to derived a second-order gradient micropolar model. In Section 5.2 , the displacement gradient is inserted into the force equation in Eq. (6) . In Section 5.3 , the particle rotation is inserted

into the torque equation in Eq. (7) . In Section 5.4 , the Cauchy stress tensor is derived. In Section 5.5 , the couple stress tensor is derived. The Cauchy stress and the couple stress tensors form the constitutive equation of the second-order gradient micropolar model. Finally in Section 5.6 , the constitutive equations of the classical, second-order gradient and Cosserat model of elasticity are retrieved from the second-order gradient micropolar model by set- ting some of the tensors of this latter model to zero.

5.1. Expansionofthedisplacementandangularrotation

The discrete displacements and rotations are transformed into continuum ﬁelds using a Taylor expansion (the choice of the order of the expansion will be explained below)

u β_i = u α_i +D βα_j

∂

u αi

∂

x j +1 2D βα j D βαk

∂

2_uα i

∂

x j

∂

x k +1₆D βα_j D βα_k D βα_l

∂

3_uα i

∂

x j

∂

x k

∂

x l +· · · w β_i = w α_i +D βα_j

∂

w α i

∂

x j + 1 2D βα j D βα k

∂

2_wα i

∂

x j

∂

x k +· · · (35)

The displacement gradient

∂

uα_i_/

∂

xj =uα_{i, j}+uα_[_{i, j]} can be decom-

posed into its symmetric part (which is the classical strain) u α_i_{, j}=

α i j =

αji= 1 2

∂

u α_i

∂

x j +

∂

u αj

∂

x i

, (36)

and its antisymmetric part u α_[_i_{, j]} =−uα[j,i]= 1 2

∂

u α_i

∂

x j −

∂

u αj

∂

x i

. (37)

The particle rotation w α_i =

α

i +

αi (38)

is composed of two components, the particle spin

α_i and the vorticity

_iα. The particle spin is independent of the displacement ﬁeld, whereas the vorticity is completely generated by the antisymmetric part of the displacement gradient according to

α i =− 1 2

ε

i jku α[j,k] (39) and reciprocally u α_[_{i, j]}=−

ε

i jk

kα. (40)

(7)

The vorticity is here the local rotation of the assembly of particles, i.e., the rigid body rotation or the background rotation. In order to have the same order c of derivatives in the model, the order of derivative of u should be one order higher than the order of derivative of w as

{

c,c− 1

}

as detailed in Suiker et al. (2001a ). In order to get the second-order gradient micropolar constitutive equations, the order of derivatives is chosen to be {3, 2}.

5.2. Relativedisplacementandappliedforce

Putting Eq. (35) in the relative displacement

δ

βα_i in Eq. (6) , yields

δ

βα i = D βαj

∂

u α_i

∂

x j +1 2D βα j D βαk

∂

2_uα i

∂

x j

∂

x k +1 6D βα j D βαk D βαl

∂

3_uα i

∂

x j

∂

x k

∂

x l +1 2

ε

i jkD βα j

2w α_k+D βα_l

∂

w α k

∂

x l +1 2D βα l D βα p

∂

2_wα k

∂

x l

∂

x p

. (41)

Inserting the identities (37) and (38) , Eq. (41) becomes

δ

βα i = D βαj

(

i jα+u α[i, j]

)

+ 1 2D βα j D βαk

∂

α i j

∂

x k +

∂

u α_[_{i, j]}

∂

x k

+1 6D βα j D βα k D βα l

∂

2

α i j

∂

x k

∂

x l +

∂

2_uα [i, j]

∂

x k

∂

x l

+1 2

ε

i jkD βα j

2

(

α_k+

α k

)

+D βαl

∂

α k

∂

x l +

∂

kα

∂

x l

+1 2D βα l D βα p

_∂

2

α k

∂

x l

∂

x p+

∂

2

α k

∂

x l

∂

x p

. (42)

With the help of Eqs. (39) and (40) , it becomes

δ

βα i = D βαj

i jα+ 1 2D βα j D βαk

∂

α i j

∂

x k +1 6D βα j D βαk D βαl

∂

2

α i j

∂

x k

∂

x l −1 12D βα j D βαk D βαl

∂

2_uα [i, j]

∂

x k

∂

x l +1 2

ε

i jkD βα j

2

α_k+D βα_l

∂

α k

∂

x l +1 2D βα l D βαp

∂

2

α k

∂

x l

∂

x p

. (43)

From Eqs. (43) and (6) , the force applied by a bead

β

on a bead

α

is written as F _iβα= K _{i j}βα

δ

βα_j =K _{i j}βαD βα_k

α_jk+1 2K βα i j D βαk D βαl

∂

α jk

∂

x l +1 6K βα i j D βαk D βαl D βαm

∂

2

α k j

∂

x l

∂

x m −1 12K βα i j D βαk D βαl D βαm

∂

2_uα [j,k]

∂

x l

∂

x m+ K _{i j}βα

ε

jklD βαk

αk +1 2K βα i j

ε

jklD βαk D βαm

∂

α l

∂

x m +1 4K βα i j

ε

jklD βαk D βαm D βαp

∂

2

α l

∂

x m

∂

x p. (44)

5.3. Relativeangularrotationandappliedtorque

Inserting Eq. (35) in the relative angular rotation

θ

βα_j in Eq. (7) , the relative angular rotation becomes

θ

βα i =D βαk

∂

w α_j

∂

x k + 1 2D βα k D βαl

∂

2_wα j

∂

x k

∂

x l. (45)

Inserting Eqs. (38) and (39) in Eq. (45) , the applied torque in Eq. (7) becomes M βα_i = G i j

D βα_k

∂

α j

∂

x k +1 2D βα k D βαl

∂

2

α j

∂

x k

∂

x l +1 2

jnmD βα k

∂

u α_[_m_,n_]

∂

x k +1 4

jnmD βα k D βαl

∂

2_uα [m,n]

∂

x k

∂

x l

. (46)

5.4.Cauchystresstensor

The Cauchy stress tensor for a particle

α

in contact with N other particles

β

can be written as ( Suiker et al., 2001a )

σ

i j = 1 V cell N β=1 F _jβαD βα_i , (47)

where Vcell is the volume of the unit cell which is composed of

two beads, see Section 6 . Combining Eqs. (44) and (47) , the stress tensor can be expressed as

σ

hi= 1 V cell N β=1

K _{i j}βαD βα_k D βα_h

α_jk+1 2K βα i j D βαk D βαl D βαh

∂

α jk

∂

x l +1 6K βα i j D βαk D βαl D βαm D βαh

∂

2

α k j

∂

x l

∂

x m −1 12K βα i j D βαk D βαl D βαm D βαh

∂

2_uα [j,k]

∂

x l

∂

x m+ K _{i j}βα

ε

jklD βαk D βαh

αk +1 2K βα i j

ε

jklD βαk D βαm D βαh

∂

α l

∂

x m +1 4K βα i j

ε

jklD βα_k D βαm D βαp D βα_h

∂

2

α l

∂

x m

∂

x p

,

σ

i j = C i jkl

kl+D i jklm

∂

_∂

_xkl m+ E i jklmn

∂

2

kl

∂

x m

∂

x n+ F i jklmn

∂

2_u [k,l]

∂

x m

∂

x n +M i jm

m+N i jmp

∂

m

∂

x p + P i jmpq

∂

2

m

∂

x p

∂

x q. (48)

The constitutive tensors are therefore deﬁned as C i jkl= 1 V cell N β=1 K βα_jk D βα_l D βα_i , (49) D i jklm= 1 2V cell N β=1 K βα_jk D βα_l D βαm D βαi , (50) E i jklmn= 1 6V cell N β=1 K βα_jk D βα_l D mβαD βαn D βαi , (51) F i jklmn=− 1 12V cell N β=1 K βα_jk D βα_l D mβαD βαn D βαi , (52) M i jm= 1 V cell N β=1 K βα_jk

ε

klmD βαl D βαi , (53) N i jmp= 1 2V cell N β=1 K βα_jk

ε

klmD βαl D βαp D βαi , (54) Pi jmpq = 1 4V cell N β=1 K βα_jk

ε

klmD βαl D βαp D βαq D βαi . (55)

(8)

For a statistically homogeneous material, which is assumed to be the case for granular materials, the constitutive tensors have to be centrally symmetrical ( Suiker et al., 2001a; Suiker and de Borst, 2005; Chang and Gao, 1995b; dell’Isola et al., 2009 ); the tensor should have an even order. The tensors Dijklm and Nijkl do not be-

have centrally symmetrical, thus we set Di jklm =Ni jkl =0 , so that

the stress tensor reduces to

σ

i j= C i jkl

kl+E i jklmn

∂

2

kl

∂

x m

∂

x n+ F i jklmn

∂

2_u [k,l]

∂

x m

∂

x n + M i jm

m+Pi jmpq

∂

2

m

∂

x p

∂

x q. (56)

The different elements of this stress tensor can be related to the different interactions between the grains. The tensors Cijkl and

E_ijklmn are related to the classical translational normal and shear interactions with the ﬁrst order and second order gradient of displacement, respectively. The tensor Fijklmn is related to the gradi-

ent of vorticity that participate to the rotational contribution to the translational surface shear interaction through the sliding resistance, i.e., the rotation activate the shear stiffness KS. The tensors

M_ijmand P_ijmpqare related to the particle rotation interaction with the ﬁrst order and second order gradient rotation, respectively, and participate to the rotational contribution to the translational surface shear interaction through the sliding resistance.

5.5.Couplestresstensor

Following the reasoning of the Cauchy stress tensor, the couple stress tensor can be written as ( Suiker et al., 2001a )

μ

i j= 1 V cell N β=1 M βα_j D βα_i . (57)

Combining Eqs. (46) and (57) , the couple stress tensor is expressed as

μ

i j = 1 V cell N β=1

G jkD βαl D βαi

∂

α k

∂

x l + 1 2G jkD βα l D βαp D βαi

∂

2

α k

∂

x l

∂

x p +1 2G jk

knmD βα l D βαi

∂

u α_[_m_,n_]

∂

x l +1 4G jk

knmD βα l D βαp D βαi

∂

2_uα [m,n]

∂

x l

∂

x p

, (58) = Q i jkl

∂

_∂

k x l + R i jkl p

∂

2

k

∂

x l

∂

x p+ S i jnml

∂

u [m,n]

∂

x l + Ti jnml p

∂

u [m,n]

∂

x l

∂

x p. (59)

The constitutive tensors are therefore deﬁned with Q i jkl= 1 V cell N β=1 G jkD βαl D βαi , (60) R i jkl p= 1 2V cell N β=1 G jkD βαl D βαp D βαi , (61) S i jnml= 1 2V cell N β=1 G jk

knmD βαl D βα i , (62) Ti jnml p= 1 4V cell N β=1 G jk

knmD βαl D βα p D βαi . (63)

Similarly to the case of the tensors Dijklm and Nijkl, the tensors

Rijklpand Tijnmlp do not behave centrally symmetrical, thus Ri jkl p =

T_{i jnml p}= 0 . The couple stress tensor then becomes

μ

i j= Q i jkl

∂

_∂

k x l + S i jnml

∂

u [m,n]

∂

x l . (64)

The tensor Q_ijkl is related to the particle rotation interaction that activates both twisting and rolling resistances. The tensor Sijnml is

related to the gradient of vorticity and also activates both twisting and rolling resistances.

5.6. Relationbetweentheclassical,secondordergradientand Cosseratmodels

The stress tensor in Eq. (56) and the couple stress tensor in Eq. (64) can be compared to those of others elastic models. It should be noticed that the different contributions to the stress tensor due to the symmetric strain

ij, the antisymmetric strain u[i,j]

and the particle spin

_i, are here discriminated and dealt with separately. The forms of the different elastic tensors involved here are different from those found when no such discrimination is performed ( Jenkins, 1990; Jenkins and Ragione, 20 01; 20 03; Ragione and Jenkins, 2009; Ragione and Magnanimo, 2012 ). The tensor Cijkl

refers to the classical elastic tensor. In order to take into account all the symmetries of the tensor of deformation, the tensor C_ijklcan be determined with ( Ashcroft and Mermin, 1976 )

C i jkl = 1 4V cell N β=1 K βα_jk D βα_l D βα_i +K _ikβαD βα_l D βα_j +K βα_jl D βα_k D βα_i +K _ilβαD βα_k D βα_j . (65)

By exluding the rotational degrees of freedom from the analysis, the tensors Mijm,Pijmpq,Qijkland Sijmnlare then set to zero, the cou-

ple stress tensor vanishes and the stress tensor reduces to

σ

i j=C i jkl

kl+E i jklmn

∂

2

kl

∂

x m

∂

x n+ F i jklmn

∂

2_u [k,l]

∂

x m

∂

x n, (66)

which is similar to the second-order gradient model ( Chang and Gao, 1995b; Suiker et al., 2001a ).

By keeping only the ﬁrst-order gradient in the analysis, which restricts the order of derivatives of u and w to {2, 1}, the tensors Eijklmn,Fijklmnand Pijmpqare then set to zero. The stress tensor and

the couple stress tensor then reduce to

σ

i j = C i jkl

kl+M i jm

m, (67)

μ

i j = Q i jkl

∂

_∂

_xk l +S i jnml

∂

u [m,n]

∂

x l , (68)

which are the constitutive relations of the Cosserat model. An important point is lying here. Considering the frictional case, where GB =GT =0 , the couple stress tensor is then excluded from

the analysis, i.e.,

μ

i j = 0 . The stress tensor in Eq. (67) can be found

using the following expansion of the discrete displacements and rotations

u β_i = u α_i +D βα_j

∂

u αi

∂

x j

, w β_i = w α_i, (69)

which restricts the order of derivatives of u and w to {1, 0}. Thus, the Cosserat model does not contain a derivative over space coordinates in the rotational contribution to translational shear interaction in the stress tensor, which is represented by the term Mijm

m

in Eq. (67) and involves the shear stiffness KS. In others words, a

dependence on space coordinates of the particle rotations is not taken into account in the prediction of the stress tensor. Although this activates the sliding resistance at the surface of the particles as it can be seen in Fig. 1 (a), it can be expected that the absence of a derivative over space coordinates can lead to an inconsistent

(9)

representation of the interaction due to the activation of the sliding resistance by particle rotations, and thus, to inconsistent results as it will be further discussed in the next section. On the opposite side, it is important to note that the couple stress tensor

μ

_ijfor the interactions due to rolling resistance, i.e., involving GB , and due

to twisting resistance, i.e., involving GT contains a derivative over

the space coordinates and remains the same in the Cosserat and in the second-order gradient micropolar model. In the next section, a Cosserat model is used to predict the bulk mode propagation in the FCC structure and is compared to the long wavelength approximations of the discrete model.

6. Dispersionrelationsofthebulkmodepropagatinginthe FCCstructureusingtheCosseratmodel:limitsofthemodel

The dispersion relations of the bulk mode propagating in the FCC structure are predicted using the stress tensor in Eq. (67) and the couple stress tensor in Eq. (68) , which correspond to the Cosserat model. The wave propagation is described by the equation of motion in translation

ρ ∂

2u i

∂

t 2 =

∂σ

ji

∂

x j =C i jkl

∂

_∂

_xkl j +M jim

∂

m

∂

x j , (70)

where

ρ

₌

ηρ

_b and

η

is the volume fraction of the structure, which is

η

=

π

√2 /6 0 .74 in the FCC structure. The equation of rotational motion is J

∂

2_w i

∂

t 2 =

∂μ

ji

∂

x j +

ε

ipq

σ

pq= Q jikl

∂

2

k

∂

x j

∂

x l +S jinml

∂

2_u [m,n]

∂

x j

∂

x l +

ε

ipq

(

C pqrs

rs+M pqt

t

)

. (71)

where J=√2 Ib/a3is the density of moment of inertia, which is de-

duced from the number n of beads in a cube with an edge-length a ( Merkel et al., 2011 ).

In order to evaluate Vcell, The number n of beads in a cube with

an edge-length a is found from the volume fraction

η

and the volume of one bead V_beadas

n =

η

a 3/V bead= √

2

π

a 3/ [6

(

4

π

(

a/ 2

)

3/ 3

)

]=√2, (72)

i.e., there are √2 beads in a cube of volume a3_{. In others words,}

the deﬁned volume of the unit cell, which should include 2 beads, is therefore Vcell =a3

√ 2 .

As already noticed in Section 5.6 , the elastic tensor Cijkl refers

to the classical elastic tensor ( Suiker et al., 2001a ), which relates the symmetric part of the strain

ij to the stress tensor

σ

ij. There-

fore, the symmetries of the classical elastic tensor are assumed to also apply here. By combination of Eqs. (5) , (6), (8) and (65) , the components of the matrix CIJ, which gives the 36 independent con-

stants of the tensor Cijkl ( Royer and Dieulesaint, 20 0 0 ), in the FCC

structure are CIJ= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ √ 2(KN+KS) a KN− KS a√2 KN− KS a√2 0 0 0 KN− KS a√2 5KN+3KS 2a√2 KN− KS 2a√2 0 0 0 KN− KS a√2 KN− KS 2a√2 5KN+3KS 2a√2 0 0 0 0 0 0 KN+KS a√2 0 0 0 0 0 0 KN+KS a√2 0 0 0 0 0 0 KN+3KS 2a√2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (73)

where ( I, J) _∈ {11, 22, 33, 12, 13, 23}. By combination of Eqs. (5) , (8) and (53) , the component of the tensor Mijmreduces to

M i jm=−

ε

i jm 2√2K S

a . (74)

The tensors Qjikl and Sjinml are more complex and are predicted

with the combination of Eqs. (5) , (8) with (60) and (62) , respectively, see Section 6.2 .

6.1. Frictionalcase

In this section, a frictional assembly is considered. Thus, GT =

GB =

μ

i j = 0 and the equation of rotational motion reduces to

J

∂

2_w i

∂

t 2 =

ε

ipq

(

C pqrs

rs+M pqt

t

)

. (75)

Therefore, this case does not contain all the relative degrees of motion of the Cosserat model and can be called a reduced Cosserat model ( Kulesh et al., 2009 ).

6.1.1. Longitudinalwavepropagationindirectionx1 =x

The axis perpendicular to the cubic layers is the four-fold axis of high symmetry. The calculations here consider only plane wave propagation with the polarization of the wave being pure. Thus, without loss of generality, the calculations do not require the use of scalar

φ

and vector

ψ

potentials such as u =

∇

·

φ

+

∇

∧

ψ

( Royer and Dieulesaint, 20 0 0 ); it can be checked that the use of these two potentials will lead to the same results. Writing the equations of motion (70) and (75) considering that the only non- zero component of displacement is the component u1 gives

ρ ∂

2u 1

∂

t 2 =C 11

∂

2_u 1

∂

x 2 1 , J

∂

2_w i

∂

t 2 =0. (76)

As expected, there is no contribution of the rotational degrees of freedom in the case of a pure longitudinal wave propagation. Con- sidering a plane wave solution u1 =U1ei(ωt−kxx), the dispersion re-

lation for the longitudinal wave is found

ω

L=

√ 2

(

K N+K S

)

a

ρ

k x=

6

(

K N+K S

)

π

a

ρ

b k x, (77)

which is equal to the long wavelength approximation of the discrete model in Eq. (17) . In order to compare the prediction from the Cosserat model with the approximations of the discrete model, the continuum mass density

ρ

is replaced by the mass density of the materials constituting the beads

ρ

band the volume fraction

η

.

6.1.2. Purerotationalwavepropagationindirectionx1 =x

Considering a pure rotational wave propagating along x1direc-

tion, where u=0 ,Eqs. (70) and (71) become M 1im

∂

m

∂

x 1 =0 ⇒ m =1, J

∂

2_w 1

∂

t 2 =

ε

1jkM jk1

1=

(

M 231− M321

)

1=−2M 321

1. (78)

where the index m of the particle spin

mshould be equal in both

equations. With wi =

i +

i =

i−

ε

i jku[j,k]/2 and u[j,k] =0 , Eq.

(78) then becomes J

∂

2

1

∂

t 2 =− 4√2 a K S

1, (79)

which is not an equation describing a wave propagation. With the plane wave substitution

1 =P1ei(ωt−kxx), the dispersion relation

for the pure rotational mode is found

ω

R=

4√2 aJ K S=

240K S

πρ

ba 3, (80)

which corresponds to Eq. (20) only for kx =0 . This dispersion re-

lation does not depend on k, the pure rotational mode does not propagate in the x1 direction, which is a consequence of Eqs.