## Abstract

This article discusses the derivation of continuum models that can be used for modelling the inhomogeneous mechanical behaviour of granular assemblies. These so-called kinematically enhanced models are of the strain-gradient type and of the strain-gradient micro-polar type, and are derived by means of homogenizing the micro-structural interactions between discrete particles. By analysis of the body wave dispersion curves, the enhanced continuum models are compared to corresponding discrete lattice models. Accordingly, it can be examined up to which deformation level the continuum models are able to accurately describe the discrete particle behaviour. Further, the boundary conditions for the enhanced continuum models are formulated, and their stability is considered. It is demonstrated how to use the body wave dispersion relations for the assessment of stability.

## 1. Introduction

The modelling of the mechanical behaviour of granular materials can be globally divided into two categories. Firstly, there are the discrete models, where the equilibrium conditions, the kinematic conditions and the constitutive behaviour are formulated for each particle contact. Secondly, there are the continuum models, where the equilibrium conditions, the kinematic conditions and the constitutive behaviour are formulated for an assembly of particles using the continuum concepts of stress and strain. The advantage of discrete models as opposed to continuum models is that the inhomogeneous effects at the particle level can be described more accurately. However, in the case of engineering applications, the number of particles is very large, which in turn results in a large number of equations to be solved for a discrete system.

In order to study the particle effects while limiting the number of equations to be solved, so-called *micro-structural continuum models* have been developed. These models draw upon the homogenization of the discrete particle behaviour, causing the model coefficients to be related to the micro-scale properties, such as contact stiffness, particle density and particle radius. The homogenization of granular media has initially led to models for regular granular packings (Duffy & Mindlin 1957; Deresiewicz 1958) while later also models for random granular packings were developed (Digby 1981; Walton 1987; Bathurst & Rothenburg 1988; Chang 1988; Jenkins 1988). All the above studies were based on the premise of a uniform strain field, and the effect of particle rotation was also ignored. During the last two decades, the micro-structural models for granular materials have been enhanced by incorporating the influence of particle rotation. This has led to formulations of the micro-polar type (Mühlhaus & Vardoulakis 1987; Chang & Liao 1990; Chang & Ma 1992; Suiker *et al*. 1999, 2001*a*). Furthermore, inhomogeneous effects by particle displacements have been included, yielding strain-gradient type formulations (Chang & Gao 1995, 1997; Suiker *et al*. 1999, 2001*a*; Askes & Metrikine 2002; Metrikine & Askes 2002). Only recently, the inhomogeneous effects by particle displacements and particle rotations have been captured within one framework, which has resulted in strain-gradient micro-polar models (Mühlhaus & Oka 1996; Suiker *et al*. 2001*a*,*b*).

A special approach within the theory of homogenization of granular materials is the so-called *micro-structural approach* (Christofferson *et al*. 1981; Walton 1987; Chang 1988; Chang & Ma 1992; Liao *et al*. 2000; Suiker *et al*. 2001*a*,*b*). The micro-structural approach requires three relationships to be formulated, which are (i) the discrete particle kinematics expressed in terms of kinematic field variables, (ii) the particle contact forces expressed in terms of the particle assembly stress, and (iii) a constitutive law describing the contact behaviour between two particles. Naturally, the macro-level constitutive model is obtained by combining these three relationships. The current study focuses on how the constitutive formulation is influenced by the relations (i) and (ii). For reasons of simplicity, relation (iii) is assumed as linear elastic, although the micro-structural approach allows for including more advanced contact laws as well (Chang 1993; Liao *et al*. 2000).

The continuum models treated in this article are of the *strain-gradient* type and of the *strain-gradient micro-polar* type. The derivation of the latter type of model has been discussed thoroughly in a previous article (Suiker *et al*. 2001*a*), but for reasons of clarity the main steps in the derivation procedure will be repeated here. After discussing the homogenization procedure, the body wave dispersion curves for the derived continuum models are compared to the body wave dispersion curves for corresponding discrete lattice models. Hereto, the deformation level can be determined at which the continuum model fails to describe the discrete particle behaviour accurately. This accuracy level depends on the number of terms included in the Taylor series used to approximate the discrete particle kinematics. From a practical point of view, the Taylor series should be as short as possible, since the homogenization procedure becomes (much) more involved each time the Taylor series is extended.

For some enhanced continuum models, the loss of accuracy at a specific deformation level may be accompanied by the initiation of instability and the loss of uniqueness. When analysing boundary value problems, the physically irrational character of this instability requires the analyst to exclude the unstable deformation terms from the computed response. Hence, it is useful to know a priori if, and at which deformation level, this instability is initiated. Apart from analysing the second-order variation of the internal work for finding a possible instability (Hill 1957), instabilities can be also found from the analysis of the body wave dispersion curves.

Although in general a deformed material may store energy in various forms, i.e. mechanical energy, thermal energy, chemical energy and electrical energy, the current study presumes the motion of the material to involve solely mechanical changes. Correspondingly, the material state is considered to be dependent on strain (gradients) only, which, for many processes in the field of granular mechanics, is indeed the most prominent state variable.

## 2. Homogenization of a granular material

In this section, the micro-structural approach is employed for the derivation of strain-gradient continuum models. The particle contact law and the particle kinematics are formulated, followed by expressing the particle assembly stress in terms of particle contact forces. Combining these relations, and assuming the granular material to consist of spherical, equal-sized particles, provides a constitutive relationship in a closed form. The outcome is a so-called *fourth-gradient model*, which can be reduced to the classic linear elastic model and a second-gradient model by requiring specific constitutive parameters to be zero.

### (a) Micro-level particle interaction

Figure 1 shows the interaction between two convex-shaped particles, designated as *n* and *m*. The particles can displace and rotate, where initially only the effect of the particle displacements is considered. Condensing the interparticle displacement to the contact point *c*, the contact displacement reads(2.1)

By means of a contact law, the contact displacements are related to energetically conjugated contact forces. In the present study the contact law is assumed as elastic, meaning that irreversible sliding, separation and rearrangement of particles is not considered. The constitutive relation thus becomes(2.2)where the elastic contact stiffness is written as(2.3)with *Kn* the contact stiffness in the direction normal to the contact plane, and *Ks* and *Kt* the contact stiffnesses in the directions tangential to the contact plane. Further, **n**^{c}, **s**^{c} and **t**^{c} are the orthonormal base vectors at the interparticle contact point in the normal direction and the two tangential directions, respectively, see figure 2. Here, the local orthonormal reference system (**n**^{c}, **s**^{c}, **t**^{c}) may be related to a global orthonormal reference system (**e**_{1}, **e**_{2}, **e**_{3}) as(2.4)with *γ* and *ϕ* the angle of inclination and the angle of rotation (Euler angles), respectively.

