# A microstructural elastoplastic model for unsaturated granular materials

Pierre-Yves Hicher, Ching S. Chang

► To cite this version:

Pierre-Yves Hicher, Ching S. Chang. A microstructural elastoplastic model for unsaturated granular materials. International Journal of Solids and Structures, 2007, 44 (7-8), pp.2304-2323. 10.1016/j.ijsolstr.2006.07.007. hal-01004848

HAL Id: hal-01004848

https://hal.science/hal-01004848

Submitted on 8 Mar 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
# A microstructural elastoplastic model for unsaturated granular materials

P.-Y. Hicher<sup>a</sup>, C.S. Chang<sup>b</sup>

The homogenization technique is used to obtain an elastoplastic stress–strain relationship for dry, saturated and unsaturated granular materials. Deformation of a representative volume of material is generated by mobilizing particle contacts in all orientations. In this way, the stress–strain relationship can be derived as an average of the mobilization behavior of these local contact planes. The local behavior is assumed to follow a Hertz–Mindlin’s elastic law and a Mohr–Coulomb’s plastic law. For the non-saturated state, capillary forces at the grain contacts are added to the contact forces created by an external load. They are calculated as a function of the degree of saturation, depending on the grain size distribution and on the void ratio of the granular assembly. Numerical simulations show that the model is capable of reproducing the major trends of a partially saturated granular assembly under various stress and water content conditions. The model predictions are compared to experimental results on saturated and unsaturated samples of silty sands under undrained triaxial loading condition. This comparison shows that the model is able to account for the influence of capillary forces on the stress–strain response of the granular materials and therefore, to reproduce the overall mechanical behavior of unsaturated granular materials.

**Keywords:** Unsaturated granular material; Stress–strain relationship; Capillary forces; Microstructure; Elastoplastic model

## 1. Introduction

Macroscopic properties of granular materials are governed by the properties of inter-particle contacts. For dry granular materials, the inter-particle forces are related to the applied external stresses. When a certain amount of water is added to the grain assembly, significant changes, such as an increase in stiffness and strength, can be observed. These changes can be explained by the formation of water menisci between neighboring particles, which creates capillary forces. The amplitude of these capillary forces depends upon the degree of saturation and the morphology at particle level (grain and pore sizes) (Hicher, 1998).

<footnote><sup>a</sup> Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes, University of Nantes, France</footnote>

<footnote><sup>b</sup> Department of Civil and Environmental Engineering, University of Massachusetts, Amherst, MA 01002, United States</footnote>
Macroscopic interpretations of the mechanical behavior of unsaturated soils have been developed mainly in the framework of elastoplasticity. Most of these models reject the concept of effective stress and consider that the strain tensor is governed by the total stress tensor $\sigma_{ij}$ (or $\sigma_{ij} - p_g\delta_{ij}$, $p_g$ being the gas pressure) and the negative pore water pressure or suction $s$ (Alonso et al., 1990; Wheeler and Sivakumar, 1995; Cui and Delage, 1996). In particular, they consider a new yield surface, called loading-collapse (LC) surface in the plane $\sigma_{ij} - p_g\delta_{ij}$, $s$ which controls the volume changes due to the evolution of the degree of saturation for a given loading condition. Other models introduce a capillary stress $\sigma_c$, function of the suction $s$, so that the following relation can be used to defined the concept of a generalized effective stress $\sigma'$ (Biarez et al., 1993; Kohgo et al., 1993; Modaressi and Abou-Bekr, 1994):

$$ \sigma' = \sigma - p_g \mathbf{I} - \sigma_c. \qquad (1) $$

The above equation can be considered as a general formulation of the equation proposed by Bishop and Blight (1963) for unsaturated soils:

$$ \sigma' = \sigma - p_g \mathbf{I} + X(S_1)(p_g - p_1)\mathbf{I}, \qquad (2) $$

where $p_1$ is the liquid pressure and $X(S_1)$ is a function of the degree of saturation $S_1$ of the liquid phase ($X=0$ for a dry material, $X=1$ for a fully saturated material).

An alternative method is to develop homogenization techniques in order to derive stress-strain relationship from forces and displacements at the particle level. The basic idea is to view the packing as represented by a set of micro systems. The inelastic behavior of each micro-system is characterized and the overall stress-strain relationship of packing is obtained from an average of the micro-systems behaviors. Micro-systems can be regarded as inter-particle planes (or mobilized planes) in the packing. Models based on inter-particle contact planes can be found in Jenkins and Strack (1993), Matsuoka and Takeda (1980) and Chang et al. (1989a), etc. An alternative way is to view the micro-systems as particle groups of different configurations. Models developed along this line can be found in Chang et al. (1992a,b), and more recently in Suiker and Chang (2004).

In order to characterize the packing structure, the microstructural plasticity model requires a mathematical or numerical description of the micro systems, which is not readily available. But, in cases of treating micro systems as contact or mobilized planes, the approach of estimating the overall behavior by orientational planes can be linked to G.I. Taylor's concept, developed long ago in the slip theory of plasticity for polycrystalline materials by Batdorf and Budianski (1949). These ideas were applied by Pande and Sharma (1982) to rocks and soils in what they called the overlay model, and to concrete by Bazant et al. (1995) in the so-called microplane model.

Along these lines, we recently developed a new stress-strain model which considers inter-particle forces and displacements. We incorporated into the model a micro-macro links between strain and inter-particle displacements for granular materials (Liao et al., 1997). By comparing predicted and measured triaxial loading results for sands of different void ratios and under different confining stresses, in both drained and undrained conditions, we demonstrated the ability of this model to reproduce the main features of fully saturated or dry sand behavior (Chang and Hicher, 2005). In this paper, we have extended the capability of the model to unsaturated states by incorporating the influence of capillary forces at the micro level.

## 2. Stress-strain model

In this model, we envision a granular material as a collection of particles. The deformation of a representative volume of this material is generated by the mobilization of contact particles in all orientations. Thus, the stress-strain relationship can be derived as an average of the mobilization behavior of local contact planes in all orientations. For a contact plane in the $\alpha$th orientation, the local forces $f_i^\alpha$ and the local movements $\delta_i^\alpha$ can be denoted as follows: $f_i^\alpha = \{f_n^\alpha, f_s^\alpha, f_t^\alpha\}$ and $\delta_i^\alpha = \{\delta_n^\alpha, \delta_s^\alpha, \delta_t^\alpha\}$, 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 (Fig. 1).

The forces and movements at the contact planes of all orientations 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 macro stress and strain must be equivalent to that expressed in
<figure><img src="image_1.png" /><figcaption>Fig. 1. Local coordinate at inter-particle contact.</figcaption></figure>

terms of microforces and movements. In such a formulation, it has usually been assumed that the microstructure is statically constrained, which means that the forces acting 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 generally preferred for 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 obviously 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).

## 2.1. Inter-particle behavior

### 2.1.1. Contact forces and capillary forces

