# Behavior of unbound granular materials - Part I: isotropic case

Rabah Bouzidi, L. Coulibaly, P. Jouve

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

## To cite this version:

Rabah Bouzidi, L. Coulibaly, P. Jouve. Behavior of unbound granular materials - Part I: isotropic case. Computers and Geotechnics, 2003, Vol 30 (2), pp.185-204. 10.1016/S0266-352X(02)00060-5. hal-01005288

HAL Id: hal-01005288

https://hal.science/hal-01005288

Submitted on 6 Dec 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

<figure><img src="image_2.png" /></figure>
# Behavior of unbound granular materials—
part I: isotropic case

Rabah Bouzidi, L. Coulibaly, P. Jouve

Laboratoire de Génie Civil de Nantes Saint-Nazaire, BP92208 44322 Nantes Cedex, France

The paper discusses the modeling of the behavior of unbound granular materials. A representative approach that highlights some salient features of the behavior is proposed. This approach is essentially based on experimental results and the study is extended to the construction of the elastic potential from test results. to complete the analysis, two no-linear elastic models involving 3 parameters are proposed. In the construction of these models, two important aspects—the accuracy and the numerical stability—are analyzed.

Keywords: Unbound granular materials; Behavior relationship; Elastic potential; Isotropy; Triaxial test

## 1. Introduction

The usual approach for the modeling of the unbound materials behavior is to search for a constitutive relationship that describes as closely as possible the laboratory test results and then to present them within a consistent thermodynamic framework. These laws are often based on observations and measurements performed during the experiments. However, this kind of approach does not answer some fundamental questions such as: does the observed kinematic field derive from a potential function? In the case of obtaining the potential expression, is it the closest one to the experimental observations? Is the behavior isotropic? The following additional questions can arise during the use of these models: is the model still valid
beyond the experimental data range? Can this model class be able to take into account the most important behavior characteristics of the material?

This paper summarizes our modeling methodology of unbound granular materials used in flexible pavements. It is an attempt to answer the above questions without choosing any a prior model. The key idea is to obtain as much information as possible from the laboratory test results. This work is only focused on the reversible aspect of the material behavior, so that no cyclic or plastic aspects are considered.

## 2. Background

Materials used in the sub-layers of pavement must satisfy two important criteria. The first criterion related to the cost, since the material must be inexpensive and available in large quantities. The second criterion is relative to the mechanical response of the material to the stresses induced by the traffic load. The materials obtained from crushed rock present a good quality-price compromise when used in pavement sub-layers. Hereafter, we refer to these materials as the unbound granular materials in contrast to those used in the upper layer like bituminous concrete. This kind of material can include a broad or a limited grading distribution curve, with the particle size ranging from 14 to 31.5 mm with 5–10% of fines (i.e. <0.075 mm) [2]. In the laboratory, we obtain these materials by combining and mixing several crushed rock particle sizes. For a detailed description of the subject, see e.g. Refs. [5,6].

### 2.1. *Experimental study of unbound granular material*

The rolling load located at the abscise $\lambda$ on the road surface induces vertical and horizontal stresses at the point M (Fig. 1). The ratio of the vertical stress to the horizontal stress $\sigma_1/\sigma_3$ depends on the position M. This particular repeated load path, due to the traffic, has been investigated by several laboratories developing specific devices to characterize the behavior of the material for these loadings. In 1980, the Laboratoire Central des Ponts et Chaussées (LCPC) developed specific devices for testing such granular materials. The tests were carried out with variable lateral confining pressures on cylindrical samples of 16 cm in diameter and 32 cm in

<figure><img src="image_3.png" /><figcaption>Fig. 1. Stress path under rolling load.</figcaption></figure>
height [3,4]. In this apparatus, initially developed in 1974 at Nottingham university by the Brown S.F. team., the confining pressure varies so that, the ratio of the deviator stress to the confining pressure remains constant (Fig. 2). The stress paths are set to reproduce the in situ loading conditions.

A cycling preconditioning before the test is performed. Over approximately 20,000 of load-unload cycles, with $σ_1/σ_3$ ratio equal to 6, are applied in order to obtain a quasi-reversible behavior, with very small plastic deformations (Fig. 3). The specimens obtained after this first stage of the tests are compact and exhibit significant cohesion due to the fines and moisture of the material.

After the preconditioning stage, the test is carried out with different values of $σ_1/σ_3$ ratio. The confining pressure of the cell $σ_3$ and the axial stress $σ_1$ are recorded in order to evaluate the mean pressure $p=(σ_1 + 2σ_3)/3$ and the shear stress $q=σ_1-σ_3$. Also, the axial and radial strains are measured to estimate the volumetric strain $ε_v = ε_1 + 2ε_3$ and shear strain $ε_q = 2(ε_1-ε_3)/3$. Fig. 4a and b shows an example of an interpolated test results, representing the volumetric and shear strain in the (p, q) space.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Stress path for triaxial variable confining stress test.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Evolution of the axial distortion according to cycles.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Interpolated results of the triaxial test. (a) Volumetric strains, (b) shear strains.</figcaption></figure>