For developing a continuum model that describes the mechanical behaviour of a large assembly of particles, the discrete kinematic degrees of freedom of the individual particles need to be replaced by continuous field variables, i.e.(2.5)where the lower bar denotes the continuous character of the field variable. Subsequently, the kinematics of an individual particle *m* is expressed in terms of the kinematics of its neighbouring particle *n* by employing a Taylor series(2.6)where, in agreement with figure 1, represents the branch vector that connects the centres of two neighbouring particles(2.7)

Similar to equation (2.5), the higher-order kinematic expressions in equation (2.6) have the form(2.8)

For the sake of convenience, the lower bar used for denoting the kinematic field variables will be dropped in the following. The accuracy level at which the discrete particle kinematics is approximated by a continuous field depends on the number of terms that is included in the Taylor series. In order to arrive at a fourth-order gradient continuum model, the Taylor expansion (2.6) needs to be truncated after the fifth-order term. Substituting the result into equation (2.1) gives for the contact displacement(2.9)

Trivially, the displacement gradient in equation (2.9) may be decomposed in a symmetric part and an antisymmetric part (2.10)where, within the frames of the linear deformation theory, the symmetric part represents the strain (2.11)and the antisymmetric part ensues from(2.12)

Combining equations (2.10)–(2.12) with equation (2.9) yields(2.13)

Since equation (2.13) incorporates strain contributions up to and including the fourth-order gradient term, it can be characterized as a fourth-order gradient continuum, or a fourth-gradient continuum. The fourth-gradient continuum reduces to the classic Boltzmann continuum when all terms apart from are ignored.

### (b) From micro-level to macro-level

In order to transmit the micro-structural properties to the macroscopic level, an expression for the macroscopic stress in terms of the particle contact forces needs to be derived. Starting with the equilibrium condition(2.14)the mean particle stress over the particle volume *V*^{n} may be computed as(2.15)with *X*_{i} the position vector in the *i*th direction. By invoking Gauss theorem, also known as the divergence theorem, equation (2.15) turns into a surface integral expression(2.16)where *n*_{k} is the unit outward, normal to the particle surface *S*^{n}. The discrete contact forces acting on the particle surface may be related to the surface integral in equation (2.16) as(2.17)where *N*^{n} is the total number of contacts for particle *n*, and denotes the position of particle *n*. Combining equations (2.16) and (2.17) yields(2.18)

As a next step, the mean stress for a representative volume of *M* particles is determined by taking the volume average of the local particle stresses(2.19)with *V* the volume of the particle assembly. It is obvious that the right-hand side of equation (2.19) has been constructed by substitution of equation (2.18). An alternative format of equation (2.19) can be obtained by using the fact that the branch vector, equation (2.7), may be rephrased as(2.20)which, together with , and using the fact that *two particles n* and *m* construct *one particle contact c*, yields(2.21)where *N* is the total number of particle contacts in the particle assembly. The well-known expression, equation (2.21), has been presented previously in other works on granular mechanics (Christofferson *et al*. 1981; Rothenburg & Selvadurai 1981; Chang 1988; Chang & Ma 1992; Cambou *et al*. 1995; Suiker *et al*. 2001*a*).

To compute an expression for the macroscopic constitutive behaviour, the stress formulation (2.21), the kinematic expression (2.13) and the local contact law (2.2) must be combined. For reasons of simplicity, it is initially assumed that the local contact law includes only the longitudinal stiffness *Kn*, which reduces equation (2.2) to(2.22)

Hence, the constitutive formulation becomes(2.23)where the bar on top of the Cauchy stress has been omitted for reasons of convenience. The constitutive relation above can be compactly written in tensor format as(2.24)where(2.25)

In equation (2.23), the particle-related superscript *n* of the kinematic field variables has been omitted, because the kinematic relation, equation (2.13), is assumed to be representative for all particles in the granular assembly. This statement is known as the *kinematic hypothesis*, and is often applied in homogenization procedures regarding standard granular continua (Digby 1981; Walton 1987; Chang 1988), as well as in homogenization procedures regarding enhanced granular continua (Chang & Ma 1992; Chang & Gao 1995, 1997; Mühlhaus & Oka 1996; Suiker *et al*. 2001*a*,*b*). The kinematic hypothesis enforces kinematic compatibility between the individual particles, which imposes a constraint that causes the macroscopic constitutive formulation to predict a somewhat stiffer behaviour than the true material behaviour, i.e. an overestimation (Cambou *et al*. 1995; Chang & Gao 1996; Liao *et al*. 2000). Instead of using the kinematic hypothesis, a so-called *static hypothesis* could have been employed, which relates the contact force at the particle contact uniquely to the average stress in the particle assembly volume (Cambou *et al*. 1995; Chang & Gao 1996; Liao *et al*. 2000).1 The static hypothesis commonly results in a somewhat softer material behaviour than the true material behaviour, i.e. an underestimation.

When supposing the material to be statistically homogeneous, the stiffness tensors in equation (2.25) will have a centrally symmetric character. Consequently, in equation (2.25) the stiffness tensors of uneven order become zero, which can be mathematically proven by using Weyl's theory on orthogonal invariant polynomial functions (Weyl 1946; Suiker & Chang 2000). Accordingly, the constitutive law (2.24) reduces to(2.26)

### (c) Macroscopic constitutive formulation

To derive closed-form expressions for the constitutive coefficients in equation (2.26), it is assumed that the granular material consists of equal-sized spherical particles, and behaves in an ideally isotropic fashion. The branch vector (2.7) then becomes(2.27)with *r* the particle radius. For granular assemblies consisting of a very large number of particles, the summation of an arbitrary quantity *F*^{c}(*γ*, *ϕ*) over the total number of particles *N* may be replaced by a corresponding integral form via the introduction of a density function. Here, *γ* and *ϕ* are the angle of inclination and the angle of rotation in a spherical coordinate system, see also figure 2. Assuming an isotropic distribution of interparticle contacts, the density function equals 1/4*π*, thereby satisfying the requirement(2.28)

With *N* the total number of contacts in a volume *V*, the number of contacts in the interval ranging from *γ* to *γ*+d*γ* and from *ϕ* to *ϕ*+d*ϕ* becomes (*N*/4*π*)sin *γ* d*γ* d*ϕ*. The summation of the quantity *F*^{c}(*γ*, *ϕ*) over the total number of contacts *N* may thus be rephrased as(2.29)