For dry samples, the contact forces are directly determined from the external stresses $\sigma$ applied on the granular assembly (Eq. (21)). In the case of wetted samples, different stages of saturation can be identified. The fully saturated regime corresponds to a two-phase material with water filling completely the voids between grains. The water pressure $p_w$ can be positive or negative (suction) but in both cases the effective stress concept (Terzaghi, 1925) can be applied and the contact forces determined by considering the effective stresses $\sigma'$ as the external stresses (de Buhan and Dormieux, 1996; Hicher, 1998):

$$ \sigma' = \sigma - p_w \mathbf{I}. \qquad (3) $$

In the case of partially saturated samples, the liquid phase is distributed in menisci located between close grains and as a consequence, capillary forces are applied on the grains and are added to the contact forces defined above. When the water content decreases inside a saturated granular sample, the air breaks through at a given point. The negative water pressure or suction corresponding to that point is called the air-entry pressure. Its value depends on the pore sizes. Afterwards, the sample becomes unsaturated with the water phase being connected inside the pores. This state is called the funicular regime, in which a constant decrease of the degree of saturation corresponds to a gentle increase in pore water pressure, which is homogenized along the continuous water phase. The menisci start to form between two grains, not necessarily in contact, and are connected to each other. The pendular regime starts when the water menisci become disconnected. The water is no longer continuous and equilibrium in the water pressure is obtained by the vapor pressure. At this stage, the capillary forces start to increase drastically according to the pore size. If the drying process continues, the water bridges begin to fail, starting with the non-contacting grains, until a completely dry state is achieved and no capillary forces are left inside the granular assembly.

Several studies have shown that the attractive capillary forces between two grains connected by a water bridge are a decreasing function of the distance between the grains until the bridge fails (see for example Taibi, 1994; Chateau et al., 2002). This function depends on the volume of liquid found between the grains. Different
mathematical expressions have been proposed for these capillary forces, which are the sum of the pressure forces exerted by the liquid on the wetted area and the surface tension forces acting along the wetted perimeter. In this study, we retained the following expression:

$$f_n^{\text{cap}} = f_{\text{max}}e^{-c(\frac{d}{R})}, \qquad (4)$$

where $f_{\text{cap}}$ is the capillary force between two neighboring grains, not necessary in contact, $f_{\text{max}}$ is the value of $f_{\text{cap}}$ for two grains in contact, and $R$ is the mean grain radius, $d$ represents the distance between two grains and is equal to $l - 2R$, $l$ being the branch length given as a distribution function of the grain size and the void ratio, $c$ is a material parameter, dependent on the grain morphology and on the water content, $f_{\text{max}}$ depends on the capillary pressure defined as the pressure jump across the liquid-air interface, on the liquid-air interface surface tension, as well as on the geometry of the menisci governed by the solid-liquid contact angle and the filling angle. In this study, a simplified approach was to consider an empirical relation between $f_{\text{max}}$ and the degree of saturation $S_r$, without taking into account the hysteresis along drying and wetting paths (Hicher, 1998)

$$f_{\text{max}} = f_0 \frac{S_r}{S_0} \quad \text{for } 0 < S_r < S_0, \\ f_{\text{max}} = f_0 \frac{S_0(1 - S_r)}{S_r(1 - S_0)} \quad \text{for } S_0 < S_r < 1, \qquad (5)$$

where $f_0$ and $S_0$ are material parameters. $f_0$ depends on the grain size distribution, $S_0$ represents the degree of saturation at which any further drying of the specimen will cause substantial breaking of the menisci in the pendular domain. $S_0$ depends on the nature of the granular material. The following empirical expression was proposed by Wu et al. (1984) for compacted granular materials prepared at various initial degrees of saturation:

$$(S_r)_{\text{opt}} = (-0.65 \log(d_{10}) + 1.5)/100, \qquad (6)$$

in which $(S_r)_{\text{opt}}$ is equivalent to $S_0$ in Eq. (5), $d_{10}$ is the effective grain size in mm. One should note that Eq. (6) could be dependent on the mode of preparation. Cho and Santamaria (2001) showed that, when the specimens were prepared in a fully saturated state and then subjected to gradual drying, $(S_r)_{\text{opt}}$ could be substantially lower than the value obtained by Eq. (6). For example, in the case of an assembly of mono-sized glass beads ($d = 0.32$ mm), they obtained a value of $(S_r)_{\text{opt}} = 0.007$, followed by a sharp decrease of the capillary effect as drying continued until $S_r = 0$. This could be due to the mode of preparation, which is different in the two cases and could lead to different specimen states at small degrees of saturation. For specimen prepared at a given low water content, menisci might not form between two wet particles because the water film is not thick enough, whereas during drying, the menisci may subsist longer before starting to fail when the water continues to evaporate.

Since the menisci are not necessarily all formed in the funicular regime, Eq. (2) may not be applicable for high degrees of saturation. However, in this first approach, we decided to extend it to the whole range of saturation, considering that the amplitudes of capillary forces were small for degrees of saturation higher that 80% and could therefore be approached with sufficient accuracy by using the same equation.

### 2.1.2. Elastic relationship

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. 1. The contact stiffness of a contact plane includes normal stiffness, $k_n^z$, and shear stiffness, $k_r^z$. The elastic stiffness tensor is defined by

$$f_i^z = k_{ij}^{ze} \delta_j^{ze}, \qquad (7)$$

which can be related to the contact normal and shear stiffness

$$k_{ij}^{ze} = k_n^z n_i^z n_j^z + k_r^z (s_i^z s_j^z + t_i^z t_j^z), \qquad (8)$$
where n, s, t are three orthogonal unit vectors that form the local coordinate system. Vector n is outward normal to the contact plane. Vectors s and t are on the contact plane. For the local coordinate shown in Fig. 1, the normal vectors are

$$ 
\begin{aligned}
n_1 &= \cos \gamma, \\
n_2 &= \sin \gamma \cos \alpha, \\
n_3 &= \sin \gamma \sin \alpha.
\end{aligned} 
\qquad (9) $$

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

$$ k_n = k_{n0} \left( \frac{f_n}{G_g l^2} \right)^n; \quad k_t = k_{t0} \left( \frac{f_t}{G_g l^2} \right)^n, \qquad (10) $$

where $G_g$ is the elastic modulus for the grains, $f_n$ is the contact force in normal direction, $l$ is the branch length between the two particles, $k_{n0}$, $k_{t0}$ and $n$ are material constants. For two spherical particles, the branch length is same as particle size $l=d$. Let $n=1/3$, and

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

Eq. (10) is equivalent to the Hertz-Mindlin's contact formulation.

### 2.1.3. Plastic relationship

The elastic sliding behavior between two particles does not contain a coupling effect (i.e., there is no shear induced normal movements). However, plastic sliding often occurs along the tangential direction of the contact plane with an upward or downward movement, thus shear dilation/contraction takes place. The stress-dilatancy is a well-known phenomenon in sand (see discussions in the work by Taylor, 1948; Rowe, 1962; Goddard and Bashir, 1990, etc.), and should be correctly modeled. The dilatancy effect can be described by