## 2.2. Modeling of unbound granular materials

The Boyce model is the most used to represent the behavior of such unbound granular material [1]. It is based on the generalized elastic law, in which Young's modulus E and Poisson ratio ν are given as functions of stress invariants. The model is isotropic with nonlinear elasticity and the strain components are obtained by:

$$ \epsilon_{ij} = \frac{1 + \nu}{E} \sigma_{ij} - \frac{3\nu}{E} p \delta_{ij} \quad \text{With:} \quad E = \frac{9KG}{3K + G} \quad \nu = \frac{3K - 2G}{6K + 2G} \quad (1) $$

The bulk modulus K and the shear modulus G are defined by the author as follows:

$$ K = K_a \left(\frac{p}{p_a}\right)^{1-n} \left(1 - \beta\eta^2\right)^{-1} , \quad G = G_a \left(\frac{p}{p_a}\right)^{1-n} \quad (2) $$

with $0 < n < 1$ and $\eta = p/q$. Thus, Young's modulus and Poisson ratio are obtained by:

$$ E = \frac{9G_a \left(\frac{p}{p_a}\right)^{1-n}}{3 + \frac{G_a}{K_a} \left[1 - \beta \eta^2\right]} \quad \nu = \frac{\frac{3}{2} - \frac{G_a}{K_a} \left[1 - \beta \eta^2\right]}{3 + \frac{G_a}{K_a} \left[1 - \beta \eta^2\right]} \quad (3) $$

From Eq. (1), we can write the volumetric and shear strains in the form:

$$ \epsilon_v = \frac{p_a}{K_a} \left(\frac{p}{p_a}\right)^n (1 - \beta \eta^2) , \quad \epsilon_q = \frac{p_a}{3G_a} \left(\frac{p}{p_a}\right)^n \eta \quad (4) $$

where $p_a$ is a reference pressure taken equal to 100 kPa. The model includes 3 parameters $K_a$, $G_a$ and $n$ together with a fourth parameter, $\beta$, which depends on the
<figure><img src="image_7.png" /><figcaption>Fig. 5. Interpolated results of the Boyce's model. (a) Volumetric strains, (b) shear strains.</figcaption></figure>

three others, so that the strains derive from an elastic potential. The condition of the existence of the elastic potential is given by the following equation:

$$ \frac{\partial \varepsilon_v}{\partial q} = \frac{\partial \varepsilon_q}{\partial p} \quad \Rightarrow \quad \beta = (1 - n) \frac{K_a}{6G_a} \qquad (5) $$

Boyce wrote the complementary elastic potential as follows:

$$ U^c = \frac{p_a}{(1+n)K_a} \left(\frac{p}{p_a}\right)^{n+1} \left[1 + \frac{(1+n)K_a}{6G_a}\eta^2\right] \qquad (6) $$

The strains obtained by the model are shown in Fig. 5 with the following fitted parameters: $K_a = 187$ MPa, $G_a = 298$ MPa and $n = 0.32$. These results, compared with those obtained from the experiments (Fig. 4), show that the model describes the material with a very satisfactory accuracy. However, its application beyond the experimental data field can lead to unrealistic results. The model shows singular expressions of Young's modulus and Poisson ratio for $\eta^2 = (1 + 3K_a/G_a)/\beta$, [Eq. (3)]. These non-suitable singularities can occur in structure computations when the model is used without plasticity bounding surface in the stresses space. In the following section, a numerical procedure, based on test results, is proposed to characterize the main properties of the unbound granular materials. No theoretical model is considered at this stage of the procedure.

### 3. Analysis of the reversible behavior

A new analysis procedure of the test results is presented in this section in order to obtain as much information as possible from the experiment. Regardless of the
type of model to be used later, this procedure provides enough details on the following:

* the existence of an elastic potential,

* the general form of this potential. Especially, the assessment of the influence of the mean pressure p and the shear ratio η,

* the validity of an isotropic behavior assumption.

### 3.1. Existence of the potential

Strains tensor is obtained from the complementary elastic potential as:

$$ \epsilon_{ij} = \frac{\partial U^c}{\partial \sigma_{ij}} \qquad (7) $$

In the case of axisymmetrical loads, the volumetric and shear strains are given by:

$$ \epsilon_v = \frac{\partial U^c}{\partial p} \qquad \epsilon_q = \frac{\partial U^c}{\partial q} \qquad (8) $$