Applying equation (2.29) to the constitutive tensors in equation (2.26) yields(2.30)

In the case of isotropic tensors, the tensor coefficients can be derived by using Weyl's theory about orthogonal invariant polynomial functions (Weyl 1946; Suiker & Chang 2000). Accordingly, the fourth-order isotropic tensor is expressed by three terms(2.31)whereas the sixth-order isotropic tensor is expressed by 15 terms(2.32)

The eighth-order isotropic tensor is expressed by 105 terms, which will not be given here for the sake of brevity. Employing the expressions for the isotropic tensors with respect to equation (2.30), leads, after an extensive derivation procedure, to the constitutive form(2.33)for which the coefficients read(2.34)where *θ* is the package density of the granular structure,(2.35)

Expression (2.33) shows that the Cauchy stress is symmetric, and that the constitutive relation contains only two independent coefficients, which are the ‘stiffness’ *θKn* and the particle radius *r*. As a result of the incorporation of only one contact stiffness *Kn*, the continuum model (2.33) includes a restriction with respect to the Lamé constants, i.e. *λ*=*μ* (or *ν*=0.25). This constraint can be released by incorporating an additional shear contact stiffnesses in the model, which will be shown in §4.

The *fourth-order gradient model* (2.33) reverts to the classic linear elastic model when requiring the coefficients *B*_{(1)} to *B*_{(5)} to be equal to zero,(2.36)

Additionally, only requiring *B*_{(3)} to *B*_{(5)} to be zero yields the second-gradient model(2.37)which has been previously discussed in other works (Chang & Gao 1995, 1997; Suiker *et al*. 1999).

To compute the potential energy density as well as the boundary conditions for the fourth-gradient model, the weak form of the equilibrium, equation (2.14), over a volume *V* is considered(2.38)in which the body forces have been neglected. Integrating by parts and using Gauss theorem results in(2.39)where *S* is the boundary surface and *n*_{i} is the unit outward, normal to *S*. Similar to the classic continuum theory, the surface integral includes the *natural boundary data n*_{i}*σ*_{ij}, which is prescribed complementary to *essential boundary data* for *u*_{j}. However, it is emphasized that the character of the natural boundary conditions differs from that in a classic continuum, since the stress *σ*_{ij} here also relates to higher-order strain derivatives, see equation (2.33). Accordingly, the traction *n*_{i}*σ*_{ij} consists of a standard contribution and a higher-order contribution. The volume integral emerging in equation (2.39) can be further elaborated by substitution of equation (2.33),(2.40)

Integrating the higher-order terms by parts and applying Gauss theorem leads to(2.41)

Subsequently, inserting equation (2.41) into equation (2.39) yields(2.42)where the boundary tractions *t*_{j}, *p*_{ji} and *r*_{jil} satisfy the relations(2.43)with *σ*_{ij} given by equation (2.33). The expressions given by equation (2.43) include all the natural boundary data for the fourth-gradient model, as prescribed in addition to essential boundary data for *u*_{j}, *u*_{j,i} and *u*_{j,il}. The necessity of providing higher-order boundary data, (*p*_{ji}, *u*_{j,i}), (*r*_{jil}, *u*_{j,il}), is intrinsic to the appearance of strain-gradient terms in the constitutive law, equation (2.33).

It is emphasized that for three-dimensional boundary value problems the variations *δu*_{j}, *δu*_{j,i} and *δu*_{j,il} are dependent on each other, in a sense that, if *δu*_{j} is known at the boundary, so are the tangential derivatives of *δu*_{j}. In fact, only the derivatives in the direction normal to the boundary can be varied independently. This difficulty can be overcome by resolving the gradients of *δu*_{j} into a tangential gradient and a normal gradient, and subsequently applying Stokes surface divergence theorem. For further details on this procedure, the reader is referred to Mindlin (1965).

In equation (2.42), the volume integral at the left-hand side can be identified as the variation of the potential energy, , with the potential energy density. Correspondingly, the potential energy density may be formulated as(2.44)

In the next section, it will be examined how the individual terms in equation (2.44) influence the stability of the model.

## 3. Strain-gradient continua versus discrete Born–Karman lattice

In order to scrutinize at which accuracy level the discrete particle behaviour is captured by the continuum models (2.33), (2.36) and (2.37), the wave propagation characteristics of these models will be compared with those of the discrete, one-dimensional Born–Karman lattice (Brillouin 1946; Born & Huang 1954). The Born–Karman lattice consists of point masses *M*, which are spaced along a straight line at a distance *d*=2*r* from one another, see figure 3. The interaction between a mass and its neighbouring masses occurs by means of linear longitudinal springs *Kn*, where the value of the spring constant is considered to be identical everywhere in the lattice. Due to this repetitive characteristic, the lattice response can be analysed by considering only the compound of a mass *M* and the two attaching springs, which is designated as a *cell*. The elementary, one-dimensional character of the Born–Karman lattice causes this model to pre-eminently lend itself for a simple and illustrative comparison with the continuum models. For preserving generality, further in this article a more advanced lattice will be considered as well.

Some of the ideas presented in this section have been discussed previously by other investigators, for instance Kunin (1982), although the homogenization method used by him differs from the one currently discussed. This, in a sense that Kunin considers approximate forms of the Fourier image of the equations of motion for a discrete lattice, which reflect the dispersion relations for corresponding, kinematically enhanced continua, see also Born & Huang (1954) and Maradudin *et al*. (1971). The degree at which the Fourier image is approximated determines the accuracy level at which the discrete kinematics of the lattice is represented, in a similar way as the Taylor series, equation (2.6), does in the current homogenization procedure. An example of the homogenization method employed by Kunin will be given later in this article.

### (a) Dispersion relation for Born–Karman lattice

The dynamic behaviour of the Born–Karman lattice may be described by the Lagrangian *L* of an individual cell with a dimensionless coordinate *m*(3.1)where is the potential energy and is the kinetic energy. The potential energy of the cell reads(3.2)with the length variations of the springs as(3.3)

Furthermore, the kinetic energy of cell *m* follows from(3.4)in which the dot on top of the displacement designates the total derivative with respect to time, and *M* represents the mass of the cell. Combining equations (3.1)–(3.4) and subsequently applying the Lagrange equation (see for example, Landau & Lifshitz 1976)(3.5)straightforwardly provides the equation of motion for the cell *m*(3.6)