$$ \frac{d\delta_n^P}{d\Delta^P} = \frac{T}{f_n} - \tan \phi_0, \qquad (12) $$

where $\phi_0$ is a material constant which, in most cases, can be considered equal to the internal friction angle $\phi_{\mu}$. This equation can be derived by equating the dissipation work due to plastic movements to the friction loss of contact planes in the same orientation. Note that the shear force $T$ and the rate of plastic sliding $d\Delta^P$ can be defined as

$$ T = \sqrt{f_s^2 + f_t^2} \quad \text{and} \quad d\Delta^P = \sqrt{(d\delta_s^P)^2 + (d\delta_t^P)^2}. \qquad (13) $$

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

$$ F(f_i, \kappa) = T - f_n\kappa(\Delta^P) = 0, \qquad (14) $$

where $\kappa(\Delta^P)$ is an isotropic hardening/softening parameter. When $F > 0$, it indicates loading, otherwise unloading. The hardening function is defined by a hyperbolic curve in $\kappa - \Delta^P$ plane, which involves two material constants: $\phi_p$ and $k_{p0}$.

$$ \kappa = \frac{k_{p0} \tan \phi_p \Delta^P}{|f_n| \tan \phi_p + k_{p0} \Delta^P}. \qquad (15) $$

The value of $\kappa$ asymptotically approaches $\tan\phi_p$. The initial slope of the hyperbolic curve is $k_{p0}$. On the yield surface, under a loading condition, the shear plastic flow is determined by a normality rule applied to the yield function. However, the plastic flow in the direction normal to the contact plane is governed by the stress-dilatancy equation in Eq. (12). Thus, the flow rule is non-associated.
### 2.1.4. Interlocking influence

One of the important elements to consider in granular modeling is the concept of critical state. Under critical state, the granular material will remain at a constant volume while it is subjected to a continuous distortion. The void ratio corresponding to this state is $e_c$.

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

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

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

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

$$\tan \phi_p = \left(\frac{e_c}{e}\right)^m \tan \phi_\mu, \qquad (17)$$

where $m$ is a 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.

### 2.1.5. Elasto-plastic relationship

With the elements discussed above, the final incremental stress-strain relations of the material can be derived that includes both elastic and plastic behavior, given by

$$\dot{f}_i^\alpha = k_{ij}^{\alpha p} \dot{\delta}_j^\alpha. \qquad (18)$$

Detailed expression of the elasto-plastic stiffness tensor is given in Chang and Hicher (2005).

## 2.2. Stress–strain relationship

### 2.2.1. Macro–micro relationship

The stress–strain relationship for an assembly can be determined from integrating the behavior of inter-particle contacts in all orientations. During the integration process, a relationship is required to link the macro and micro variables. Using the static hypotheses proposed by Liao et al. (1997), we obtain the relation between the macro strain and inter-particle displacement (finite strain condition not being considered here)

$$\dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \dot{\delta}_{j}^{\alpha} l_{k}^{\alpha}, \qquad (19)$$

where $\dot{\delta}_j$ is the relative displacement between two contact particles and the branch vector $l_k$ is the 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 contact orientations. The variables $\dot{\delta}_j^\alpha$ and $l_k^\alpha$ are defined respectively as the averaged values of $\dot{\delta}_j$ and $l_k$ for all contacts belong to the $\alpha$th orientation. The fabric tensor in Eq. (19) is defined as

$$A_{ik} = \sum_{\alpha=1}^{N} l_{i}^{\alpha} l_{k}^{\alpha}. \qquad (20)$$

If the magnitudes of the branch vectors $l_i^\alpha$ are different for each orientation $\alpha$, then the fabric tensor $A_{ik}$ is anisotropic. During the deformation of the assembly, the particles rearranged and the lengths of the branch vectors are changed into a new set of numbers, which is a way to allow the evolution of the strain induced anisotropy to be considered. One of the major advantage of the model is its capacity to take into account the initial and induced structural anisotropy.
Using both the principle of energy balance and Eq. (19), the mean force on the contact plane of each orientation is

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

In Eq. (21), the stress increment $\dot{\sigma}_{ij}$ can be obtained through the contact forces and branch vectors for contacts in all orientations (Christofferson et al., 1981; Rothenburg and Selvadurai, 1981). Since $l_k^\alpha$ represents the mean branch vector for the $\alpha$th orientation including both contact and non-contact particles, the value of $\dot{f}_j^\alpha$ in Eq. (21) represents the mean of contact force in the $\alpha$th orientation.

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

When the defined contact force is applied in Eqs. (21) and (22) is unconditionally satisfied. Using the definition of Eq. (22), the stress induced by capillary forces can be computed and is termed as capillary stress, given by

$$ (\dot{\sigma}_{ij})^{\text{cap}} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha(\text{cap})} l_{i}^{\alpha}. \qquad (23) $$

It is noted that this term is not analogous to the usual concept of capillary pressure or suction which represents the negative pore water pressure inside the unsaturated material. In this model, the capillary stress depends on the geometry of the pores and is a tensor rather than a scalar. Only for an isotropic distribution of the branch lengths $l_k^\alpha$, the capillary stress can be reduced to an isotropic tensor. This can be the case for an initially isotropic structure during isotropic loading, but during deviatoric loading, an induced anisotropy is created and the capillary tensor is no longer isotropic.

### 2.2.2. Computation scheme

The problem is defined as follows:

Initially, we know the global variables ($\sigma_{ij}$ and $\epsilon_{ij}$) for the assembly and the local variables ($f_j^\alpha$ and $\delta_j^\alpha$) for each contact orientation. For any given loading increment, which can be stress control, strain control or mixed mode, 6 out of the 12 variables ($\Delta\sigma_{ij}$ and $\Delta\epsilon_{ij}$) are unknown. The objective is therefore to determine all global variables ($\sigma_{ij}$ and $\epsilon_{ij}$) and local variables ($f_j^\alpha$ and $\delta_j^\alpha$) at the end of load increment. For a system with $N$ inter-particle orientations, the number of unknown is $3N$ for $f_j^\alpha$ and $3N$ for $\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. (18). Since there are three equations for each contact plane orientation, the total number of equations is $3N$: $N$ being the total number of inter-particle orientations.

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

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

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

Using Eqs. (18), (19), and (21), the following relationship between stress and strain can be obtained:

$$ \dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp}; \quad \text{where} \quad C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{ep})^{-1} l_k^{\alpha} l_n^{\alpha}. \qquad (24) $$

## 3. Summation and integration

The summation of a quantity $F^\alpha$ over all contacts can be replaced by an integral over all orientations, provided that the quantity can be expressed by a continuous function $F(\gamma, \beta)$. In granular packing, this is always the case

