# Model for granular materials with surface energy forces

Ching Chang, Pierre-Yves Hicher

► To cite this version:

Ching Chang, Pierre-Yves Hicher. Model for granular materials with surface energy forces. Journal of Aerospace Engineering, 2009, 22 (1), pp.43-52. 10.1061/(ASCE)0893-1321(2009)22:1(43) . hal-01007054

HAL Id: hal-01007054

https://hal.science/hal-01007054

Submitted on 15 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Model for Granular Materials with Surface Energy Forces

C. S. Chang<fnref n="1" /> and P.-Y. Hicher<fnref n="2" />

In light of environmental differences (such as gravitational fields, surface temperatures, atmospheric pressures, etc.), the mechanical behavior of the subsurface soil on the Moon is expected to be different from that on the Earth. Before any construction on the Moon can be envisaged, a proper understanding of soil properties and its mechanical behavior in these different environmental conditions is essential. This paper investigates the possible effect of surface-energy forces on the shear strength of lunar soil. All materials, with or without a net surface charge, exhibit surface-energy forces, which act at a very short range. Although, these forces are negligible for usual sand or silty sand on Earth, they may be important for surface activated particles under extremely low lunar atmospheric pressure. This paper describes a constitutive modeling method for granular material considering particle level interactions. Comparisons of numerical simulations and experimental results on Hostun sand show that the model can accurately reproduce the overall mechanical behavior of soils under terrestrial conditions. The model is then extended to include surface-energy forces between particles in order to describe the possible behavior of lunar soil under extremely low atmospheric pressure conditions. Under these conditions, the model shows that soil has an increase of shear strength due to the effect of surface-energy forces. The magnitude of increased shear strength is in reasonable agreement with the observations of lunar soil made on the Moon's surface.

Granular material; Stress-strain relations; Anisotropy; Elastoplasticity; Moon.

## Introduction

The soil on the Moon formed by space weathering processes and dynamic impacts of micrometeorites is composed of aggregates of minerals, rocklets, and glasses, welded together into agglutinates. The glass, which is abundant in lunar soil, is present in the form of sharp, abrasive, interlocking, fragile shards and fragments. The lunar subsurface soil can be classified as dry silty sand (Perko et al. 2001). However, observations and measurements conducted during Surveyor (1966–1968), Apollo (1968–1972), and Luna (1959–1976) missions indicate that lunar soil has an unusually high cohesion in comparison to the same type of soils tested under the Earth’s conditions. According to Buzz Aldrin, astronaut on Apollo 11, the soil on the Moon “has the cohesive property that wet sand would have” (Costes et al. 1970). A stable vertical trench excavated during the Surveyor 7 mission is shown in Fig. 1(a) and a footprint left by Surveyor 3 demonstrates the fineness and cohesiveness of the lunar soil [Fig. 1(b)].

The fine particles of lunar soil are the products of the continual impact on the surface by meteoroids, which smash and grind rocks into soil and weld soil into new rocks. The writers we are interested in the behavior of soils in a typical regolith of the Moon—a surface layer that consists of loose sand and rock fragments, which overlies the solid rock. The range of grain size distribution is summarized in a soil gradation curve given in Fig.

2 (Costes et al. 1970; Mitchell 1972; Scott 1973) based on the core samples obtained from the surface layer of the Moon (Apollo 11 and 12). According to the particle size distribution, the lunar soil can be classified as *silty sand* following the Unified Soil Classification System, ASTM standards.

The writers’ knowledge on engineering properties of lunar soil are derived mainly from operation performed during the Surveyor and Appolo programs, including bulk and core samples returned to the Earth, trench excavation, ground penetration, and some simple physical experiments. In laboratory studies, it was found that soils tested under a chamber with ultrahigh vacuum and high temperature had an increase up to 13° in angle of internal friction and 1.1 kPa in cohesion (Bromwell 1966; Nelson 1967; Perko et al. 2001). This increase of soil strength is consistent with the observed stable vertical trend excavation in dry silty sand on the Moon.

For soil on the Earth, cohesion is originated from electrostatic forces and surface-energy forces, which includes van der Waals forces (intermolecular potential energy). The surface-energy forces act at very short range. Comparing to the gravitational forces, the surface-energy forces are negligible on the Earth for soils with particle sizes larger than around 0.06 mm. Thus, terrestrially, the forces are negligible for sand or silty sand. However, as particle size decreases to the size of clay, the surface-energy forces have a greater influence on its strength. The effect of surface-energy forces on friction and adhesion properties of small particles can be found in the literature of powder technology, particulate technology, and interface sciences (Derjaguin 1934; Derjaguin et al. 1975; Johnson 1971; Molerus 1978; Thorton and Ning 1998; Tomas 2001).

Due to various environmental differences on the Moon (such as the gravitational fields, temperature extremes, low to practically nonexistent atmospheric pressure, plus the space weathering, which includes surface activating effects like UV radiation, the solar wind and galactic particle fluxes), it seems plausible to

<footnote n="1">¹Dept. of Civil and Environmental Engineering, Univ. of Massachusetts, Amherst, MA 01002. E-mail: chang@ecs.umass.edu</footnote>

<footnote n="2">²Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes–Univ. of Nantes, France.</footnote>
<figure><img src="image_1.png" /><figcaption>(a) Trench excavated by Surveyor 7 Spacecraft in Lunar subsurface soil (adapted from Perko et al., 2001).</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>(b) The footprint of Spacecraft Surveyor 3–Apollo 12 (adapted from NASA Image ID number: AS12-H-48-7110).</figcaption></figure>

**Fig. 1.** (a) Trench excavated by Surveyor 7 Spacecraft in Lunar subsurface soil (Perko et al. 2001, ASCE); (b) the footprint of Spacecraft (Surveyor 3–Apollo 12) (NASA Image ID number AS12-H-48-7110)