In order to check the two relations given above, the values of the function $U^c(p, q)$ will be determined using strain measurements recorded during a triaxial compression test. Given the values of $(p, q)$ used during the tests, it is assumed that the complementary elastic potential is a function of the independent variables $p$ and $\eta$. Relation (8) then becomes:

$$ \epsilon_v = \frac{\partial}{\partial p} U^c(p, \eta) - \frac{\eta}{p} \frac{\partial}{\partial \eta} U^c(p, \eta) \qquad \epsilon_q = \frac{1}{p} \frac{\partial}{\partial \eta} U^c(p, \eta) \qquad (9) $$

Hereafter, we first present a methodology for relating $\epsilon_v$ and $\epsilon_q$ to grid nodal values of the potential $U_c$ in $(p, \eta)$ space. Next, we proceed to the calculation of the nodal values of $U_c$ by identifying the measured values of the strains: $\epsilon_v^m, \epsilon_q^m$ and those calculated from the potential: $\epsilon_v^c$ and $\epsilon_q^c$

#### 3.1.1. Local potential interpolation

In the square $abcd$ in Fig. 6, the slope of the complementary elastic potential $U^c(p, \eta)$ is assumed to be a bilinear function with respect to the variables $p$ and $\eta$.

At any point in local co-ordinates $(x, y)$, this slope can be written as:

$$ U_{,x}(x, y) = \alpha_1 + \beta_1 \frac{x}{h} + \gamma_1 \frac{y}{k} + \delta_1 \frac{xy}{hk} \qquad (10) $$

This expression must satisfy the slope calculated inside the square $a, b, c$ and $d$.
$$U_{,x}(a) = \alpha_1 - \frac{\beta_1}{2} = \frac{U_2 - U_1}{h}$$

$$U_{,x}(b) = \alpha_1 + \frac{\beta_1}{2} = \frac{U_3 - U_2}{h}$$

$$U_{,x}(c) = \alpha_1 - \frac{\beta_1}{2} + \gamma_1 - \frac{\delta_1}{2} = \frac{U_5 - U_4}{h}$$

$$U_{,x}(d) = \alpha_1 + \frac{\beta_1}{2} + \gamma_1 + \frac{\delta_1}{2} = \frac{U_6 - U_5}{h}$$

The equation system obtained gives the interpolation factors. By substituting the expressions of these parameters in Eq. (10), it follows that:

$$U_{,p}^{c} = \sum_{i=1}^{6} a_i U_i \qquad (11)$$

with

$$a_1 = \frac{1}{h} \left( -\frac{1}{2} + \frac{x}{h} + \frac{y}{2k} - \frac{xy}{hk} \right)$$

$$a_2 = \frac{1}{h} \left( -\frac{2x}{h} + \frac{2xy}{hk} \right)$$

$$a_3 = \frac{1}{h} \left( \frac{1}{2} + \frac{x}{h} - \frac{y}{2k} - \frac{xy}{hk} \right)$$

$$a_4 = \frac{1}{h} \left( -\frac{1}{2} \frac{y}{k} - \frac{xy}{hk} \right)$$

$$a_5 = \frac{1}{h} \left( -\frac{2xy}{hk} \right)$$

$$a_6 = \frac{1}{h} \left( \frac{1}{2} \frac{y}{k} + \frac{xy}{hk} \right)$$

Relation (11) gives the slope $U_{,p}^c(x,y)$ at any point in the square $abcd$ according to its co-ordinates in the local reference $(o, x, y)$ and the values of the complementary elastic potential at the six nodes. The same previous methodology could be adopted to establish the expression of $U_{,\eta}^c$. By using similar interpolation inside square $efgh$ in Fig. 6b it can be shown that $U_{,\eta}^c$ is written, in the local reference $(o, x, y)$, as follows:

$$U_{,\eta}^{c} = \sum_{i=1}^{6} b_i U_i \qquad (12)$$
<figure><img src="image_8.png" /><figcaption>Fig. 6. Interpolation fields of $U_{,p}$ and of $U_{,\eta}$.</figcaption></figure>

with

$$b_1 = \frac{1}{k} \left( -\frac{1}{2} + \frac{x}{2h} + \frac{y}{k} - \frac{xy}{hk} \right)$$

$$b_2 = \frac{1}{k} \left( -\frac{x}{2h} + \frac{xy}{hk} \right)$$

$$b_3 = \frac{1}{k} \left( -\frac{2y}{k} + 2\frac{xy}{hk} \right)$$

$$b_4 = \frac{1}{k} \left( -2\frac{xy}{hk} \right)$$

$$b_5 = \frac{1}{k} \left( \frac{1}{2} - \frac{x}{2h} + \frac{y}{k} - \frac{xy}{hk} \right)$$

$$b_6 = \frac{1}{k} \left( \frac{x}{2h} + \frac{xy}{hk} \right)$$

### 3.1.2. Evaluation of $U_{,p}^c$ and $U_{,\eta}^c$ from measurements