$$ \sum_{\alpha=1}^{N_c} F^{\alpha} = \frac{N_c}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) F(\gamma, \beta) \sin \gamma \, d\gamma \, d\beta. \qquad (25) $$
$N_c$ is the number of contact orientations $E(\gamma, \beta)$ is the density function of contacts. When the function $F(\gamma, \beta) = 1$, the following consistency condition must be satisfied for the density function of contacts

$$1 = \frac{1}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) \sin \gamma \, d\gamma \, d\beta. \qquad (26)$$

The density function can be obtained from a spherical harmonic expansion of the true distribution of contacts in discrete orientations. For a natural deposit, the contact distribution is symmetric about the major axis, thus the density function is only a function of $\gamma$ (see Chang and Misra, 1990) given by:

$$E(\gamma, \beta) = \left(1 + \frac{a_0}{4}(3 \cos 2\gamma + 1)\right). \qquad (27)$$

Eq. (27) is a truncated form of the spherical harmonic expansion up to second-order terms. For example, the integration form of the fabric tensor given in Eq. (20) is:

$$\sum_{\alpha=1}^{N_c} l_i^{\alpha} l_k^{\alpha} = \frac{N_c}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) l_i^{\alpha} l_k^{\alpha} \sin \gamma \, d\gamma \, d\beta. \qquad (28)$$

Assuming that the density function is in the form of Eq. (27), and that the branch lengths are equal in all directions, will therefore: $l_i^{\alpha} = ln_{i}^{\alpha}$. Using Eq. (9) for the normal vectors, the fabric tensor is:

$$[A_{ij}] = \frac{N_c l^2}{3} \begin{bmatrix} 1 + a_0 & 0 & 0 \\ 0 & 1 - \frac{a_0}{2} & 0 \\ 0 & 0 & 1 - \frac{a_0}{2} \end{bmatrix}. \qquad (29)$$

The integral of Eq. (25) can usually 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 non-linear nature of the local constitutive equations, a numerical calculation with an iterative process is necessary to carry out the summation. In order to facilitate the numerical calculation for the spherical integral shown in Eq. (25), it is convenient to discretize the domain into a selected set of Gauss integration points in a spherical coordinate, each point having a location $(\gamma, \beta)$ on the unit sphere, and having its associated weighting $w^{\alpha}$. Summing up these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over all orientations. We found that the results were more accurate by using a set of fully symmetrical integration points. From a study of the performance of using different numbers of orientations, we found 37 points to be adequate. Spherical integration can then be evaluated by

$$\int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) F(\gamma, \beta) \sin \gamma \, d\gamma \, d\beta = 4\pi \sum_{\alpha=1}^{37} E^{\alpha} F^{\alpha} w^{\alpha}, \qquad (30)$$

where superscript $\alpha$ represents the orientation of the $\alpha$th integration point. Thus, the values of $E^{\alpha} = E(\gamma, \beta)$ and $F^{\alpha} = F(\gamma, \beta)$ for each orientation can be easily evaluated. When $E(\gamma, \beta) = 1$ and $F(\gamma, \beta) = 1$ for all orientations, the integral represents the surface area of a unit sphere (i.e., $2\pi$). Thus, the weighting factors must satisfy

$$\sum_{\alpha=1}^{37} w^{\alpha} = \frac{1}{2}. \qquad (31)$$

Using these concepts, it is possible to make an equivalence between the summation of all contacts to the summation of the 37 orientations of the integration points:

$$\sum_{i=1}^{N_c} F^{\alpha} = 2N_c \sum_{i=1}^{37} E^{\alpha} F^{\alpha} w^{\alpha}. \qquad (32)$$

For a strain-controlled test, Eq. (24) 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. The detailed procedure is not included here.
### 3.1. Summary of parameters

One can summarize the material parameters as:

- Normalized contact number per unit volume: $Nl^3/V$;

- Mean particle size, $2R$;

- Inter-particle elastic constants: $k_{n0}$, $k_{t0}$ and $n$;

- Inter-particle friction angle: $\phi_\mu$ and $m$;

- Inter-particle hardening rule: $k_{p0}$ and $\phi_0$;

- Critical state for packing: $\lambda$ and $\Gamma$ or $e_{\text{ref}}$ and $p_{\text{ref}}$;

- Capillary force equation: $f_0$, $S_0$ and $c$.

Other than critical state parameters, all other parameters are inter-particles. Standard values for $k_{p0}$ and $\phi_0$ are the following: $k_{p0} = k_n$ and $\phi_0 = \phi_\mu$ and a typical ratio $k_{t0}/k_{n0} = 0.4$ can generally be assumed (Hicher and Chang, 2006). Therefore, for dry or saturated samples, only six parameters have to come from experimental results and these can all be determined from the stress-strain curves obtained from drained or undrained compression triaxial tests. For unsaturated states, three more parameters need to be determined, using specific triaxial tests on partially saturated samples.

The initial geometrical state (initial distribution of $l_i^\alpha$) is assumed to be isotropic for randomly packed granulates. Along a loading path, the computed length of branch vector $l_i^\alpha$ changes differently for each orientation, which is contributed from inter-particle deformation due to contractancy or dilatancy. Thus the evolution of the distance $d$, which depends on $l$, can be taken into account.

## 4. Results of numerical simulation

Based on triaxial test results, numerical simulations of experimental testing on two different granular materials were performed. Fig. 2 presents the grain size distribution curves for the tested materials, two silty sands. For each material, the mean particle size, $2R$, was considered as equal to the $d_{50}$ value.

To determine the capillary forces, Eq. (4) includes two material parameters $c$ and $d$. Very little information could be gathered on these two parameters by simply using the test results. The distance between two grains, $d$, was taken equal to $l - 2R$, $l$ being the branch length. An initial value of $l$, equal in all directions, was considered, so that the initial value of $d$ was equal to $0.01 d_{50}$, which leads to an averaging of the capillary forces inside the volume. The evolution of $d$ during loading follows therefore the evolution of $l$, which was previously discussed. The value of $c$ depends on the water content and the grain morphology. Several theoretical and experimental results allow us to indirectly determine a value of $c$ in the case of two identical spheres (see

<figure><img src="image_2.png" /><figcaption>Fig. 2. Grain size distribution curves.</figcaption></figure>
for example Chateau et al., 2002). It is however difficult to extrapolate these results for grains of various
shapes and sizes. Since most of the following simulations were performed at constant water content, on spec-
imens with degrees of saturation higher than 50%, a constant value of c = 8 was assumed, considered as rep-
resentative of the amount of water present in the volume during the tests. This hypothesis is less acceptable for
the last example, which consists of a wetting test. We assumed, however, that retaining this hypothesis would
not affect in a significant way the results obtained by the model. These values of c and d give a standard
evolution for the capillary forces, function of the distance between particles, as well as an initial isotropic
distribution of these forces.

4.1. Undrained triaxial tests. Influence of the confining stress