consider surface-energy forces as a dominant cause for altering the mechanical behavior of lunar soil with the size of silty sand. An issue that may be addressed is to what order of magnitudes the surface-energy forces can influence the strength and deformation of lunar soils.

In this paper, a model is developed that accounts for the surface-energy forces between particles and is capable of modeling its effect on the shear strength of the soil assembly. For this purpose, a microstructural modeling approach is adopted (Chang and Hicher 2005; Hicher and Chang 2005). Models using this approach can also be found in the work by Jenkins (1988), Walton (1987), Rothenburg and Selvadurai (1981), Chang (1988), Emericault and Cambou (1996), Liao et al. (2000), Kruyt and Rotherburg (2002), among others. Elastoplastic models using this approach can be found in Chang et al. (1992a,b), Suiker and Chang (2004).

It is noted that the above-mentioned microstructural continuum approach (Chang and Hicher 2005) is not a “complete” micromechanical approach. As the complete microstructural details for a particle assembly are not known, the complex deformation behavior of the particle assembly cannot be derived directly from the first-principles theory based solely on particle-level properties. Two major simplified assumptions have been applied.

1. Besides the parameters of particle-level properties, the present model adopts some macroscale parameters, which are related to the degree of interlocking and porosity for the soil assembly. These macroscale parameters are adopted because interparticle behavior is not exclusively a local phenomenon; it is considerably influenced by the degree of interlocking and porosity for a soil assembly. In the present model, these parameters are based on the critical void ratio concept, and are empirically determined from testing on soil specimens.

2. Particle rotation is neglected in the model derivation. Particle rotation is a significant factor that influences the stress–strain

behavior of assembly. A model considering particle rotation would result in the use of higher-order stresses and strain such as micropolar strain and couple/Cosserat stress (Chang and Liao 1990; Chang and Ma 1992; Chang and Gao 1996, Luding et al. 2001; Chang and Kuhn 2005; Tejchman and Bauer 2005). Currently, there is no stress–strain model with explicit consideration of particle rotation that can successfully simulate complex behavior of soil. The difficulties of establishing such a continuum model mainly arise from the lack of experimental results on measured particle rotation field and related couple stress for an assembly of soil particles.

With the previous discussion on restrictions and simplified assumptions, the present model should be regarded as a model of semiempirical nature. It is necessary to calibrate parameters from soil tests in order to predict complex behavior with good quantitative agreement. The main advantage of the present model over conventional continuum constitutive models is its microscale consideration, which allows one to conveniently extend the model for new phenomenon at the particle level, such as the surface-energy forces.

However, adding the surface-energy term from first-principle physics into an otherwise semiempirical model may elicit the following concern: how can one part of the physics be treated from first principles if the other parts of the physics are treated only empirically? The answer is that even though the surface-energy forces are added from first principle physics, it does not alter the empirical nature of the model. The model is aimed to predict the “order of magnitude” of the surface-energy effect.

In what follows, we first describe the formulation of this model, which takes account of the surface-energy forces between particles. The model performance is then demonstrated through the results of a simulation of fine Hostun sand tested under terrestrial conditions (i.e., without surface-energy forces). The fine sand has a mean particle diameter of 0.5 mm, loaded to failure in triaxial testing cells under different confining stresses. Finally, the effect of surface-energy forces between particles is introduced to predict the stress–strain–strength behavior of lunar soil.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Approximated gradation range for Apollo 11 and 12 soil core samples</figcaption></figure>
## Stress-Strain Model

In this model, we envision a granular material as a collection of particles. The deformation of a representative volume of the material is generated by the mobilization of contact particles in all contacts. Thus, the stress-strain relationship can be derived as an average of the mobilization behavior of local contact planes in all contacts. For a contact plane in the αth orientation, the local forces $f_i^α$ and the local movements $δ_i^α$ can be denoted as follows: $f_i^α = \{f_n^α, f_s^α, f_t^α\}$ and $δ_i^α = \{δ_n^α, δ_s^α, δ_t^α\}$, where the subscripts n, s, and t represent the components in the three directions of the local coordinate system. The direction normal to the plane is denoted as n; the other two orthogonal directions, s and t, are tangential to the plane. Rotation of particles is not considered here.

The forces and movements at contact planes of all contacts are suitably superimposed to obtain the macroscopic stress strain tensors. The macroscopic stiffness tensor is obtained on the condition that the rate of energy dissipation expressed in terms of the macrostress and strain must be equivalent to that expressed in terms of microforces and movements. In such formulation, it has usually been assumed that the microstructure is statically constrained, which means that the forces on each contact plane are assumed equal to the resolved components of the macroscopic stress tensor. Another equally simple possibility is to assume a kinematically constrained microstructure, in which the movements, rather than the forces on a contact plane, are the resolved components of the macroscopic strain tensor.

The kinematically constrained models are more popularly used in concrete models. The overriding reason for using this approach is that, in case of strain softening, it is easier to construct a stable model using a kinematic rather than a static constraint. However, the kinematic constraint gives a more restrictive deformation pattern, thus providing a stiffer result, especially when granular materials are subjected to high deviatoric stresses (Chang and Misra 1990). Further discussion on this issue can be found in the work by Chang and Gao (1996), Kruyt and Rothenburg (2002), and Kruyt (2003).

## Interparticle Behavior

### van der Waals Forces

The surface-energy forces pull soil particles together, thus increase the shear strength of soil. The physics of the surface-energy force between two particles is reviewed in this section. The adhesive forces between two solid grains mainly result from electrostatic force and van der Waals forces. The electrostatic component in lunar soil is assumed to be negligible (Perko et al. 2001). Therefore, in this paper, the surface-energy forces are calculated from van der Waals energy fields. There are several situations that can alter the surface properties and surface energy, such as monolayers of water or passivation layers. For simplicity, lunar soil grains are represented by spheres of equal radii that are separated by thin layers of adsorbed air molecules.