During triaxial compression tests with variable confining pressure, a set of $m$ measurement points $(p, q, \varepsilon_v, \varepsilon_q)$ is recorded. Given a range of experimental data $p$ and $\eta$, a regular rectangular grid with $n$ nodes can be formed as we can see in Fig. 7.

The triangular dots refer to the measurement points obtained during the tests. From the interpolation systems [Eqs. (11) and (12)], a system of $2m$ equations with $n$ unknowns, is obtained for the whole measurement points.

$$\{U_{,p}^c\} = [A]\{U\} \quad \{U_{,\eta}^c\} = [B]\{U\}$$
The elements of the matrices [A] and [B] depend on the position of the measure-
ment points (p, η). From Eq. (9), the volumetric and shear strains at the measure-
ment points, can be obtained from the elastic potential:

$$
\{\epsilon_v^c\} = \{U_{,p}^c\} - \left\{\frac{\eta}{p} U_{,\eta}^c\right\} \quad \{\epsilon_q^c\} = \left\{\frac{1}{p} U_{,\eta}^c\right\}
$$

it follows in matrix notations:

$$
\{\epsilon_v^c\} = [M_v]\{U\} \quad \{\epsilon_q^c\} = [M_q]\{U\} \tag{13}
$$

$$
\text{with } (M_v)_{ij} = A_{ij} - \frac{\eta_i}{p_i} B_{ij} (M_q)_{ij} = \frac{1}{p_i} B_{ij}
$$

If the condition 2*m* > *n* is satisfied, an approximate solution can be obtained for the equation system (13). In practice, this condition should be taken into account when forming the grid of interpolation. The solution of this problem is given by optimizing the potential using the least squares method. The quantity to be minimized is the square difference between the measured strains ε^m and the calculated strains ε^c from the potential [Eq. (13)].

$$
\Phi = \sum_{i=1}^{m} \left[ \frac{1}{\vartheta_v} (\varepsilon_{vi}^{c} - \varepsilon_{vi}^{m})^2 + \frac{1}{\vartheta_q} (\varepsilon_{qi}^{c} - \varepsilon_{qi}^{m})^2 \right]
$$

The variances of the measurements $\vartheta_v$ and $\vartheta_q$ are defined as follows:

$$
\vartheta_v = \sum_{i=1}^{m} (\varepsilon_{vi}^{m} - \bar{\varepsilon}_v) \quad \vartheta_q = \sum_{i=1}^{m} (\varepsilon_{qi}^{m} - \bar{\varepsilon}_q)
$$

The variances are used here to give comparable weights to the volumetric and
shear strains. Symbols $\bar{\epsilon}_v$ and $\bar{\epsilon}_q$ denote the averages of the measurements. The
minimization is achieved by writing the optimal condition $\partial\Phi/\partial U_i$ at each node of
the grid which results in a system of $n$ equations with $n$ unknowns $\{U\}$.

$$
[K]\{U\} = \{F\} \tag{14}
$$

with

$$
[K] = \frac{1}{\vartheta_v} [M_v]^T [M_v] + \frac{1}{\vartheta_q} [M_q]^T [M_q]
$$

and

$$
\{F\} = \frac{1}{\vartheta_v} [M_v]^T \{\epsilon_v^m\} + \frac{1}{\vartheta_q} [M_q]^T \{\epsilon_q^m\}
$$
<figure><img src="image_9.png" /><figcaption>Fig. 7. Example of the potential grid and the positions of the measurement points.</figcaption></figure>

Since the matrix [K] is singular, the condition $U_1=0$ must be added in order to obtain a unique solution. The quality of the adjustment of the two functions is evaluated by a correlation coefficient $\rho$ defined as follows:

$$ \ln(\rho) = -\sqrt{\frac{1}{2}\left[\frac{1}{\vartheta_v}(\varepsilon_{vi}^c - \varepsilon_{vi}^m)^2 + \frac{1}{\vartheta_q}(\varepsilon_{qi}^c - \varepsilon_{qi}^m)^2\right]} \quad (15) $$

The coefficient $\rho$ varies from zero to one. It tends toward one when the difference between the calculated value and the measured one vanishes.

### 3.2. Application to the test results

The method for calculating the potential values at the grid nodes [Eq. (14)], is applied to laboratory tests. The materials are reconstituted gravels of different kinds [4], typically Ecuelles and Poulmar'H ones which had been combined in four different way. From 10 up to 25 measurement points ($p_i, q_i, \varepsilon_{vi}, \varepsilon_{qi}$) are obtained from the tests. The efficiency of the method proposed here is tested on hypothetic "experimental" results obtained using Boyce's law. These examples were based on real values of the loadings ($p, q$). The theoretical strains are calculated with Boyce's model and we obtain a table of values ($p_i, q_i, \varepsilon_{vi}, \varepsilon_{qi}$) similar to that given by the experimental tests. The correlation coefficients obtained for these tests are greater than 96%. From a theoretical point of view, the correlation must be equal to 100%. Nevertheless, the results of the calculations indicate an acceptable correlation. The marginal inaccuracy mentioned above is mainly due to the approximation of the derivative evaluation of the potential [Eqs. (11) and (12)]. However, the nodal
approximation of the potential can be considered acceptable. In the next step, this method is applied to experimental results.