In order to satisfy equation (3.6), a solution of the harmonic form is sought(3.7)in which *ω* is the angular frequency, is the wave amplitude and *k*_{x} is the wavenumber in the *x*-direction. The wavenumber is related to the wavelength *Λ*_{x} and to the phase velocity *c*_{x} via, respectively(3.8)

Combining equation (3.7) with equation (3.6) leads to(3.9)

Replacing the exponential terms in the above equation by a corresponding sinusoidal function, followed by some re-ordering, yields the dispersion relation(3.10)in which *c*^{∞} is the velocity of the compression wave at infinite wavelength(3.11)

### (b) Dispersion relation and stability aspects for strain-gradient continua

For an adequate comparison of the wave propagation characteristics of the one-dimensional Born–Karman model with those of the continuum models, the balance of linear momentum in the continuum needs to be confined to the *x*-direction, where stress contributions related to the *y* and *z*-directions are ignored(3.12)with *ρ* the density and ,*t* the partial derivative with respect to time. For the fourth-gradient model (2.33), the normal stress in the *x*-direction reads(3.13)where insertion of equation (3.13) into equation (3.12), and using the kinematic relation (2.11), yields the equation of motion(3.14)

To satisfy this sixth-order differential equation, a wave of the harmonic type is considered(3.15)

Substitution of this expression into equation (3.14) provides the algebraic equation(3.16)which, after some rephrasing, leads to the following dispersion relation(3.17)with *c*^{P} the compression wave velocity of a classic linear elastic material,(3.18)

Invoking the constitutive coefficients, equation (2.34), turns the dispersion relation for the fourth-gradient model, equation (3.17), into(3.19)

Because for arbitrary real wavenumbers *k*_{x} the term under the square root is larger than zero, i.e.(3.20)the dispersion relation (3.17) always represents waves of the harmonic type. Consequently, equation (3.15) may be formulated as(3.21)which clearly exhibits its harmonic character with respect to time. When *B*_{(1)} to *B*_{(5)} are set to zero, equation (3.17) indeed reduces to the frequency expression for a classic linear elastic continuum(3.22)

Furthermore, neglecting the terms *B*_{(3)} to *B*_{(5)} yields the dispersion relation for the second-gradient model(3.23)

Similar to equation (3.19), this expression may be formulated as(3.24)

However, for the term under the square root it now follows that(3.25)so that equation (3.24) needs to be partitioned as(3.26)

In the equations above, the parameter *ξ* is real. Correspondingly, the parameter *ω** is imaginary, which is signified by the asterisk. Thus, the latter parameter may not be interpreted as the angular frequency. In actual fact, this yields the following partitioning of the general form, equation (3.15),(3.27)

Due to its exponential character with respect to time, equation (3.27*b*) goes to infinity when time unboundedly increases, i.e.(3.28)

The exponential growth with respect to time characterizes an unstable response, which naturally is an unrealistic property of the granular model. Consequently, when using the second-gradient model in the analysis of dynamic boundary value problems, the analyst has to be convinced that the generated deformation pattern contains only terms that satisfy the requirement of structural stability. In other words, the minimum wavelength appearing in the response needs to be larger than the critical wavelength that initiates structural instability. When following this requirement, in dynamic boundary value problems the use of a potentially unstable, granular continuum model may nevertheless provide a stable response (Chang & Gao 1997).

Instead of employing a dispersion analysis for finding the critical wavelength, this can be done also by analysing the second-order variation of the internal work, *δ*^{2}*W*^{int}, simply called the second-order internal work. The second-order internal work needs to be positive in order to guarantee structural stability (Hill 1957), i.e.(3.29)in which *δ* denotes the variation of an admissible kinematic field. For a one-dimensional bar with length *L*^{b}, a combination of equations (2.44) and (3.29) yields(3.30)

Subsequently, the following six boundary conditions are applied at the bar ends(3.31)

Here, equation (3.31*a*) reflects the classical, essential boundary condition, while equation (3.31*c*) reflects one of the two higher-order essential boundary conditions (see discussion below equation (2.43)). Further, equation (3.31*b*) is the natural boundary condition for the higher-order traction in equation (2.43*b*). In correspondence with equation (2.43*b*), for a one-dimensional bar the two boundary conditions presented by equation (3.31*b*) turn into(3.32)

When combining equation (3.32) with equation (3.31*c*), it follows that(3.33)

A displacement variation that satisfies the boundary conditions (3.31*a*,*c*) and (3.33) is2(3.34)in which *n* represents an arbitrary integer. With the use of equation (3.34), the stability condition (3.30) becomes(3.35)

In equation (3.35), the term between the parentheses governs the sign of the expression, and thus determines whether or not the stability requirement is fulfilled. For an *arbitrary bar length*, 0<*L*^{b}<∞, stability is warranted if(3.36)

Substitution of the constitutive coefficients (2.34) into equation (3.36) indeed yields the same stability requirement as equation (3.20). It is emphasized, however, that the above consideration concerns a linear elastic system. For other material characteristics, the agreement between equations (3.20) and (3.36) does not necessarily hold. In view of our previous findings, equation (3.36) illustrates that the coefficients with a minus sign, *B*_{(1)} and *B*_{(2)}, have a *destabilizing* effect on the structural response, while the other coefficients, relating to a positive sign, have a *stabilizing* effect. Furthermore, when the order of the included gradient terms increases, stabilizing and destabilizing contributions appear in an alternating manner, i.e. the zeroth-order term is stabilizing, the second-order term is destabilizing, the fourth-order term is stabilizing, and so on. The appearance of stabilizing and destabilizing terms in granular material models has also been reported by Mühlhaus & Oka (1996), who discuss the derivation of an enhanced continuum formulation by means of a Taylor approximation of the wave equations of a representative discrete system.

For the derivation of equation (3.36), a deformation pattern has been considered which consists of only one harmonic component, cf. equation (3.34). However, in general a harmonic response is represented by a multiple number of harmonic components. If the response is composed of a multiple number of harmonic components, and, at a certain deformation stage, one harmonic component violates the condition (3.36), it might be possible that the total response remains stable. Nonetheless, under most circumstances the instability of an individual harmonic component will initiate an overall instability.

### (c) Discussion of dispersion curves