van der Waals forces between two bodies are derived from the dispersion interaction energy between two identical atoms or molecules (Israelachvili 1992):

$$W(r) = -\frac{C}{r^6} \qquad (1)$$

where r=distance between the two atoms and C=London dispersion coefficient.

Then, with the assumption of additivity of these attractive forces, the interaction energy of a molecule located at a distance D of the planar surface of a solid made up of like molecules is the sum of its interaction with all the molecules in the body:

$$W(D) = -\frac{\pi C \rho}{6D^3} \qquad (2)$$

Where ρ=number of atoms per unit volume. The corresponding van der Waals force is

$$F = \frac{\delta W(D)}{\delta D} = \frac{\pi C \rho}{2D^4} \qquad (3)$$

It is then possible to calculate the interaction energy between two solids as, e.g., two spheres:

$$W = -\frac{AR_1R_2}{6D(R_1 + R_2)} \qquad (4)$$

D=distance between two spheres; $R_1$ and $R_2$=radii of the two spheres; $A=\pi^2 C \rho_1 \rho_2$=Hamaker constant; and $ρ_1$ and $ρ_2$=number of atoms per unit volume for the two bodies. For two flat surfaces, one obtains the interaction energy per unit area:

$$W = -\frac{A}{12\pi D^2} \qquad (5)$$

The van der Waals force between two solids in contact can then be computed from the interaction energy. One obtains the following two cases;

Between two spheres:

$$f = 2\pi W(D) \frac{R_1 R_2}{R_1 + R_2} \qquad (6)$$

Eq. (6) is called the Derjaguin approximation (Israelachvili 1992). In the case of two identical spheres, Eq. (6) becomes

$$f = \frac{AR}{12D^2} \qquad (7)$$

Between two flat surfaces, per unit area:

$$f = \frac{A}{6\pi D^3} \qquad (8)$$

Let us now consider two soil grains as two identical elastic spheres. If subjected to an external force f, the two spheres will deform and create a flat circular contact area with a radius a

$$a = \left( \frac{3(1 - v_p^2)R}{4E_p} \right)^{1/3} f^{1/3} \qquad (9)$$

where $E_p$ and $v_p$=Young's modulus and Poisson's ratio of the particles, respectively. The derivation of Eq. (8) can be found in Derjaguin et al. (1975), Johnson (1971), Dahneke (1972), and Valverde et al. (2001). It is noted that for simplicity, we consider only the case of elastic contact flattening [Eq. (9)]. For the case of elastic-plastic flattening, refer to Molerus (1978), Thornton and Ning (1998), and Tomas (2001).

Under these simplified conditions, van der Waals forces acting on the two particles can be considered as the sum of two terms, one due to the interaction between two flat surfaces of area S =πa² and one along the remaining surface of the two spheres. Using the Derjaguin approximation for the second term, we obtain the expression of the van der Waals force between two particles:
$$f = \frac{A}{6D^3}a^2 + \frac{AR}{12D^2} \qquad (10)$$

where A=Hamaker coefficient; D=thickness of molecules layer between two particles; and R=radius of particles. Hamaker constant A was estimated to be $4.3 \times 10^{-20}$ J for lunar soil and $1.5 \times 10^{-20}$ J for terrestrial quartz sand (Perko et al. 2001).

The thickness of molecules layer between two particles D are highly dependent on the atmospheric pressure and composition. On the Moon, the atmospheric pressure is nearly zero, which can lead to a very thin layer of molecules between two particles compared to that under a terrestrial environment. Therefore, according to Eq. (10), it is reasonable to expect that the surface-energy forces between particles are much higher than those between particles under a terrestrial environment.

It is noted that if the radius a of the contact area in Eq. (10) increases with confining stress of a specimen. Thus, the surface-energy force also increases with confining stress, which indicates that the surface-energy force will contribute to the shear strength not only on the cohesive component but also on the frictional component.

The orientation of a contact plane between two particles is defined by the vector perpendicular to this plane. On each contact plane, an auxiliary local coordinate can be established as shown in Fig. 3. The surface-energy force f given in Eq. (10) represents only the force magnitude. Its direction is normal to the interparticle contact plane given by $n_i^\alpha$. Thus in a vector form, the surface-energy force $f_i^{\alpha(\text{SE})} = fn_i^\alpha$.

Between two particles, the interparticle force $f_i^\alpha$ can be decomposed into two components: (1) due to applied load on the boundary of the soil assembly $f_i^{\alpha(A)}$, and (2) due to surface-energy forces between particles $f_i^{\alpha(\text{SE})}$. Thus,

$$f_i^\alpha = f_i^{\alpha(A)} + f_i^{\alpha(\text{SE})} \qquad (11)$$

The surface-energy force is a function of contact area as shown in Eq. (10). However, the contact area is, in turn, a function of interparticle force $f_i^\alpha$ [see Eq. (9)], thus the expression of Eq. (11) is implicit and nonlinear in nature.

## Interparticle Force-Displacement Relationship

### Elastic Part

The contact stiffness of a contact plane includes normal stiffness, $k_n^\alpha$, and shear stiffness, $k_r^\alpha$ (assuming $k_r^\alpha = k_t^\alpha = k_s^\alpha$). The elastic stiffness tensor is defined by

$$f_i^\alpha = k_{ij}^{\alpha e} \delta_j^{\alpha e} \qquad (12)$$

which can be related to the contact normal and shear stiffness

$$k_{ij}^{\alpha e} = k_n^\alpha n_i^\alpha + k_r^\alpha (s_i^\alpha s_j^\alpha + t_i^\alpha t_j^\alpha) \qquad (13)$$

where n, s, and t=three orthogonal unit vectors that form the local coordinate system (see Fig. 3). The vector n is outward normal to the contact plane. Vectors s and t are on the contact plane.

The value of the stiffness for two elastic spheres can be estimated from Hertz-Mindlin's formulation (Mindlin 1969). For sand grains, a revised form was adopted (Chang et al., 1989), given by