Well-documented experimental results on a loose decomposed granitic soil were presented by Ng and Chiu (2003). The grain size distribution is shown in Fig. 2. The soil can be described as a silty, clayey, gravelly sand. Triaxial samples were prepared by the wet tamping method, compacted with a water content of 19% and dry densities ranging from 74% to 84% of the Proctor maximum, which corresponded to a void ratio ranging from 0.85 to 1.1 and degrees of saturation ranging from 45 to 58%. Undrained triaxial tests on saturated samples were performed at different confining pressures and initial void ratios. Figs. 3 and 4 present the results obtained on three representative specimens. After preparation by compaction, the samples were saturated and then isotropically compressed to confining pressures ranging from 25 to 200 kPa. The void ratios before shearing ranged from 0.857 to 0.928.

The mean particle size 2R was selected equal to d₅₀ = 1 mm. The elastic parameters were not directly acces-
sible from these test results. From previous studies (Hicher, 1996; Hicher and Chang, 2006), typical values
were adopted in this study. The inter-particle elastic constant kₚ₀ was assumed to be equal to 60 000 N/mm
and kₜ₀/kₚ₀ = 0.4. The normalized contact number per volume Nl³/V could be obtained by using an empirical
relation with the void ratio. The test results presented in Figs. 3 and 4 and several other ones on saturated
samples at different void ratios and different confining pressures made it possible to determine the critical state
parameters: M = 1.55 (φμ = 38°), λ = 0.11 and eₐref = 1.08 for pₐref = 10 kPa. Table 1 summarizes the set of
parameters used for the modeling of the saturated sample behavior.

The comparison between experimental and numerical results for the three selected specimens is presented in Figs. 3 and 4. The specimens compressed at low confining pressure had an initial state located on the left side of the critical state line in the e.p' plane and showed therefore a tendency to dilate, which led to an increased effective mean stress. The deviatoric stress increased to up to the ultimate strength at large deformations, corresponding to the critical state. The specimen compressed at higher confining stress had an initial state located on the right side of the critical state and showed a tendency to contract, which led to a decrease of the effective mean stress and to a peak in the deviatoric stress, followed by a decrease toward a constant value at large

<figure><img src="image_3.png" /><figcaption>Fig. 3. Undrained tests on saturated samples of decomposed granitic soil: stress–strain relationship.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Undrained tests on saturated samples of decomposed granitic soil: effective stress path.</figcaption></figure>

<table><caption>Table 1
Model parameters for the decomposed granitic soil</caption><thead><tr><td>e<sub>ref</sub></td><td>p<sub>ref</sub> (kPa)</td><td>λ</td><td>φ<sub>μ</sub> (°)</td><td>φ<sub>0</sub> (°)</td><td>m</td></tr></thead><tbody><tr><td>1.08</td><td>10</td><td>0.11</td><td>38</td><td>38</td><td>1</td></tr></tbody></table>

deformations. The numerical simulations show that the model is able to match the general trend of the experimental behavior and that the quantitative results are in reasonably good agreement with the experimental ones.

After applying an initial suction of 40 kPa, unsaturated samples were isotropically compressed and then sheared at a constant water content. Figs. 5 and 6 present the results obtained for three different confining pressures, 50 kPa, 100 kPa and 150 kPa on specimens with the same initial void ratio $e_0 = 1.09$. The results show a constant increase of the deviatoric stress up to a constant value for axial strain over 25%. The deviatoric stress increased with the applied mean stress.

The three specimens exhibited a contractive volumetric behavior, with higher contraction for the specimen subjected to higher confining stress. A stabilization of the volume change could also be observed at large deformations. This contractant behavior led to an increase of the degree of saturation for all samples since the water content was kept constant. The initial degree of saturation was equal to 43% for the three specimens.

<figure><img src="image_5.png" /><figcaption>Fig. 5. Constant water content tests on unsaturated samples of decomposed granitic soil: stress-strain relationship.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 6. Constant water content tests on unsaturated samples of decomposed granitic soil: volumetric-axial strain relationship.</figcaption></figure>

and the final degree of saturation reached 48% to 61%, according to the amount of volume change during the tests. Measured suctions showed a decrease of 8 kPa to 15 kPa for the most contractive specimen (Fig. 9).

For numerical simulations, determining the material parameters $S_0$ and $f_0$ was the necessarily first step. $S_0$ was determined by using Eq. (6), considering that the mode of preparation by compaction was similar to the one used by Wu et al. (1984). $f_0$ was determined from the test result at the confining stress of 100 kPa, in order to obtain a computed strength identical to the experimental one. From this analysis, we retained that $S_0 = 0.2$ and $f_0 = 0.01$ N. The other parameters were kept identical to the ones used in the modeling of the saturated samples. The results of the numerical simulations are presented together with the experimental ones in Figs. 5 and 6. The stress-strain relationship shows a good agreement between experimental and numerical results for the three confining pressures. In particular, the values of the maximum strength at large deformations are computed with good accuracy. The computed volume changes correspond also to a contractant behavior of the three specimens. The comparison with experimental results shows that the amplitudes of volume change given by the model are less dependent on the confining pressure than the measured ones, but experimental and numerical results remain reasonably close.

The final states of the three specimens were plotted in the total stress plane p,q together with the critical state line obtained from the tests on saturated samples (Fig. 7). One can see that the state at the end of each test on unsaturated specimens is located above the critical state line. This demonstrates the increase in strength

<figure><img src="image_7.png" /><figcaption>Fig. 7. Maximum strength in a total stress plane.</figcaption></figure>
due to the capillary forces. It can be noted that the final points are closer to the critical state line for higher applied net mean stress. This is due to the volumetric contraction, which increases with the confining pressure. As a consequence, the degree of saturation increases during shearing, causing a decrease of the capillary forces, which is more pronounced at a higher confining pressure. If we extend the maximum strength locus for the unsaturated samples toward the origin of the axes, we obtain an intercept $q_{unsat} = 10$ kPa for $p = 0$, which reflects the contribution of the capillary forces on the maximum strength. If we transfer these results to a Mohr-Coulomb plane, we obtain an apparent cohesion for the unsaturated samples $c_{ap} = 5.3$ kPa and a friction angle similar to the one obtained for saturated samples: $\phi = 37.7^\circ$.

Using Eq. (23), we can calculate the capillary stress tensor and its evolution during triaxial testing. Before shearing, the material is assumed to have an isotropic structure and therefore, the capillary stress tensor is also isotropic. The shearing phase at constant water content induces a structure anisotropy as well as a change in the degree of saturation, leading to an evolution of the capillary stress tensor, which is no longer isotropic. We define a generalized effective stress tensor $\sigma^*$ by adding the applied and the capillary stresses:

$$ \sigma^* = \sigma + \sigma^{\text{cap}}. \qquad (33) $$

We can now plot the final states of the unsaturated samples in a generalized effective stress plane $p^*, q^*$. The results are presented in Fig. 8. One can see that, in this plane, the critical state is unique for saturated and unsaturated samples.