Table 1 shows that the values of the correlation coefficient are very satisfactory. The shape of the surface of the complementary elastic potential is similar to that presented in Fig. 8. The correlation coefficients given by the adjustment calculations are good since the values are higher than 79%. Therefore, we conclude that the behavior of the materials considered here can be modeled by laws—to be defined—which derive from elastic potential. Moreover, it should be noticed from Fig. 8, that when the variable *p* or *η* is constant, the curves obtained on the potential surface seem to be homothetical and regular. Thus, we can assume the existence, for each material, of two main curves, one depending only on *p* and the other only on *η*. In other words, the elastic potential can be expressed by the product of two independent functions given as follows: $U^c(p, \eta) = f(p) g(\eta)$.

### 3.3. Variables separation

The purpose of this section is to check whether the complementary potential associated with the measurements can be written as a product of two independent functions $f(p)$ and $g(\eta)$. In the grid considered in Fig. 8, there are *k* values of $U^c$ corresponding to *k* values of *p* by line of *η* constant and there are *r* values of $U^c$ corresponding to *r* values of *η* by column of *p* constant. At each node of the grid, the complementary potential is obtained by the product of $f_i(p)$ by $g_j(\eta)$, with $i=1..k$ and $j=1..r$. In addition, $f_i$ is assumed equal to zero, so that a single decomposition can be performed. The functions $f_i$ and $g_j$ are given using the nonlinear method of the least squares which consists in minimizing the following quantity:

$$ \Psi = \sum_{i,j} (f_i g_j - U_{ij}^c)^2 $$

The minimization of $\Psi$ is obtained by solving a nonlinear system formed by the following equations:

$$ \frac{\partial \Psi}{\partial f_i} = 0 \quad \Rightarrow \quad f_i = \frac{\sum_{j=1}^{r} g_j U_{ij}^c}{\sum_{j=1}^{r} (g_j)^2} \quad i = 1..k $$

$$ \frac{\partial \Psi}{\partial g_j} = 0 \quad \Rightarrow \quad g_j = \frac{\sum_{i=1}^{k} g_j U_{ij}^c}{\sum_{i=1}^{k} (g_j)^2} \quad j = 1..r $$

An iterative process gives the solution of this system. The value $f_i$ are computed

<table><caption>Table 1
Correlation coefficients of the elastic potential surface</caption><thead><tr><th>Material</th><th>Ecu1</th><th>Ecu2</th><th>Ecu3</th><th>Ecu4</th><th>Poul1</th><th>Poul2</th><th>Poul3</th><th>Poul4</th></tr></thead><tbody><tr><th>ρ (%)</th><td>83</td><td>93</td><td>79</td><td>90</td><td>91</td><td>90</td><td>90</td><td>81</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Complementary elastic potential obtained for material Ecuelles2.</figcaption></figure>

using those of $g_j$, and $g_j$ values are obtained from those of $f_i$. The initial values of $g_j$ are taken equal to 1. In order to assess the quality of the calculation process given above, a use is made of the correlation coefficient $\rho$ defined by:

$$ \ln(\rho) = -\sqrt{\frac{1}{2} \left( \frac{1}{\vartheta_u} \sum_{i,j} (f_i g_j - U_{ij}^c)^2 \right)} $$

with the variance $\Sigma_{i,j}(U_{ij}^c - \bar{U}^c)^2$ and the average $\bar{U}^c = \frac{1}{r+k} \sum U_{ij}^c$

### 3.4. Application to the test results

The application of the decomposition method of the potential to two independent functions $f(p)$ and $g(\eta)$, gives very satisfactory results. Table 2 summarizes the values of the correlation coefficients obtained for each tested material.

These results validate the assumption of the separation of the variables $p$ and $\eta$. Fig. 9a and b depicts functions $f(p)$ and $g(\eta)$ for the three materials.

The function $f(p)$ always shows the same monotonous and increasing shape. On the other hand, the function $g(\eta)$ generally presents a parabolic shape but shows

<table><caption>Table 2
Correlation coefficients of the adjustments of $f(p)$ and $g(\eta)$</caption><thead><tr><th>Material</th><th>Ecu1</th><th>Ecu2</th><th>Ecu3</th><th>Ecu4</th><th>Poul1</th><th>Poul2</th><th>Poul3</th><th>Poul4</th></tr></thead><tbody><tr><td>ρ (%)</td><td>98</td><td>96</td><td>95</td><td>98</td><td>95</td><td>98</td><td>98</td><td>97</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Curve of f(p) and g(η) for the different materials.</figcaption></figure>