<figure><img src="image_4.png" /><figcaption>Fig. 3. Local coordinate at interparticle contact</figcaption></figure>

$$k_n = k_{n0} \left( \frac{f_n}{G_g l^2} \right)^n, \quad k_r = k_{r0} \left( \frac{f_n}{G_g l^2} \right)^n \qquad (14)$$

where $G_g$=elastic modulus for the grains; $f_n$=contact force in normal direction; $l$=branch length between the two particles; $k_{n0}$, $k_{r0}$, and $n$=material constants. For two spherical particles, the branch length is the same as the particle size $l=d$. If the Hertz-Mindlin's contact formulation is used, the value of $k_{n0}$ in Eq. (14) can be expressed as

$$k_{n0} = G_g \frac{d}{2} \left( \frac{\sqrt{12}}{1 - v_g} \right)^{2/3} \qquad (15)$$

### Plastic Part

The elastic sliding behavior between two particles does not have a coupling effect (i.e., the sliding direction is on the contact plane; no movement component is in the normal direction.). However, the plastic sliding often occurs with an upward or downward movement, thus the shear dilation/contraction takes place. The stress dilatancy is a well-recognized phenomenon in sand [see discussions in the work by Taylor (1948), Rowe (1962), and Goddard and Bashir (1990)], and should be correctly modeled. The dilatancy effect has been described by (Chang and Hicher 2005)

$$\frac{\dot{\delta}_n^p}{\dot{\Delta}^p} = \frac{f_r}{f_n} - \tan \phi_0 \qquad (16)$$

where $\phi_0$=material constant which, in most cases, can be considered equal to the interparticle friction angle $\phi_\mu$. Note that the shear force T and the rate of plastic sliding $\dot{\Delta}^p$ are defined as

$$f_r = \sqrt{f_s^2 + f_t^2}, \quad \dot{\Delta}^p = \sqrt{(\dot{\delta}_s^p)^2 + (\dot{\delta}_t^p)^2} \qquad (17)$$

The yield function is assumed to be of Mohr-Coulomb type,

$$F(f_i, \kappa) = f_r - f_n \kappa (\Delta^p) = 0 \qquad (18)$$

where $\kappa(\Delta^p)$=isotropic hardening/softening parameter. When $F > 0$, it indicates loading, otherwise unloading. The hardening function, defined as a hyperbolic relationship between $\kappa$ and $\Delta^p$, involves two material constants: $\phi_p$ and $k_{p0}$ (Chang and Hicher 2005).

$$\kappa = \frac{k_{p0} \tan \phi_p \Delta^p}{|f_n| \tan \phi_p + k_{p0} \Delta^p} \qquad (19)$$

The value of $\kappa$ asymptotically approaches $\tan \phi_p$. The initial slope of the hyperbolic curve is $k_{p0}$. On the contact plane, under a yield condition, the direction of plastic shear sliding $\dot{\Delta}^p$ follows the associated flow rule, thus, it is perpendicular to the yield surface. However, the plastic movement in the direction normal to
the contact plane is governed by the stress-dilatancy equation shown in Eq. (16). Thus the overall flow rule is nonassociated.

## Interlocking Influence

One of the important elements to be adopted in granular modeling is the critical state concept. Under critical state, the granular material remains at a constant volume whereas it is subjected to a continuous distortion. The void ratio corresponding to this state is denoted as $e_c$.

The critical void ratio $e_c$ is a function of the mean stress. The relationship has traditionally been written as follows:

$$e_c = \Gamma - \lambda \log(p') \quad \text{or} \quad e_c = e_{\text{ref}} - \lambda \log\left(\frac{p'}{p_{\text{ref}}}\right) \quad (20)$$

where $\Gamma$ and $\lambda$=two material constants; $p'=\text{mean stress of the packing, and } (e_{\text{ref}}, p_{\text{ref}})=\text{reference point on the critical state line.}$

The interparticle friction angle $\phi_\mu$ is a constant for the material. However, the peak friction angle, $\phi_p$, on a contact plane between two particles is dependent on the degree of interlocking of neighboring particles, which can be related to the state of packing void ratio $e$ by

$$\tan \phi_p = \left(\frac{e_c}{e}\right)^m \tan \phi_\mu \quad (21)$$

where $m$=material constant (Biarez and Hicher 1994).

For dense packing, the peak frictional angle $\phi_p$ is greater than $\phi_\mu$. When the packing structure dilates, the degree of interlocking and the peak frictional angle are reduced, which results in a strain-softening phenomenon.

## Elastoplastic Force-Displacement Relationship

With the elements discussed previously, the relations between the rate of force displacement for two particles can be derived that includes both elastic and plastic behavior, given by

$$\dot{f}_i^\alpha = k_{ij}^{\alpha p} \delta_j^\alpha \quad (22)$$

The detailed expression of $k_{ij}^{\alpha p}$ can be found in Chang and Hicher (2005).

## Stress-Strain Relationship

### Macro- and Microrelationships

The stress-strain relationship for an assembly can be determined from integrating the behavior of all interparticle contacts. During the integration process, a relationship is required to link the macro- and microvariables. Using the static hypotheses proposed by Liao et al. (2000), we obtain the relation between the macrostrain rate and interparticle displacement rate (here, we do not consider the finite-strain condition)

$$\dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \delta_{j}^{\alpha} l_{k}^{\alpha} \quad (23)$$

where $\dot{\delta}_j$=relative displacement rate between two contact particles and the branch vector $l_k$=vector joining the centers of two contact particles. It is noted that contact particles include both direct contact and indirect contact of neighboring particles associate with a Voronoi polyhedron as discussed by Cambou et al. (2000). For convenience, we let $N$ be the total number of contacts.

The variables $\dot{\delta}_j^\alpha$ and $l_k^\alpha$ are defined, respectively, as the values of $\dot{\delta}_j$ and $l_k$ associated with the $\alpha$th contact. The fabric tensor $A_{ik}$ in Eq. (23) is defined as