In figure 4, the dispersion relations (3.10), (3.19), (3.22) and (3.24) have been visualized by using the parameter set *λ*=*μ*=80 MPa, *ρ*=1800 kg m^{−3}, *d*=2*r*=50 mm. Because the dispersion curves are symmetric with respect to *k*_{x}=0 and *ω*=0, only the positive wavenumber axis and the positive frequency axis have been plotted. In the long-wave limit, *k*_{x}→0, it is required that the continuum models and the discrete Born–Karman model give equivalent results. The velocity at infinite wavelength of the Born–Karman model, equation (3.11), has therefore been set equal to the compression wave velocity *c*^{P} of the classic elastic continuum, equation (3.18).3 It can be noticed that for an increasing wavenumber (or a decreasing wavelength), the inhomogeneous effect by the particle size become more pronounced. In fact, the granular models tend to behave more dispersively, where the dispersion curves start to deviate from the secant slope that reflects the compression wave velocity *c*^{P} of the classic elastic continuum. The intensity of this dispersion depends on the extent to which the phase velocity *c*_{x} of the harmonic waves(3.37)differs from the group velocity of the harmonic waves(3.38)

In correspondence with equations (3.37) and (3.38), in figure 4 the phase velocity is reflected by the secant slope of the dispersion curve, while the group velocity is reflected by the tangential slope. For un-absorbed travelling harmonic waves, the group velocity (3.38) equals the velocity of propagation of the wave energy (Brillouin 1946; Achenbach 1973).

The dispersion curves for the second-gradient continuum and the Born–Karman lattice show a clear similarity, in a sense that both curves first distinguish an upward trend, which after passage of a maximum frequency proceeds into a downward trend. When the downward branch of the second-gradient dispersion curve encounters the frequency axis *ω*=0, the stable harmonic response turns into an unstable, exponentially increasing response, see also equation (3.26*b*). The wavenumber at which this instability occurs equals , and thus can be called *critical*. Although the dispersion curve for the Born–Karman lattice also meets the axis *ω*=0, the sinusoidal periodicity of the dispersion relations (3.10), which characterizes the so-called *Brillouin zones* with *n* an integer (Brillouin 1946), prevents the frequency from becoming imaginary. So, the response of the Born–Karman lattice is unconditionally stable.

The appearance of a frequency maximum in the second-gradient model and in the Born–Karman lattice model implies that these media act as a granular filter, where only relatively low frequencies are transmitted (=low-pass filter). The wave corresponding to the maximum frequency may be interpreted as a *standing wave*, since the energy of this wave does not propagate, . The downward branch that follows after the frequency maximum relates to a negative group velocity, . As outlined in Brillouin (1946), incorporating this downward branch into a dispersion analysis introduces a non-uniqueness into the relation between the wavenumber and the frequency. For discrete systems, this non-uniqueness actually is insignificant, since waves with *a wavelength smaller than two times the particle diameter cannot be registered by a chain of discrete particles*. In other words, only wavelengths that lie inside the *first Brillouin zone k*_{x}*d*∈[−*π*, *π*], have a physical meaning. Hence, the interpretation of the dispersion curves of discrete models should be limited to the first Brillouin zone.