some differences on the original slope according to the type of material. The effects of the average stress p and the shearing stress η can be separated in the complementary elastic potential expression. The expressions of functions f(p) and g(η) will be determined in the next section.

## 4. A generic expression for the elastic potential

### 4.1. Isotropy condition

When a material behavior is assumed isotropic and elastic linear or non-linear, the volumetric and shear strains can be obtained using Hooke's law:

$$ \epsilon_v = \frac{p}{K} \quad \epsilon_q = \frac{q}{3G} $$

the behavior derives from an elastic potential, the strains can be written as:

$$ \begin{gathered} \epsilon_v(p, \eta) = \frac{\partial U^c(p, \eta)}{\partial p} - \frac{\eta}{p} \frac{\partial U^c(p, \eta)}{\partial \eta} \qquad \epsilon_q(p, \eta) = \frac{1}{p} \frac{\partial U^c(p, \eta)}{\partial \eta} \\ \epsilon_q(p, 0) = \frac{1}{p} \frac{\partial U^c(p, 0)}{\partial \eta} = 0 \end{gathered} $$

This equation is a condition of isotropy which means that when a material is defined by an isotropic elastic behavior, the corresponding function g(η) must have a zero slope at the origin on Fig. 9b. Some materials considered here present non-zero slopes at the origin which denote orthotropic behavior at the hydrostatic stress state. For a good interpretation of the results, this fact should be considered in the analysis of this type of materials behavior.
## 4.2. Expression of the elastic potential

The choice of the expressions of the two functions $f(p)$ and $g(\eta)$ must be based on two different criteria. Obviously, the first criterion is related to the representative-ness of the model for the experimental, and the second criterion relates to the use of these laws in the computer programs. For this last case, it is significant to make sure that the model works correctly beyond the experimental data. Among several functions which can be used to represent the functions $f(p)$ and $g(\eta)$; we preferred the following expressions :

$$f(p) = \frac{p_a^{1-n} p^{1+n}}{(1+n)K_a} \quad g(\eta) = \frac{1+\alpha\eta^2}{\sqrt{1+\gamma\eta^2}} \quad (1+n) \quad \text{with } 0 < n < 1$$

These two expressions, which stem from the curves in Fig. 9a and b, define the model entirely. It can be seen that $g(\eta)$ function satisfy the conditions $g(0)=0$ and $g'(0)=0$. Thus, the model is isotropic and needs four parameters $K_a, n, \alpha$ and $\gamma$, whereas $p_a$ is taken equal to 100 kPa. The resulting volumetric and shear strains are given as:

$$\epsilon_v(p, q) = \frac{\partial U^c(p, \eta)}{\partial p} = \frac{p}{K(p, q)} \quad \epsilon_q(p, q) = \frac{\partial U^c(p, q)}{\partial q} = \frac{q}{3G(p, q)}$$

$$K(p, q) = K_a \left(\frac{p}{p_a}\right)^{1-n} \left(\frac{\sqrt{1 + \gamma\eta^2}}{1 + \alpha\eta^2}\right)^{1+n} \frac{(1 + \alpha\eta^2)(1 + \gamma\eta^2)^2}{1 + (2\gamma - \alpha)\eta}$$

$$G(p, q) = \frac{K_a}{3} \left(\frac{p}{p_a}\right)^{1-n} \left(\frac{\sqrt{1 + \gamma\eta^2}}{1 + \alpha\eta^2}\right)^{1+n} \frac{(1 + \alpha\eta^2)(1 + \gamma\eta^2)}{2\alpha - \gamma + \alpha\gamma\eta^2}$$

For large values of $\eta$, the bulk and shear moduli can be written as follows:

$$K(p, q) \approx K_a \left(\frac{p}{p_a}\right)^{1-n} \frac{\gamma^{\frac{n+3}{2}}}{\alpha^n (2\gamma - \alpha)} \eta^{1-n} \quad (16)$$

$$3G(p, q) \approx K_a \left(\frac{p}{p_a}\right)^{1-n} \frac{\gamma^{\frac{n+1}{2}}}{\alpha^{n+1}} \eta^{1-n} \quad (17)$$

and

$$\frac{G}{K} \approx \frac{2\gamma - \alpha}{3\alpha\gamma} \quad (18)$$

It can be verified by using this $G/K$ ratio and from the Eq. (1) that for an infinite value of $\eta$, Poisson ratio and Young's modulus remain finite. This result will be
shown in the next subsection. In the case of low values for $\eta$, the second order Taylor expansion gives the following volumetric and shear strains:

$$ \epsilon_v \approx \frac{p_a^{1-n} p^n}{K_a} (1 - \beta \eta^2) \quad \text{with} \quad \beta = (1 - n) \frac{2\alpha - \gamma}{2} \qquad (19) $$

$$ \epsilon_q \approx \frac{p_a^{1-n} p^n}{K_a} (2\alpha - \gamma)\eta \qquad (20) $$

From Boyce's law, we have $\beta=(1-n)K_a/6G_a$. Replacing this result in Eq. (19), it ensures:

$$ \frac{2\alpha - \gamma}{K_a} = \frac{1}{3G_a} $$

This ensures that our model is still close to Boyce's one for the low values of the shear ratio $\eta$.

### 4.3. Formulation of E and $\nu$

Using Eq. (1), the Young's modulus and the Poisson ratio can be written as:

$$ E = 9K_a \left(\frac{p}{p_a}\right)^{1-n} \left(\frac{\sqrt{1 + \gamma\eta^2}}{1 + \alpha\eta^2}\right)^{1+n} \frac{(1 + \alpha\eta^2)(1 + \gamma\eta^2)}{(9(2\alpha - \gamma) + 1) + (2\gamma - \alpha + 9\alpha\gamma)\eta^2} \qquad (21) $$

$$ \nu = \frac{(9(2\alpha - \gamma) - 2) + (9\alpha\gamma - 2(2\gamma - \alpha))\eta^2}{(18(2\alpha - \gamma) + 2) + (18\alpha\gamma + 2(2\gamma - \alpha))\eta^2} \qquad (22) $$

The model presented above has 4 independent parameters: $K_a$, $\alpha$, $n$ and $\gamma$. The parameter $G_a$ is obtained from Eq. (18) as:

$$ \alpha = \frac{1}{2} \left( \gamma + \frac{K_a}{3G_a} \right) \qquad (23) $$

The presence of a horizontal asymptotic line for both E and $\nu$ [Eq. (21)] indicates that the formulation adopted offers a less brutal variation of these functions for high values of the $q/p$ ratio. Consequently, the model formulated does not present numerical instabilities. However, from the triaxial compression test results, the determination of the parameters becomes difficult or even impossible. It is then necessary to determine the parameters by other experimental results. In the following section, relations between the parameters are proposed to reduce their number.

## 5. Models with three parameters

Because of the limitation described above, we limit the model to 3 parameters so that they can be identified from triaxial test results. We considered two cases with a simple relationship between the parameters $\alpha$ and $\gamma$.
### 5.1. Model with constant Poisson ratio

As a first possibility, let $\gamma = \alpha$. The volumetric and shear strains can be rewritten as follows:

$$ \epsilon_v = \frac{p_a^{1-n} p^n}{K_a} (1 + \alpha \eta^2)^{\frac{n-1}{2}} \quad \epsilon_q = \frac{p_a^{1-n} p^n}{3G_a} (1 + \alpha \eta^2)^{\frac{n-1}{2}} \eta $$

Eq. (23) leads to $\alpha = K_a/(3G_a)$. Eq. (16) gives $G/K = \epsilon_v\eta/(3\epsilon_q = 1/(3\alpha))$, Eq. (3) yields to:

$$ \nu = \frac{9\alpha - 2}{2(9\alpha + 1)} \qquad (24) $$

Obviously, the Poisson ratio is constant. To get a positive Poisson ratio, the condition $\alpha > 2/9$ must be satisfied and therefore we obtain $\nu \le 1/2$, which guarantees the condition of elasticity.

### 5.2. Model with a variable Poisson ratio

As a second possibility, we set $\gamma = \alpha/2$ and therefore we get from Eq. (23), $\alpha = 2K_a/9G_a$. The volumetric and shear strains are given as follows:

$$ \epsilon_v = \frac{p_a^{1-n} p^n}{K_a} \left( \frac{1 + \alpha \eta^2}{\sqrt{1 + \frac{\alpha}{2} \eta^2}} \right)^{1+n} \frac{1}{(1 + \alpha \eta^2)(1 + \frac{\alpha}{2} \eta^2)} $$

$$ \epsilon_v = \frac{p_a^{1-n} p^n}{K_a} \left( \frac{1 + \alpha \eta^2}{\sqrt{1 + \frac{\alpha}{2} \eta^2}} \right)^{1+n} \frac{3 + \alpha \eta^2}{(1 + \alpha \eta^2)(1 + \frac{\alpha}{2} \eta^2)} \eta $$

The ratio $G/K$ becomes:

$$ \frac{G}{K} = \frac{\epsilon_v}{3\epsilon_q} \eta = \frac{2}{3\alpha(3 + \alpha\eta^2)} $$