$$A_{ik} = \sum_{\alpha=1}^{N} l_i^{\alpha} l_k^{\alpha} \quad (24)$$

Using the principle of energy balance and using Eq. (23), the mean force rate on the $\alpha$th contact is

$$\dot{f}_j^\alpha = \dot{\sigma}_{ij} A_{ik}^{-1} l_k^\alpha V \quad (25)$$

In Eq. (25), the stress increment $\dot{\sigma}_{ij}$ can be obtained by the contact forces and branch vectors for all contacts (Christofferson et al. 1981; Rothenburg and Selvadurai 1981):

$$\dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha} l_{i}^{\alpha} \quad (26)$$

Applying the defined contact force in Eq. (25), Eq. (26) is unconditionally satisfied. Because of its approximation nature, Eq. (25) can be viewed as an averaged solution, in which the interparticle force can be regarded as the mean value for forces on all contact planes of the same orientation. For convenience, we also regard the branch length as the mean value for all contact planes of the same orientation.

The interparticle force in of Eq. (26) consists of two parts: $f_i^{\alpha(A)}$ due to applied load on the boundary of the soil assembly, and $f_i^{\alpha(SE)}$ due to surface-energy forces between particles [see Eq. (11)]. Thus the stress can be separated into two parts:

$$\dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} f_j^{\alpha(A)} l_i^{\alpha} + \frac{1}{V} \sum_{\alpha=1}^{N} f_j^{\alpha(SE)} l_i^{\alpha} \quad (27)$$

The second part of Eq. (27) represents the stress induced by surface-energy forces denoted as

$$(\dot{\sigma}_{ij})^{\text{SE}} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha(\text{SE})} l_{i}^{\alpha} \quad (28)$$

It is noted that this term is not analogous to the usual concept of cohesion for bulk material. In the present model, the surface-energy stress depends on the packing structure and is a tensor rather than a scalar. Only for an isotropic distribution of the branch lengths $l^\alpha$, the surface-energy stress can be reduced to an isotropic tensor.

As mentioned earlier, the surface-energy force is a function of contact area, which is in turn a function of interparticle force. Thus all equations in the second section are nonlinear in nature. The set of nonlinear equations for the stress-strain relationship is discussed in next section.

## Computation Scheme

The problem is defined by the following. Initially, we know the global variables ($\sigma_{ij}$ and $\varepsilon_{ij}$) for the assembly and the local variables ($f_j^\alpha$ and $\delta_j^\alpha$) for each contact orientation. For a given loading increment, which can be stress control, strain control, or mixed mode, out of the 12 variables ($\dot{\sigma}_{ij}$ and $\dot{\varepsilon}_{ij}$), 6 of them are unknown. The objective is to determine all global variables ($\sigma_{ij}$ and $\varepsilon_{ij}$) and local variables ($f_j^\alpha$ and $\delta_j^\alpha$) at the end of load increment. For a system with $N$ interparticle orientations, the number of unknown is $3N$ for $f_j^\alpha$ and $3N$ for $\dot{\delta}_j^\alpha$. The total number of unknown is $3N+3N+6$.
The following constraints must be satisfied:

1. The local constitutive equation, i.e., Eq. (22). As there are three equations for each contact plane orientation, the total number of equations is $3N$, $N$ being the total number of interparticle orientations.

2. Static hypothesis between global stress and local forces, i.e., Eq. (25): the number of equations is $3N$.

3. Strain definition between global strain and local displacement, i.e., Eq. (23). The number of equations is 6 (strain is symmetric).

The total number of unknowns is the same as the total number of equations. Therefore, the solution can be determined.

Using Eqs. (22), (23), and (25), the following relationship between stress and strain can be obtained:

$$ \dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp} \quad (29) $$

where

$$ C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{ep})^{-1} l_{k}^{\alpha} l_{n}^{\alpha} $$

The summation in Eq. (29) can be replaced by an integral over orientations. The integral can lead to a closed-form solution for the elastic modulus of randomly packed equal-size particles (Chang et al. 1995). However, in the elastic plastic behavior, due to the nonlinearity nature of the local constitutive equations, a numerical calculation with iterative process is necessary to carry out the summation in Eq. (29). In order to facilitate the numerical calculation, the orientations are selected to coincide with the locations of Gauss integration points in a spherical coordinate. Summation over these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over orientations. Results were found to be more accurate by using a set of fully symmetric integration points. From a study of the performance of using different numbers of orientations, we found $N \ge 74$ to be adequate (Chang and Hicher 2005).

For a strain-controlled test, Eq. (29) is not useful especially at the after-peak range of strain softening. In this case, a method of “elastic predictor-plastic corrector” was adopted to obtain the solution. For a mix-mode loading condition, an additional process of distributing the unbalance stresses was also needed (Chang and Hicher 2005).

## Summary of Parameters

One can summarize the material parameters as:

1. Normalized contact number per unit volume: $Nl^3/V$.

2. Mean particle size, $d$:

* Interparticle elastic constants: $k_{n0}$, $k_{r0}$, and $n$;

* Interparticle friction angle: $\phi_{\mu}$ and $m$;

* Interparticle hardening rule: $k_{p0}$ and $\phi_0$; and

* Critical state for the packing: $\lambda$ and $\Gamma$ or $e_{ref}$ and $p_{ref}$.

Besides critical state parameters, all other parameters are for interparticles. Standard values for $k_{p0}$ and $\phi_0$ are the following: $k_{p0}=k_{n0}$ and $\phi_0=\phi_{\mu}$ and a typical ratio $k_{r0}/k_{n0}=0.4$ can be generally assumed. Therefore, only six parameters have to be derived from experimental results and they can all be determined from the stress-strain curves obtained from drained compression triaxial tests.