There is no direct relationship between the measured suction and the capillary stress tensor, but their evolution during the shearing phase is similar and controlled by the evolution of the degree of saturation due to the volume change. Figs. 9 and 10 show the evolution of, respectively, the suction and the mean capillary stress, with the axial strain for the three unsaturated specimens. One can see a very similar evolution, even if the absolute values are different.

## 4.2. Undrained triaxial tests. Influence of the degree of saturation

Hadiwardoyo (2002) reported results on an Indonesian sand, a material classified as a silty sand with spread distribution (see grain size distribution curve in Fig. 2). The samples were all prepared following the same procedure: the soil was mixed with a given quantity of water and then compacted to a chosen water content and density. Undrained triaxial tests were performed on saturated samples. Figs. 11 and 12 present the results obtained on three saturated specimens. After compaction, the samples were saturated and then compressed isotropically to three different effective pressures: 100 kPa, 200 kPa and 250 kPa. The initial void ratios before shearing were respectively: 0.59, 0.545 and 0.535. For each specimen, one could see an initial tendency to contract, which led to a decrease of the mean effective stress $p'$, followed by a tendency to dilate, leading to an increase of $p'$. The deviatoric stress increased continuously in the three tests and reached a constant value for axial strains of about 10%.

<figure><img src="image_8.png" /><figcaption>Fig. 8. Maximum strength in a generalized effective stress plane.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 9. Relationship between change of suction and axial strain during constant water content triaxial tests.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10. Relationship between mean capillary stress and axial strain during constant water content triaxial tests.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11. Undrained triaxial tests on silty sand. Deviatoric stress versus axial strain.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 12. Undrained triaxial tests on silty sand. Effective stress paths.</figcaption></figure>

The mean particle size $2R$ was selected equal to $d_{50} = 0.5$ mm. The elastic parameters were not directly accessible from these test results. As for the previous example, typical values were adopted in this study. The inter-particle elastic constant $k_{\mu 0}$ is assumed to be equal to 60000 N/mm and $k_{r0}/k_{n0} = 0.4$. The normalized contact number per volume $Nl^3/V$ can be obtained by using an empirical relation with the void ratio. The test results presented in Figs. 11 and 12 on saturated samples at different confining pressures allowed the determination of the critical state parameters: $M = 1.75$ ($\phi_{\mu} = 42.6^{\circ}$), $\lambda = 0.057$ and $e_{\text{ref}} = 0.64$ for $p_{\text{ref}} = 100$ kPa. From the evolution of the pore water pressure, we were able to determine the slope of the characteristic line, separating the contractive and the dilative domains, which led to a value of $\phi_0 = 40^{\circ}$, slightly different of $\phi_{\mu}$. Table 2 summarizes the set of parameters used for the modeling of the saturated sample behavior. Comparison between experimental and numerical results are presented in Figs. 11 and 12. A very good agreement between experimental and numerical results could be obtained, for the stress–strain relationship as well as for the effective stress paths.

Figs. 13 and 14 present results of constant water content triaxial tests on unsaturated samples. All the specimens were subjected to the same net confining stress $p_0 = 78$ kPa. The initial void ratio varied between 0.48 and 0.54 and the water content varied between 10.5% and 17.8%. The initial degree of saturation took place between 54% and 98%. For high degrees of saturation, the deviatoric stress increased with the axial strain up to a constant value reached for axial strains around 5% and then remained approximately constant. For low degrees of saturation, a peak was observed in the stress–strain curves at axial strains occurring between 4% and 6% (Fig. 14).

This peak could be the consequence of the dilatant behavior of the samples, which can be seen on the volume change evolution in Fig. 15. According to the author, however, strain localization was the mode of failure in most of tests, which could induce a strain softening after the peak. The author also questioned the reliability of the volume change measurements, since they were made through the volume change of water surrounding the specimen in the triaxial cell. Whenever we examine carefully the different curves, we cannot determine a clear relationship between volume change amplitude, degree of saturation or void ratio. Therefore, we gave more importance to the stress–strain relationship, trying only to reproduce on the whole a dilatant behavior, without focusing on the volume change amplitudes.

In order to determine the two parameters $S_0$ and $f_0$ controlling the amplitude of the capillary forces, we adopted at first the same procedure as in the previous section. Eq. (6) led to a value of $S_0 = 0.26$. We tried afterwards to determine a value of $f_0$, which could lead to a good agreement between experimental and

<table><caption>Table 2
Model parameters for the Indonesian silty sand</caption><thead><tr><td>$e_{\text{ref}}$</td><td>$p_{\text{ref}}$ (kPa)</td><td>$\lambda$</td><td>$\phi_{\mu}$ (°)</td><td>$\phi_0$ (°)</td><td>$m$</td></tr></thead><tbody><tr><td>0.64</td><td>100</td><td>0.057</td><td>42.6</td><td>40</td><td>1.5</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Fig. 13. Constant water content triaxial tests on silty sand: stress–strain relationship.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 14. Constant water content triaxial tests on silty sand: volumetric-axial strain relationship.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 15. Computed constant water content triaxial tests on silty sand: stress–strain relationship.</figcaption></figure>
computed results. This goal was not possible to achieve, due in particular to the dilatant behavior of the sam-
ples, which produced a decrease of the degree of saturation during shearing. As a consequence, the capillary
forces inside the specimen increased and created a too strong increase of the material strength. Therefore, we
disregarded the validity of Eq. (6) and tried to determine at the same time the values of the two parameters S₀
and f₀, using the series of test results on unsaturated samples. The following values were finally retained:
S₀ = 0.01 and f₀ = 1 N. These values are very different from the ones determined for the granitic soil, when
the two silty sands have similar grain size distributions (even if the Indonesian silty sand has a bigger percent-
age of fines). Tests on unsaturated granitic soil were performed at the same degree of saturation S_r = 43%.
Thus, this case did not allow us to study the evolution of the capillary forces with the degree of saturation.
Different couples S₀, f₀ can lead to the same amplitude of the capillary forces for a given degree of saturation.
As an example, if we consider S₀ = 0.01 for the granitic soil, its leads to a value f₀ = 0.25 N for equal capillary
force amplitude at S_r = 43%, and therefore to similar computed stress-strain curves for this material. In theses
conditions, there is a good coherence between the parameter values obtained for the two soils, since the higher
amount of fine particles in the case of the Indonesian sand should give higher capillary stress amplitude, i.e.
higher value for f₀ if S₀ has the same value for both soils. The validity of the empirical equation proposed by
Wu et al. was already discussed by Cho and Santamaria (2001), whose results did not agree with the equation.
More studies are certainly required on this topic.