For the second-gradient model, a standing wave emerges at *k*_{x}≈2.05/*d*. This value has been computed by combining the condition ∂*ω*/∂*k*_{x}=0 with the dispersion relation (3.24). Likewise, for the fourth-gradient model it can be derived that a standing wave arises at *k*_{x}≈3.13/*d*. Figure 4 shows that for the fourth-gradient model the resemblance with the discrete lattice model in the first Brillouin zone is better than for the second-gradient model. This is, because the Taylor series, equation (2.6), used in the derivation of the fourth-gradient model contains terms of a higher order than the Taylor series used in the derivation of the second-gradient model. Outside the first Brillouin zone, the dispersion curve of the fourth-gradient model increases continuously, implying that higher frequencies are not filtered. Since this part of the dispersion curve is strongly affected by the truncation error of the Taylor series, equation (2.6), it should be interpreted as a purely mathematical phenomenon. Nonetheless, the monotonically increasing character of this dispersion curve confirms the unconditionally stable behaviour of the fourth-gradient model.

## 4. Higher-order continuum that includes particle rotation

The continuum models discussed previously include only the effect of particle translation. Laboratory measurements on granular media have demonstrated that at local, inhomogeneous deformation patterns, i.e. shear bands, apart from particle translations, particle rotations may also have a relevant contribution to the response (Oda *et al*. 1997). Particle rotations can be accounted for by extending equation (2.1) as(4.1)

Herein, *e*_{ijk} is the well-known permutation symbol, and is the contact rotation, which is composed by the difference in rotation of the particles *n* and *m*, see also figure 1. The rotation at the contact point *c* is energetically conjugated to the contact moment . The contact moment can be interpreted as a *point* representation of the non-uniform distribution of multiple contact forces between two arbitrary-shaped particles (Suiker *et al*. 2001*a*). Accordingly, it may account for the effect of angularity and flatness of particles. The introduction of a contact moment requires the contact law (2.2) to be extended with the relation(4.2)where the linear elastic rotational stiffness reads(4.3)

In the above expression, *Gn* is the twisting resistance, and *Gs* and *Gt* are the rolling resistances. Approximating the discrete rotations by a continuous field, the relation between the rotations of two neighbouring particles can be expressed as(4.4)where the continuous rotation in equation (4.4) may be decomposed into a spin contribution and a rigid body contribution (4.5)

The rigid body rotation is related to the anti-symmetric part of the first gradient of displacement,(4.6)or, in a reverse form,(4.7)

In the following, the lower bar denoting the continuous behaviour of a variable will be omitted. The effects of particle displacement and particle rotations are combined by truncating the Taylor series for the displacement, equation (2.6), after the third-order derivative, and the Taylor series for the rotation, equation (4.4), after the second-order derivative. Both the strain and rotational contributions are then selected up to and including the second-order derivative. Substituting equations (2.6) and (4.4) into equation (4.1) and subsequently invoking equations (2.11), (2.12), (4.5) and (4.6) leads to(4.8)which characterizes the kinematics of a so-called *second-gradient micro-polar* continuum. Apparently, equation (4.8) contains only particle displacement (derivatives) and particle spin (derivatives), which are independent degrees of freedom.

For a granular Cosserat medium, equation (4.8), reduces to (Chang & Ma 1992)(4.9)

When defining the strain variables *γ*_{ji} and *κ*_{ji} as(4.10)equation (4.9) can be recast as(4.11)

For micro-polar media, apart from the stress expression, equation (2.21), an additional stress expression needs to be derived that relates the couple stress *μ*_{ij} to local contact couples. The equilibrium condition with respect to couple stress is (Eringen 1968)(4.12)

Following the procedure proposed by Chang & Liao (1990), a polar stress *M*_{ij} is introduced(4.13)where *X*_{k} is the position vector in the *k*th direction. With equation (4.13), the equilibrium condition, equation (4.12), turns into(4.14)

Employing the above expression, the mean polar stress for a particle with volume *V*^{n} may be written as(4.15)

By invoking Gauss theorem, equation (4.15) turns into a surface integral expression(4.16)where *n*_{k} is the unit outward, normal to the particle surface *S*^{n}. Subsequently, the surface integral in equation (4.16) may be related to the contact moments at the particle surface as(4.17)where *N*^{n} is the number of contacts along the surface of particle *n*, and and are the contact force and the contact moment transmitted by neighbouring particles, respectively. Combining equations (4.16) and (4.17) leads to(4.18)

Similar to equation (2.19), the mean polar stress for a representative volume *V* of *M* particles can be computed by the volume average of the corresponding local particle stresses(4.19)

Using equation (2.20), the above expression can also be cast into a summation over the total number of particle contacts in the assembly, *N*, i.e.(4.20)

Here, use has been made of the fact that . Comparing the format of equation (4.20) with that of equation (4.13), the stress and the couple stress can be identified as(4.21)

Notice that equation (4.21*a*) is equal to the stress expression for the micro-structural model without particle rotation, equation (2.21). For notational convenience, the superimposed tilde for the stress parameters *σ*_{ij} and *μ*_{ij} will be omitted in the remainder. Combining the stress expressions (4.21), the kinematic relations (4.8) and the contact laws (2.2) and (4.2) leads, after an extended mathematical procedure, to the following constitutive form (Suiker *et al*. 2001*a*)(4.22)with the coefficients(4.23)and the package density *θ* given by equation (2.35). Note that the tangential stiffnesses *Kt* and *Gt* do not appear in equation (4.23), which is because the supposition of equal-sized spheres enables to replace *Kt* and *Gt* by *Ks* and *Gs*, respectively. It can be noticed that the Cauchy stress *σ*_{ij} in equation (4.22) relates to both symmetric strain (gradient) terms and anti-symmetric rotational (gradient) terms, whereas the couple stress *μ*_{ij} relates to only anti-symmetric rotational gradient terms. Essentially, the stiffness contributions *θKn* and *θKs* in equation (4.23) capture the homogeneous part of the material behaviour, where equations (4.23*a*) and (4.23*b*) comprise the following dependencies on the Lamé constants *λ* and *μ*(4.24)

When invoking equation (4.24), all coefficients in equation (4.23), except for *C*_{(5)} and *C*_{(6)}, can be expressed in terms of Lamé constants and the particle radius *r*. Apart from the parameters *θKn*, *θKs* and *r*, two additional independent parameters can be recognized, namely *θGn* and *θGs*. The constitutive relations (4.21) thus contain the tractable number of five independent material parameters. Experimental calibration of these parameters requires the back-fitting of inhomogeneous deformation patterns. Hereto, the applicability of measuring techniques used in the field of fracture mechanics may be examined, where inhomogeneous deformation patterns are monitored by placing a regular grid of local sample points onto the test specimen (Geers 1997).

The second-gradient micro-polar model (4.22) contains some simpler constitutive relationships, which are the linear elastic model (*C*_{(1)} to *C*_{(6)} equal zero),(4.25)the second-gradient model (*C*_{(3)} to *C*_{(6)} equal zero)(4.26)and the Cosserat model (*C*_{(1)}, *C*_{(2)} and *C*_{(4)} equal zero)(4.27)

Apparently, the Cosserat model is expressed in terms of the kinematic variables *Π*_{i}, *u*_{[i,j]} and *ϵ*_{ij}, although the use of the strain measures *γ*_{ji} and *κ*_{ji}, given by (4.10), is more common (see for example, Eringen 1968; Mühlhaus & Vardoulakis 1987; De Borst & Sluys 1991; Chang & Ma 1992; Sluys 1992). Using these strain measures, equation (4.27) becomes(4.28)

As a result of the isotropy of the underlying microstructure, the constitutive coefficients related to the strain terms *κ*_{kk} and *κ*_{ji} are identical, i.e. both are equal to 2*C*_{(5)}. Nevertheless, in general these two strain terms may be related to different constitutive coefficients.

For a medium that includes micro-polarity, the weak form of the equilibria with respect to stress and couples stress, equations (2.14) and (4.12), over a volume *V* is given by(4.29)

Integration by parts and employing Gauss theorem yields(4.30)

Invoking the strain variables, equation (4.10), the volume integrals in equation (4.30) can be rephrased as(4.31)

In order to further extend equation (4.31), first the constitutive form, equation (4.22), needs to be rephrased in terms of the strain measures given by equation (4.10), i.e.(4.32)

Inserting the above expression into equation (4.31), followed by integrating by parts, and applying Gauss theorem, yields(4.33)

Combining equation (4.33) with equation (4.30) leads to(4.34)with the tractions *t*_{j}, *m*_{j} and *q*_{ij} at the boundary surface as(4.35)where *σ*_{ij} and *μ*_{ij} are presented by equation (4.32). The natural boundary data given by the tractions *t*_{j}, *m*_{j} and *q*_{ij} is complementary to essential boundary data for *u*_{j}, *ω*_{j} and *γ*_{ij}, respectively. As mentioned in §2*c*, independent boundary conditions are related to the kinematic derivatives in the direction normal to the boundary. A procedure to retrieve the independent boundary conditions for a three-dimensional boundary value problem can be found in Mindlin (1965).

The left-hand side of equation (4.34) can be recognized as the variation of the potential energy, . Accordingly, the potential energy density for the second-gradient micro-polar model may be expressed as(4.36)

Using the knowledge obtained from the energy considerations in §3*b*, it can be concluded that in equation (4.36) the terms with a minus sign, i.e. the coefficients *C*_{(1)} and *C*_{(2)} have a *destabilizing* effect on the structural response. Thus, the Cosserat model, equation (4.27) (or equation (4.28)), which does not contain strain terms related to *C*_{(1)} and *C*_{(2)}, is unconditionally stable if the microstructural material parameters, equation (4.23) are used with realistic choices for the contact stiffnesses, i.e. *Kn*≥*Ks*≥0, *Gn*≥*Gs*≥0. This is due to the fact that all kinematic descriptors in this model (i.e. strain, first gradient of rotation) are energetically conjugated to corresponding stresses (i.e. stress, couple stress). In order to avoid the presence of destabilizing terms in the second-gradient micro-polar formulation, the second-gradient strain terms in equation (4.22) should have been energetically coupled to additional, conjugated stress terms, comparable to the concept of ‘double stress’ introduced in the pioneering work of Mindlin (1964). However, the avoidance of destabilizing terms can have an unfavourable effect on the agreement between the responses of a kinematically enhanced continuum model and corresponding discrete lattice models. More specifically, enhanced continuum models with only stabilizing terms often reveal P-wave and S-wave dispersion curves which, in the long wavelength domain of the *ω–k* plane, lie *above* the P-wave and S-wave dispersion curves for a classic elastic continuum model (see for instance, Fig. 3 in Mindlin 1964). This may cause them to be in a very limited agreement with the P-wave and S-wave dispersion curves for corresponding discrete lattices, which generally lie *below* the P-wave and S-wave dispersion curves for a classic elastic continuum. In other words, when avoiding destabilizing terms in a micro-mechanically based enhanced continuum model, the P-wave and S-wave phase velocity and group velocity related to the underlying microstructure may be considerably overpredicted.

## 5. Strain-gradient micro-polar continuum versus hexagonal lattice

In this section, the wave propagation characteristics of the second-gradient micro-polar model (4.22) and the Cosserat model (4.27) will be compared to the wave propagation characteristics of a two-dimensional, 7-cell hexagonal lattice, depicted in figure 5. The purpose of this comparison is to exhibit up to which deformation level the enhanced continuum models accurately represent the inhomogeneous effects caused by particle translation and particle rotation. Figure 5 shows that the hexagonal lattice consists of a cell with dimensionless coordinates (*m*, *n*), which is connected to six neighbouring cells at a distance *d*, by means of longitudinal springs *Kn*, shear springs *Ks* and rotational springs *G*. Each cell in the lattice has three degrees of freedom, namely two translations and and a rotation . The lattice characteristics are thus equivalent to the micro-structural characteristics of the second-gradient micro-polar model.

### (a) Dispersion relations for hexagonal lattice

The equations of motion of the hexagonal lattice can be derived in a straightforward manner by taking the Lagrangian *L* of the lattice cell (*m*,*n*) as a point of departure(5.1)where the potential energy follows from(5.2)with Δ*l*_{(i)}, Δ*s*_{(i)} and Δ*r*_{(i)} the variations in length of the *i*th spring in the longitudinal direction, the shear direction and the rotational direction, respectively. Expressions for these length variations can be found in Appendix A. The kinetic energy for the cell (*m*, *n*) is formulated by(5.3)with *M* the mass and *I* the moment of inertia of the cell. The dot on top of the degrees of freedom denotes the total derivative with respect to time. Combining equations (5.1)–(5.3), and applying the Lagrange equations (Landau & Lifshitz 1976)(5.4)provides the equations of motion for the inner cell (*m*, *n*)(5.5)

In order to analyse the body propagation through the lattice, a solution to equation (5.5) will be sought in the form of plane, harmonic waves,(5.6)with , and the wave amplitudes, *ω* the angular frequency, and *k*_{x} and *k*_{z} the wavenumbers in the *x*- and *z*-directions, respectively. Substituting (5.6) into (5.5) yields a set of three homogeneous algebraic equations(5.7)which reflects the dispersion curves for the compression wave, the shear wave and the micro-rotation wave. Before computing these dispersion curves, first a similar set of equations will be derived for the second-gradient micro-polar continuum and the Cosserat continuum.

### (b) Dispersion relations for second-gradient micro-polar continuum and Cosserat continuum

The derivation of the equations of motion for continua that incorporate micro-polar effects starts with the formulation of the balance of linear momentum and the balance of angular momentum. For two-dimensional micro-polar continua that include the Cauchy stress *σ*_{ij} and the couple stress *μ*_{ij}, these equations have the form (Eringen 1968)(5.8)where, in correspondence with the lattice model, a restriction has been made to the *x**–z* coordinate system. Furthermore, *J* is the moment of inertia per unit volume. Inserting the constitutive relations for the second-gradient micro-polar model, equation (4.22), into equation (5.8), and invoking the kinematic relations, equations (2.11), (2.12), (4.5) and (4.6), leads to(5.9)where . In order to satisfy the equations of motion, equation (5.9), the following harmonic expressions are substituted(5.10)which yields(5.11)

The above set of algebraic equations characterizes the propagation of the compression wave, the shear wave and the micro-rotation wave in the second-gradient micro-polar continuum. For the Cosserat model a similar set of equations can be obtained, by putting *C*_{(1)}, *C*_{(2)} and *C*_{(4)} to zero, which renders(5.12)

The propagation of free body waves in the Cosserat continuum, as captured by equation (5.12), has been treated previously by other authors (Eringen 1968; De Borst & Sluys 1991). Nevertheless, in their contributions the Cosserat model has been expressed in macroscopic stiffness parameters and a phenomenological length scale parameter *l*, while here the material parameters have been related to micro-scale material properties.

### (c) Discussion of dispersion curves

It is required that in the long-wave limit *k*→0, the dispersion relations for the lattice model, equation (5.7), and the enhanced continuum models, equations (5.11) and (5.12), have equivalent characteristics. This condition can be satisfied by considering the lattice dispersion relations, equation (5.7), taking a second-order Taylor approximation (=long-wave approximation) with respect to the normalized wavenumbers *k*_{x}*d* and *k*_{z}*d*, and projecting this approximation onto the Cosserat dispersion relations (5.12). For a more thorough discussion on this homogenization procedure, see Kunin (1982). Accordingly, the following equalities between the macro-level continuum parameters and the micro-level lattice parameters are found(5.13)

From these equalities, equation (5.13*e*) is actually redundant, since the lattice shear stiffness *Ks* is also determined by the Lamé constants *λ* and *μ*. Note that the micro-macro relations (5.13*c–f*) are different from those presented by equations (4.23*a*,*b*,*e*,*h*), respectively. This is due the fact that a different homogenization procedure has been used, as well as to the fact that the hexagonal lattice reflects a two-dimensional, regular particle configuration, while the second-gradient micro-polar model, equation (4.22), has been derived from a three-dimensional, random particle configuration.

The dispersion curves for the hexagonal lattice and the enhanced continuum models can be determined by numerically solving the corresponding set of algebraic equations, equations (5.7), (5.11) and (5.12). A non-trivial solution can be found if, for a given set of elastic parameters, the determinant Δ(*ω*, *k*_{x}, *k*_{z}) of the specific system of equations is equal to zero. The parameters used for the computation are: *λ*=55.5 MPa, *μ*=83.3 MPa, *d*=2*r*=50 mm, *C*_{(6)}=69.5 kN, *ρ*=1800 kg m^{−3} and *J*=0.5625 kg m^{−1}. The value for *C*_{(6)} relates to a micro-level stiffness ratio *G*/*Ks*=*d*^{2}, whereas the value for the moment of inertia per unit volume *J* reflects a medium of equal-sized circular particles, for which the moment of inertia of an individual particle ensues from *I*=*Md*^{2}/8.