Regarding the effect of the particle size, Eq. (15) shows that the interparticle stiffness $k_{n0}$ is influenced by the mean particle size. As in this model, the tangential stiffness $k_{r0}$ and the plastic stiffness $k_{p0}$ are directly related to $k_{n0}$, they are also particle size dependent.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Gradation curve for fine Hostun sand</figcaption></figure>

Considering an assembly of particles under an externally applied stress, according to Eqs. (14) and (15) and using the expression of contact force in terms of external stress, the interparticle stiffness is proportional to the mean particle size (see Liao et al. 2000):

$$ k_n \propto G_g^{1-n} d \left(\frac{Nl^3}{V}\right)^{-n} \sigma^n \quad (30) $$

Although the magnitude of the interparticle stiffness is dependent on the particle size, the elastic shear modulus of a packing is not particle-size dependent. This can be derived by the fact that the shear modulus of a packing is proportional to $k_n/d$, thus it leads to (Liao et al. 2000)

$$ G \propto G_g^{1-n} (\sigma)^n \left(\frac{Nl^3}{V}\right)^{1-n} \quad (31) $$

Note that the dimensionless density factor $Nl^3/V$ is a function of the void ratio only. Thus, the shear modulus of the packing is not particle-size dependent.

## Results of Numerical Simulation

### Triaxial Tests on Terrestrial Soil without Surface Energy Forces

A series of drained triaxial tests on fine Hostun Sand were performed by Al Mahmoud (1999). The gradation curve tested sand is shown in Fig. 4. The mean size of the particles for fine Hostun sand is $d=0.5$ mm. It is classified as uniform fine sand. The minimum and maximum void ratio for this sand is $e_{min}=0.575$ and $e_{max}=0.943$, respectively.

It is assumed here that surface-energy forces can be neglected under terrestrial conditions. This assumption will be discussed in the next section. The model needs a small number of input parameters, such as mean particle size, particle stiffness, interparticle friction, initial porosity tensor, and an initial degree of interlocking. The interparticle elastic constant $k_{n0}$ is assumed as 61,000 N/mm. Considering the grain size distribution curve, the contact numbers per volume $Nl^3/V$ are 2 for dense sand and 0.9 for loose sand (Hicher and Chang 2005). The evolution of $Nl^3/V$ due to new contact generation with neighboring particles was not considered.

The value of $k_{r0}/k_{n0}$ is commonly about 0.4, corresponding to a Poisson's ratio for Hostun Sand $\nu=0.2$ and the exponent $n=0.5$ (Biarez and Hicher 1994). From test results, we were able to
<table><caption>**Table 1.** Model Parameters for Fine Hostun Sand</caption><thead><tr><th>e<sub>ref</sub></th><th>P<sub>ref</sub><br>(MPa)</th><th>λ</th><th>φ<sub>μ</sub><br>(deg)</th><th>φ<sub>0</sub><br>(deg)</th><th>m</th><th>k<sub>p0</sub>/k<sub>n</sub></th></tr></thead><tbody><tr><td>1.05</td><td>0.01</td><td>0.14</td><td>34</td><td>34</td><td>1.5</td><td>0.02</td></tr></tbody></table>

derive the values of the two parameters corresponding to the position of the critical state in the $e-p'$ plane: $λ=0.14$ and $P_{ref}=0.01$ MPa for $e_{ref}=1.05$. The friction angle $φ_μ$ was also determined from the stress state corresponding to the critical state: $φ_μ=34°$. The equation governing the dilatancy rate requires the determination of the parameter $φ_0$. This parameter represents the concept of “phase transformation” as defined by Ishihara and Towhata (1983) or “characteristic state” as defined by Luong (1980). A value of $φ_0=φ_μ$ was retained.

The peak friction angle is not an intrinsic parameter, but varies with the void ratio according to Eq. (21). A value of $m=1.5$ was determined from the test results. The values of $k_{p0}$ are directly connected to the elastic properties through the relation: $k_{p0}=k_n$. The set of parameters for fine Hostun sand is presented in Table 1.

The tests were performed at different confining pressures on samples prepared at different initial void ratios. Typical results are presented in Figs. 5 and 6, which show the triaxial testing results for both dense and loose specimens made of Hostun sand. The stress–strain curves are plotted for dense sand in Fig. 5, and loose sand in Fig. 6. The model performance can be demonstrated by comparing the predicted and measured macrobehavior.

One can see the combined influence of the initial void ratio and of the mean effective stress on the stress–strain curves and the volumetric change. The stress–strain curve has a peak, which increases with confining stress. For dilatant materials, the deviatoric stress reduces after peak and converges toward a constant state of stress corresponding to the critical state. In practice, it is difficult to reach this state because of strain localization, especially in dense materials. The critical state was estimated from the results on loose specimens.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Stress–strain curves for loose Hostun sand</figcaption></figure>

## Triaxial Tests with Consideration of Surface Energy Forces

In order to evaluate the possible effect of lunar environmental conditions on the stress–strain response of soils, the influence of the surface-energy forces for different values of interparticle distance $D$ was examined. The stress–strain curves were predicted using the same soil parameters given in Table 1. Given the introduction of van der Waals forces, shear strength is expected to be higher. Fig. 7 shows the influence of the distance $D$ between the particles on the shear strength for soil specimens under 20 kPa confining stress, where $q_0$ represents the shear strength for soil without surface-energy forces and $q$ represents the shear strength with the effect of these forces.

It is noted that when the distance $D$ is greater than 2 nm, the effect of the van der Waals forces is practically negligible. This distance $D$ is a function of the amount of molecules, which can be adsorbed on the solid surface. The adsorbate thickness may be estimated by means of potential theory (Adamson 1990; Perko et al. 2001). Molecule adsorption on a solid surface is conditioned by temperature, gas pressure, and atmospheric composition. On Earth, adsorption conditions are easily fulfilled due to the presence of a high atmospheric pressure. Therefore, the distance $D$ can easily exceed several nanometers such that the effect of the

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Stress–strain curves for dense Hostun sand</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Effect of distance between particles on shear strength</figcaption></figure>