The results from the numerical simulations, using the set of parameters determined from the tests on saturated samples and the values of S₀ and f₀ determined above, are presented in Figs. 15 and 16. An overall good agreement could be achieved, in particular the evolution of the maximum strength as a function of the degree of saturation was well obtained by the model. It should be noted that the experimental results obtained on the specimen corresponding to w = 15.4% and e = 0.48 did not agree with the one corresponding to w = 15% and e = 0.49, whose initial condition was almost identical and which gave a different stress-strain relationship. In the simulations, we replaced these two specimens by another one, having the following initial state: w = 15% and e = 0.5, which allowed us to illustrate in a better way the overall influence of the degree of saturation comprised in the interval 0.5–1. The computed maximum strength was located midway between the values obtained experimentally on the two specimens. The computed volume change shows a dilatant behavior of the unsaturated samples, in agreement with the experimental results. The volumetric strain amplitudes depend both on the degree of saturation and on the initial void ratio. Considering the remarks made by the author on the reliability of the measures, they appear to be reasonably well predicted by the model.

4.3. *Behavior during a wetting phase*

Ng and Chiu (2003) performed also several tests, consisting of wetting specimens of the same silty sand at a
constant stress state. After an isotropic consolidation up to 100 kPa, a specimen was sheared at constant net

<figure><img src="image_16.png" /><figcaption>Fig. 16. Computed constant water content triaxial tests on silty sand: volumetric-axial strain relationship.</figcaption></figure>
mean stress and constant suction, $s = 40$ kPa, until it reached a net stress ratio $\eta = q/p = 1.4$. Then, the specimen was wetted by reducing the suction while the deviatoric stress and the net mean stress were kept constant. Fig. 17 shows the evolution of the axial strain and of the volumetric strain during the wetting phase as a function of the suction. A continuous increase of the axial strain can be noted during the whole process. At the beginning of the wetting, this increase remains limited, toward the end of the test, the rate of increase accelerates. The same kind of evolution can be observed for the volumetric strain. At first, a very limited volumetric contraction takes place, but, as the suction continues to decrease, toward the end of the test the rate of increase in volumetric contraction also increases. The authors attribute this increase in the rate of volumetric strain to the collapse of the soil skeleton.

A numerical simulation of the same wetting test was undertaken. The shearing phase was conducted at a constant capillary stress in order to reproduce, as closely as possible, the condition of constant suction. The numerical results are presented in Fig. 18. An increase of both the axial and the volumetric strain was obtained, as in the experiment, with a continuous increase in the strain rate when the capillary stress decreased. However, the evolution was more progressive than in the experimental results, the specimen starting to strain as early as the beginning of the wetting phase. This could be due to the choice made for the relationship between capillary forces and degree of saturation (Eq. (5)) as well as the way the parameters $S_0$ and $f_0$ were determined, as already discussed in the previous section. The magnitude of the induced volumetric strain remained smaller than the one obtained during the constant water content test, for identical values of axial

<figure><img src="image_17.png" /><figcaption>Fig. 17. Relationship between axial strain, volumetric strain and suction during a wetting test ($p_0$ = 100 kPa).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 18. Computed axial strain and volumetric strain evolution during a wetting phase ($p_0$ = 100 kPa).</figcaption></figure>
strains. This could be explained by the high stress level under which the wetting test was conducted, which reduces the amount of contraction that could be developed by the sample (Eq. (12)). The model is therefore capable of reproducing the phenomenon called the collapse of unsaturated soils during a wetting phase under a deviatoric stress. During the wetting phase, the reduction of the capillary forces with the increase of the degree of saturation weakens the mechanical strength of the unsaturated specimen, which, at a certain point, cannot sustain any longer the applied deviatoric stress.

**5. Summary and conclusion**

A microstructural approach is conceptually appealing, because it dispenses in having to consider relations between the stress and strain tensors. Rather, we need to only describe a simple relation between the vectors of forces and relative displacements on a contact plane, which thus requires fewer material parameters. The introduction of the capillary forces between neighboring grains, function of the degree of saturation of the assembly, allows the model to integrate non saturated as well as dry or fully saturated materials. The stress and strain tensors are obtained by integrating over all spatial orientations.

In the model, a simple elastic-plastic behavior was assumed for each contact plane. The elastic part is based on the Hertz–Mindlin contact formulation, while the plastic part is based on a Mohr–Coulomb friction law with an isotropic hardening assumption and a non associated flow rule. For the entire packing, a critical state behavior is assumed at large deformations and the friction angle on each plane is related to the actual void ratio, compared to the critical void ratio at the same state of stress. On the whole, the model requires only a limited number of parameters, which can easily be determined from conventional triaxial testing.

The ability of the model to reproduce the main features of an unsaturated granular material behavior has been demonstrated. The numerical simulations included constant water content triaxial tests at a given degree of saturation and various net confining pressure, constant water content triaxial tests at a given net confining pressure and various degrees of saturation, and a wetting test at constant deviatoric stress and constant net mean stress. For constant water content triaxial tests, the results showed an increase in the mechanical properties whenever the degree of saturation decreased. The results also demonstrate that the collapse of an unsaturated specimen could be obtained during wetting, provided that a sufficiently high stress ratio was applied during the whole process.

The integration of the capillary forces over all spatial orientations permits to define a capillary stress tensor. This tensor is isotropic, i.e., can be reduced to a capillary pressure, only if the particulate structure is itself isotropic. A generalized effective stress tensor can therefore be determined by adding the total stress and the capillary stress tensors. Numerical simulations demonstrated that the maximum strength criteria, expressed in terms of generalized effective stresses, is unique for dry, partially and fully saturated samples.

The capillary forces were computed as a function of the pore morphology and of the water content, using a simple empirical relationship. Further improvement of the model will require a better understanding of the capillary force amplitude and evolution, based more closely on the physics of the phenomenon. In particular, the transition between saturated and unsaturated states, at high degrees of saturation, as well as the transition between unsaturated and dry states, at low degrees of saturation, need to be better investigated.

**References**

Alonso, E.E., Gens, A., Hight, D.W., 1990. A constitutive model for partially saturated soils. Geotechnique 40 (3), 405–430.

Batdorf, S.B., Budianski, B., 1949. A mathematical theory of plasticity based on concept of slip. NACA Tech Note TN, 1871.

Bazant, Z.P., Xiang, Y., Ozbolt, J., 1995. Nonlocal microplane model for damage due to cracking. Proceedings of Engineering Mechanics 2, 694–697.

Biarez, J., Fleureau, J.-M., Taibi, S., 1993. Mechanical constitutive model for unsaturated granular media. Second International Conferences on Micromechanics of Granular Media. Balkema, Thorton, pp. 51–58.

Biarez, J., Hicher, P.-Y., 1994. Elementary Mechanics of Soil Behaviour. Balkema, p. 208.

Bishop, A.W., Blight, G.E., 1963. some aspects of effective stress in saturated and partly saturated soils. Geotechnique 13 (3), 177–197.

Cambou, B., Dedecker, F., Chaze, M., 2000. Relevant local variables for the change of scale in granular materials. In: Dimitrios, Kolymbas (Ed.), Constitutive Modelling of Granular Materials. Springer, Berlin, pp. 275–290.