$$ \nu = \frac{9\alpha^2\eta^2 + 27\alpha - 4}{2(9\alpha^2\eta^2 + 27\alpha + 2)} \qquad (25) $$

In the present case, the Poisson ratio varies with $\eta$ as shown in Fig. 10.

### 5.3. Comparison between the proposed models and the test results

In the previous section, two models of 3 parameters have been presented. The Poisson ratio is assumed to be constant in the first model, and variable in the second
<figure><img src="image_12.png" /><figcaption>Fig. 10. Evolution of the Poisson ratio.</figcaption></figure>

one. In this section, the verification and validation of the two models will be pre-sented. Numerical adjustments of the parameters of the models were conducted to fit experimental results for different unbound granular materials (Ecuelles, Soreze, Poulmarc'H and Garonne). The results of the adjustments on the Ecuelles and Poulmarc'h materials are presented in Figs. 11, 12 and 13. Also, we can see in these

<figure><img src="image_13.png" /><figcaption>Fig. 11. Evolution of the $K_a$ parameter for different models.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. Evolution of the $G_a$ parameter for different models.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Evolution of the correlation coefficients for different models.</figcaption></figure>

figures the position of Boyce's model parameters. The results of the adjustments indicate that:

* The correlation coefficients, shown in Fig. 14, are greater than 60% (Fig. 14) for the majority of the materials tested except for Ecuelles3. This ensures that the two models are satisfactory for these tests.

* The correlation coefficients obtained with our nonlinear models are not very different from those obtained with Boyce's model. However, in the next section, we will show that beyond the experimental data field, the models are different from each other.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Comparison of the evolution of the Poisson ratio, case of Ecuelles1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Evolution of the Poisson ratio for different models, case of Ecuelles1.</figcaption></figure>

## 5.4. Numerical stability of the proposed models

Fig. 15 shows the agreement and disagreement regions between the proposed nonlinear models and the Boyce's model. For low values of $\eta$ (up to 3), the results indicate good agreement between the two models for all the materials tested. Beyond the experimental data, Poisson ratio tends towards 0.5 for the model with $\gamma = \alpha/2$, and the value is constant for the second model. However, for Boyce's model, Poisson ratio varies significantly and reaches unrealistic values for large values of $q/p$. It is important to notice that in reality, when the ratio $q/p$ becomes large this kind of material undergoes plastic deformations which are not taken into account in this
study. The comments previously made on Poisson ratio also apply to the elastic
modulus. The proposed nonlinear models and Boyce's one compare favorably rela-
tively to the experimental test results. However, beyond this experimental data
range, the results obtained by Boyce's model appear unrealistic.

**6. Conclusions**

This study has shown the ability of a model that deriving from an elastic potential to represents the recoverable strain behavior of unbound granular materials. The method has been based on the determination of the complementary potential values from strain's measurements. in all analyzed cases, the results has shown that the strains can be explained by a complementary potential with a good correlation coefficient. Moreover, the analysis has also shown that it is possible to express the complementary potential as a product of two functions with one depending only on the average pressure *p* and the other on the shear stress *q/p* ratio.

Nonlinear elastic models with three parameters have been proposed to model the
behavior of unbound granular materials. For q/p < 3, these models have shown
similar numerical results to those obtained by the Boyce's model. However, for lar-
ger values of shear stress ratio, both the elastic modulus and the Poisson ratio reach
finite values, which avoid the singularities.

**Acknowledgements**

This study is a part of the second author's thesis accomplished at the Laboratory LGCNSN. During this work we used the experimental results obtained by Mr. J.-L. Paute at the LRPC Laboratory of Saint-Brieuc, France.

**References**

[1] Boyce HR. A non-linear model for the elastic behaviour of granular materials under repeated loading. In: International symposium on soils under cyclic and transient loading, Balkema; 1980. p. 285-94.

[2] Marignier J, Paute JL, Vidal B. Le triaxial à chargement répété lpc pour l'étude des graves non traitées. Bulletin de Liaison des Laboratoires des Ponts et Chaussées 1994;190.

[3] Hornych P, Paute JL, Benaben JP. Comportement mécaniquedes graves non traitées. Bulletin de Liaison des Laboratoires des Ponts et Chaussées 1994;190.

[4] J.L. Paute, P. Hornych and J.P. Benaben. Repeated load triaxial testing of granular materials in the French network of Laboratoires des Ponts et Chaussées. In: Gomes Correia, editor. Proceeding flexibles pavements, Balkema; 1996. p. 53-64.

[5] SETRA-LCPC. Complément de recommendations: réalisation des assises de chaussées en grave non traitées. Technical report, Laboratoire Central des Ponts et Chaussées, Paris; 1980.

[6] SETRA-LCPC. Guide technique: conception et dimensionnement des structures de chaussées. Technical report, Laboratoire Central des Ponts et Chaussées, Paris; 1994.