van der Waals force is negligible. Under lunar atmospheric con-
dition, the thickness of the adsorbed molecule layer is likely to be
very thin. For the case of 0.3 nm thickness (D=0.6 nm), the com-
parisons of the stress–strain curves in the Earth’s and the Moon is
atmosphere environments are given in Fig. 8 for dense soil and
Fig. 9 for loose soil under a 20 kPa confining stress. Assuming
the thickness D is 0.3–0.6 nm, the computed shear strength for

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Predicted stress–strain curves for lunar and terrestrial dense soils</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Predicted stress–strain curves for lunar and terrestrial loose soils</figcaption></figure>

<figure><img src="image_14.png" /></figure>

lunar soils is 12–15% higher than that of the Earth’s soil. The stiffness of the lunar soil is also higher. By plotting the maximum shear strength envelope in the Mohr plane for specimens under different confining stresses, the effect of van der Waals forces for the dense sand would result in an increase of 2–3 kPa in cohesion, whereas an increase of 0.5 kPa in cohesion is sufficient to hold a 1 m trench cut in lunar gravity condition.

**Summary and Conclusion**

In the present study we introduce surface-energy forces between
particles and estimate its effect on the shear strength of soil. The
microstructural approach is adopted because it allows us to model
interactions at particle level.

In the model, a simple elastic-plastic behavior was assumed on
each contact plane. The elastic part is based on the Hertz-
Mindlin's contact formulation, whereas the plastic part is based
on a Mohr-Coulomb friction law with an isotropic hardening as-
sumption and a nonassociated flow rule. The interaction among
interparticle planes is assumed to be dependent on the degree of
interlocking and porosity of the assembly of soil. The effect of
porosity is modeled by a phenomenological approach using the
concepts of critical state such that strain softening behavior can be modeled for dense materials. On the whole, the model requires a limited number of parameters, which can easily be determined from conventional triaxial tests.

Ability of the model to reproduce the main features of sand behavior has been demonstrated. Model simulations were compared with drained triaxial tests results at different initial void ratios and different confining stresses leading to contractant or dilatant behavior of the sand specimens. The comparison demonstrated that the model is capable of reproducing the general trend for both loose and dense sands.

The model was applied to predict the shear strength of lunar soil under extremely low atmospheric pressure. Model predictions indicate that soil under extremely low atmospheric pressure has an increase of shear strength by 12–15% higher than the one which would be present under Earth’s atmospheric pressure. The magnitude is in the same order as the measured increase of shear strength for lunar soil simulants tested under the usual atmospheric pressure and under a chamber with ultrahigh vacuum. However, the predicted shear strength for lunar soil due to surface-energy forces is only a hypothesis, which must await further testing on the Moon for validation/calibration.

## References

Al Mahmoud, M. (1999). "Etude en laboratoire du comportement des sables sous faibles contraintes." *Thèse de doctorat, l'Université des Sciences et Technologies de Lille, France* (in French).

Adamson, A. W. (1990). *Physical chemistry of surfaces*, 5th Ed., Wiley, New York.

Biarez, J., and Hicher, P. Y. (1994). *Elementary mechanics of soil behaviour, Balkema, Rotterdam, The Netherlands*.

Bromwell, L. G. (1966). "The friction of quartz in high vacuum." *Res. in Earth Physics Phase Rep. No. 7, R66-18*, Dept. of Civ. Engrg, Massachusetts Institute of Technology, Cambridge, Mass.

Cambou, B., Dedecker, F., and Chaze, M. (2000). "Relevant local variables for the change of scale in granular materials." *Constitutive modelling of granular materials*, D. Kolymbas, ed., Springer, Berlin, 275–290.

Chang, C. S. (1988). "Micromechanical modeling of constructive relations for granular material." *Micromechanics of granular materials*, M. Satake and J. T. Jenkins, eds., El Sevier Science Publications, Amsterdam, 271–279.

Chang, C. S., Chao, S. C., and Chang Y. (1995). "Estimates of mechanical properties of granulates with anisotropic random packing structure." *Int. J. Solids Struct.*, 32(14), 1989–2008.

Chang, C. S., and Gao, J. (1996). "Kinematic and static hypotheses for constitutive modelling of granulates considering particle rotation." *Acta Mech.*, 115(1–4), 213–229.

Chang, C. S., and Hicher, P.-Y. (2005). "An elasto-plastic model for granular materials with microstructural consideration." *Int. J. Solids Struct.*, 42(14), 4258–4277.

Chang, C. S., Kabir, M., and Chang, Y. (1992). "Micromechanics modelling for the stress strain behavior of granular soil—II: Evaluation." *J. Geotech. Engrg.*, 118(12), 1975–1994.

Chang, C. S., and Kuhn, M. R. (2005). "On virtual work and stress in granular media." *Int. J. Solids Struct.*, 42(13), 3773–3793.

Chang, C. S., and Liao, C. (1990). "Constitutive relations for particulate medium with the effect of particle rotation." *Int. J. Solids Struct.*, 26(4), 437–445.

Chang, C. S., and Ma, L. (1992). "Elastic material constants for isotropic granular solids with particle rotation." *Int. J. Solids Struct.*, 29(8), 1001–1018.

Chang, C. S., and Misra, A. (1990). "Application of uniform strain theory to heterogeneous granular solids." *J. Eng. Mech.*, 116(10), 2310–2328.

Chang, C. S., Misra, A., and Acheampon, K. (1992b). "Elastoplastic deformation of granulates with frictional contacts." *J. Eng. Mech.*, 118(8), 1692–1708.

Chang, C. S., Sundaram, S. S., and Misra, A. (1989). "Initial moduli of particulate mass with frictional contacts." *Int. J. Numer. Analyt. Meth. Geomech.*, 13(6), 626–641.