Figure 6 depicts the dispersion curves of the compression wave, the shear wave and the micro-rotation wave. The direction of propagation of the waves corresponds to . For the discrete lattice, this causes the first Brillouin zone to be represented by *kd*∈[−*π*, *π*]. At the end of the first Brillouin zone, all three waves in the discrete lattice attain a frequency maximum, where the group velocity of the waves is zero, *c*^{g}=0. As mentioned before, these waves may be characterized as standing waves. For the second-gradient micro-polar model, the compression wave and the shear wave also reveal a frequency maximum, though these maxima and the corresponding wavenumbers are smaller than those for the lattice model. The micro-rotation wave in this model does not have a frequency maximum, because the polynomial, equation (4.4), describing the contact rotation only contains terms up to and including the second order. If this Taylor series had been truncated after the fourth-order term, the micro-rotation wave most likely would also have exhibited a frequency maximum. For the Cosserat model, the Taylor series regarding the contact displacement (2.6) and the contact rotation (4.4) are both truncated after the first order term, which is the reason why the dispersion curves of this model do not exhibit frequency maxima at all.

Up to a wavenumber *k*≈1.5/*d*, which relates to a wavelength *Λ*≈4*d*, the overall agreement between the second-gradient micro-polar continuum and the lattice is good, in a sense that for the dispersion curves the relative difference in frequency is less than 10%. For the Cosserat model, the accuracy is somewhat less, due to a lesser extent of the Taylor series (2.6) and (4.4).

In §4, it has been discussed that the constitutive terms *C*_{(1)}, *C*_{(2)} and *C*_{(4)} have a destabilizing effect on the response of the second-gradient micro-polar model, equation (4.22). Although the presence of destabilizing terms not necessarily induces overall instability, see §3*b*, the second-gradient micro-polar continuum nevertheless becomes unstable for *k*>2.6/*d*, i.e. when the wavelength of the deformation pattern becomes smaller than 2.4 times the particle diameter. At this stage the dispersion curve of the compression wave intersects with the frequency axis *ω*=0, see figure 6, after which the frequency becomes imaginary and the response increases exponentially with time. When using the second-gradient micro-polar model in the analysis of dynamic boundary value problems, wavelengths smaller than the critical wavelength should be excluded from the response to avoid the appearance of this physically unrealistic instability.

Instead of using a 7-cell hexagonal lattice to reveal the accuracy level of higher-order continua, other lattice configurations could have been chosen as well, such as the 9-cell square lattice described in Suiker *et al*. (2001*c*). The different topology of the 9-cell square lattice causes the micro–macro relationships to be different from equations (5.13). Moreover, the anisotropy and inhomogeneity characteristics resulting from the discrete nature of the lattice are different, although this becomes manifest only for very short deformation patterns, i.e. when the wavelength is less than 6 times the particle diameter *d* (Suiker *et al*. 2001*c*).

## 6. Discussion

The derivation of continuum formulations for isotropic granular materials consisting of equal-sized spherical particles has been discussed. By employing a homogenization method known as the ‘micro-structural approach’, several kinematically enhanced continuum models have been derived. The constitutive coefficients in these enhanced continua have been expressed in terms of micro-structural properties, such as particle radius, contact stiffness and package density. The accuracy level at which the discrete structure is simulated by the continuum models depends on the truncation level of the Taylor series used for describing the particle kinematics. This has been illustrated by comparing the dispersion curves of continuum models with those of two discrete lattices, which are the one-dimensional Born–Karman lattice and the two-dimensional 7-cell hexagonal lattice.

As exemplified, including only the effect of particle displacement leads to a strain-gradient model, while adding the effect of particle rotation leads to a strain-gradient micro-polar model. By analysing the dispersion curves for body wave propagation, it has been shown that enhanced continua and discrete lattices may retain waves with frequencies lying above a specific threshold frequency. This filtering effect is the result of including the particle size into the model formulation. The maximum frequency transmitted thereby relates to a minimum transmittable wavelength, which for a discrete lattice equals two times the particle diameter.

The stability of the continuum models has also been examined. It has been shown that some enhanced continuum models may reveal a physically spurious instability when the wavelength of the imposed deformation pattern becomes smaller than a specific *critical wavelength*. Nonetheless, this occurs at deformation levels where the Taylor approximation reflecting the particle kinematics has already lost its accuracy. The critical wavelength at which the instability emerges can be found either from examination of the dispersion relations for body wave propagation, or from examination of the second-order work. Because of the physically irrational character of the instability, the use of these models in dynamic boundary value problems requires the analyst to exclude deformation components from the response with a wavelength smaller than the critical wavelength. Only then, a stable response can be computed. Discrete lattices do not suffer from such a deficiency, since these models describe the particle kinematics in an ‘exact manner’. In fact, for discrete models wavelengths smaller than two times the particle diameter will be naturally excluded from the solution to a boundary value problem, as they cannot be uniquely registered by the particles.

## Acknowledgments

The authors highly appreciate the discussions on granular mechanics and wave propagation with Prof. Ching S. Chang of the University of Massachusetts, Amherst, United States and Dr Andrei V. Metrikine of the Russian Academy of Sciences, Nizhny Novgorod, Russia. Their knowledge and experience in these research fields has proven to be extremely valuable for the preparation of this article.

## Footnotes

One contribution of 7 to a Theme Issue ‘Thermodynamics in solids mechanics’.

↵Apart from explicitly posing the static hypothesis, it also can be derived by postulating that the mean displacement field in a granular assembly is the

*best-fit*for the actual displacement field (Liao*et al*. 1997, 2000).↵Prescribing the parameters

*u*_{x},*u*_{x,xx}and*u*_{x,xxxx}at the boundaries implies that the variation of these parameters at the boundaries is equal to zero.↵This requirement in fact implies that the

*random*granular structure in the continuum model and the*regular*granular structure in the lattice model have an equivalent package density and macroscopic stiffness.- © 2005 The Royal Society