Chang, C.S., Gao, J., 1996. Kinematic and static hypotheses for constitutive modeling of granulates considering particle rotation. Acta Mechanica 115 (1–4), 213–229.
Chang, C.S., Hicher, P.-Y., 2005. An elastoplastic model for granular materials with microstructural consideration. International Journal of Solids and Structures 42 (14), 4258-4277.

Chang, C.S., Misra, A., 1990. Application of uniform strain theory to heterogeneous granular solids. Journal of Engineering Mechanics, ASCE. 116 (10), 2310-2328.

Chang, C.S., Misra, A., Weeraratne, S.P., 1989a. A slip mechanism based constitutive model for granular soils. Journal of Engineering Mechanics, ASCE 115 (4), 790-807.

Chang, C.S., Sundaram, S.S., Misra, A., 1989b. initial moduli of particulate mass with frictional contacts. In: International Journal for Numerical and Analytical Methods in Geomechanics. John Wiley and Sons, vol. 13(6), pp. 626-641.

Chang, C.S., Kabir, M., Chang, Y., 1992a. Micromechanics modelling for the stress strain behavior of granular soil-II: evaluation. Journal of Geotechnical Engineering, ASCE 118 (12), 1975-1994.

Chang, C.S., Misra, A., Acheampon, K., 1992b. Elastoplastic deformation of granulates with frictional contacts. Journal of Engineering Mechanics, ASCE 118 (8), 1692-1708.

Chang, C.S., Chao, S.C., Chang, Y., 1995. Estimates of mechanical properties of granulates with anisotropic random packing structure. International Journal of Solids and Structures 32 (14), 1989-2008.

Chateau, X., Moucheront, P., Pitois, O., 2002. Micromechanics of Unsaturated Granular Media. Journal of Engineering Mechanics, ASCE 128 (8), 856-863.

Cho, G.C., Santamaria, J.C., 2001. Unsaturated particulate materials – Particle-level studies. Journal of Geotechnical Engineering, ASCE 127 (1), 84-96.

Christofferson, J., Mehrabadi, M.M., Nemat-Nassar, S., 1981. A micromechanical description on granular material behavior. ASME Journal of Applied Mechanics 48, 339-344.

Cui, Y.J., Delage, P., 1996. Yielding and plastic behaviour of an unsaturated compacted silt. Geotechnique 46 (2), 291-311.

de Buhan, P., Dormieux, L., 1996. On the validity of the effective stress concept for assessing the strength of saturated porous materials: a homogenization approach. Journal of the Mechanics and Physics of Solids 44 (10), 1649-1667.

Goddard, J.D., Bashir, Y.M., 1990. On Reynolds dilatancy. In: De Kee, D., Kaloni, P.N. (Eds.), Recent Development in Structured Continua, vol. 2. Longman's, London, pp. 23-35.

Hadiwardoyo, S., 2002. Caractérisation des matériaux granulaires. Thèse de doctorat de l'Ecole Centrale de Paris.

Hicher, P.-Y., 1996. Elastic properties of soils. Journal of Geotechnical Engineering, ASCE 122 (8), 641-648.

Hicher, P.-Y., 1998. Experimental behavior of granular materials. In: Cambou, B. (Ed.), Behavior of Granular Materials. Springer, Wien New York, pp. 1-97.

Hicher, P.-Y., Chang, C.S., 2006. Anisotropic non linear elastic model for particulate materials. Journal of Geotechnical and Environmental Engineering, ASCE 132 (8), 1052-1061.

Jenkins, J.T., Strack, O.D.L., 1993. Mean-field inelastic behavior of random arrays of identical spheres. Mechanics of Material 16, 25-33.

Kohgo, Y., Nakano, M., Miyazaki, T., 1993. theoretical aspects of constitutive modeling for unsaturated soils. Soils and Foundations 33 (4), 49-63.

Kruyt, N.P., 2003. Statics and kinematics of discrete Cosserat-type granular materials. International Journal of Solids and Structures 40, 511-534.

Kruyt, N.P., Rothenburg, L., 2002. Micromechanical bounds for the effective elastic moduli of granular materials. International Journal of Solids and Structures 39 (2), 311-324.

Liao, C.L., Chang, T.P., Young, D., Chang, C.S., 1997. Stress-strain relationship for granular materials bases on hypothesis of best fit. International Journal of Solids and Structures 34 (31-32), 4087-4100.

Matsuoka, H., Takeda, K., 1980. A stress-strain relationship for granular materials derived from microscopic shear mechanisms. Soils and Foundation 20 (3), 45-58.

Mindlin, R.D., 1969. Microstructure in linear elasticity. Archive for Rational Mechanics and Analysis 16, 51-78.

Modaressi, A., Abou-Bekr, N., 1994. constitutive model for unsaturated soils; validation on a silty materialProceedings of Numerical Methods in Geomechanics, NUMOG V, 1. Balkema.

Ng, C., Chiu, A., 2003. Laboratory study of loose saturated and unsaturated decomposed granitic soil. Journal of Geotechnical Engineering, ASCE 129 (6), 550-559.

Pande, G.N., Sharma, K.G., 1982. Multi-laminate model of clays – a numerical evaluation of the influence of rotation of the principal stress axis. In: Desai, C.S., Saxena, S.K. (Eds.), Proceedings of Symposium on Implementation of Computer Procedures and Stress-Strain Laws in Geotechnical Engineering. Acorn Press, Chicago, Durham, NC, pp. 575-590.

Rothenburg, L., Selvadurai, A.P.S., 1981. Micromechanical definitions of the Cauchy stress tensor for particular media. In: Selvadurai, A.P.S. (Ed.), Mechanics of Structured Media. Elsevier, Amsterdam, pp. 469-486.

Rowe, P.W., 1962. The stress-dilatancy relations for static equilibrium of an assembly of particles in contact. Proceedings of the Royal Society London A 269, 500-527.

Suiker, A.S.J., Chang, C.S., 2004. Modelling failure and deformation of an assembly of spheres with frictional contacts. Journal of Engineering Mechanics, ASCE 130 (3), 283-293.

Taibi, S., 1994. Comportement mécanique et hydraulique des sols partiellement saturés. Thèse de doctorat, Ecole Centrale de Paris, p. 374.

Taylor, D.W., 1948. Fundamentals of Soil Mechanics. John Wiley and Sons, New York, NY.

Terzaghi, K., 1925. Principles of soil mechanics, a summary of experimental results of clay and sand. Engineering News Record, 3-98.

Wheeler, S.J., Sivakumar, V., 1995. An elasto-plastic critical state framework for unsaturated soils. Géotechnique 45 (1), 35-53.

Wu, S., Gray, D.H., Richart Jr., F.E., 1984. capillary effects on dynamic modulus of sands and silts. Journal of Geotechnical Engineering ASCE 110 (9), 1188-1203.