Christoffersen, J., Mehrabadi, H. M., and Nemat-Nasser, S. (1981). "A micromechanical description of granular material behavior." *J. Appl. Mech.*, 48(2), 339–344.

Costes, N. C., Carrier, W. D., Mitchell, J. K., and Scott, R. F. (1970). "Apollo 11: Soil mechanics results." *J. Soil Mech. and Found. Div.*, 96(6), 2045–2080.

Dahneke, B. (1972). "Influence of flattening on adhesion of particles." *J. Colloid Interface Sci.*, 40(1), 1–13.

Derjaguin, B. (1934). "Analysis of friction and adhesion. Part IV. The theory of the adhesion of small particles." *Kolloid-Z.*, 69(2), 155–164.

Derjaguin, B. V., Muller, V. M., and Toporov, Y. P. (1975). "Effect of contact deformations on adhesion of particles." *J. Colloid Interface Sci.*, 53(2), 314–326.

Emeriault, F., and Cambou, B. (1996). "Micromechanical modelling of anisotropic non-linear elasticity of granular medium." *Int. J. Solids Struct.*, 33(18), 2591–2607.

Goddard, J. D., and Bashir, Y. M. (1990). "On Reynolds dilatancy." *Recent development in structured continua*, D. De Kee and P. N. Kaloni, eds., Vol. II, Longman's, London, 23–35.

Hicher, P.-Y., and Chang, C. S. (2005). "Evaluation of two homogenization techniques for modelling the elastic behaviour of granular materials." *J. Eng. Mech.*, 131(11), 1184–1194.

Ishihara, K., and Towhata, I. (1983). "Cyclic behavior of sand during rotation of principal axes." *Mechanics of granular materials*, Elsevier, New York, 55–73.

Israelachvili, J. (1992). *Intermolecular and surface forces*, 2nd Ed., Academic, San Diego.

Jenkins, J. T. (1988). "Volume change in small strain axisymmetric deformations of a granular material." *Micromechanics of granular materials*, M. Satake and J. T. Jenkins, eds., El Sevier Science Publications, Amsterdam, 143–152.

Johnson, K. L. (1971). "Surface energy and contact of elastic solids." *Proc. R. Soc. London, Ser. A*, 324(1558), 301–313.

Kruyt, N. P. (2003). " Statics and kinematics of discrete Cosserat-type granular materials." *Int. J. Solids Struct.*, 40, 511–534.

Kruyt, N. P., and Rothenburg, L. (2002). "Micromechanical bounds for the effective elastic moduli of granular materials." *Int. J. Solids Struct.*, 39(2), 311–324.

Liao, C. L., Chan, T. C., Suiker, A. S. J., and Chang, C. S. (2000). "Pressure-dependent elastic moduli of granular assemblies." *Int. J. Numer. Analyt. Meth. Geomech.*, 24(3), 265–279.

Luding, S., Latzel, M., Volk, W., Diebels, S., and Herrmann, H. J. (2001). "From discrete element simulations to a continuum model." *Comput. Methods in Appl. Mech. Eng. Elsevier*, 191(1–2), 21–28.

Luong, M. P. (1980). "Stress-strain aspects of cohesionless soils under cyclic and transient loading." *Int. Symp. on Soils under Cyclic and Transient Loading*, 353–376.

Mindlin, R. D. (1969). "Microstructure in linear elasticity." *Arch. Ration. Mech. Anal.*, 16(1), 51–78.

Mitchell, J. K. (1972). "Soil mechanical properties at the Apollo 14 site." *J. Geophys. Res.*, 77(29), 5641–5664.

Molerus, O. (1978). "Effect of interparticle cohesive forces on flow Behavior of Powders." *Powder Technol.*, 20, 161–175.

Nelson, J. D. (1967). "Environmental effects on engineering properties of simulated lunar soils." Ph.D. thesis, Illinois Institute of Technology, Chicago.

Perko, H., Nelson, J., and Sadeh, W. (2001). "Surface cleanliness effect on lunar soil shear strength." *J. Geotech. Geoenviron. Eng.*, 127(4), 371–383.
Rothenburg, L., and Selvadurai, A. P. S. (1981). "Micromechanical definitions of the Cauchy stress tensor for particular media." *Mechanics of structured media*, A. P. S. Selvadurai, eds., Elsevier, Amsterdam, The Netherlands, 469–486.

Rowe, P. W. (1962). "The stress-dilatancy relations for static equilibrium of an assembly of particles in contact." *Proc. R. Soc. London, Ser. A*, 269(1339), 500–527.

Scott, R. F., (1973). "Lunar soil mechanics." *Proc., 8th Int. Congress On Soil Mechanics and Foundations Engineering*, Springer, New York, 177–190.

Suiker, A. S. J., and Chang, C. S. (2004). "Modelling failure and deformation of an assembly of spheres with frictional contacts." *J. Eng. Mech.*, 130(3), 283–293.

Taylor, D. W. (1948). *Fundamentals of soil mechanics*, Wiley, New York.

Tejchman, J., and Bauer, E. (2005). "Modeling of a cyclic plane strain compression-extension test in granular bodies within a polar hypoplasticity." *Granular Matter*, 7(4), 227–242

Thornton, C., and Ning, Z. M. (1998). "A theoretical model for the stick/ bounce behaviour of adhesive, elastic-plastic spheres." *Powder Technol.*, 99(2), 154–162.

Tomas, J. (2001). "Assessment of mechanical properties of cohesive particulate solids. Part 1: Particle contact constitutive model." *Part. Sci. Technol.*, 19(2), 95–110.

Valverde, J. M., Castellanos, A., and Watson, P. K. (2001). "The effect of particle size on interparticle adhesive forces for small loads." *Powder Technol.*, 18(3), 236–241.

Wallton, K. (1987). "The effective elastic moduli of a random packing of spheres." *J. Mech. Phys. Solids*, 35(22), 213–226.