# Resolution of elastoplastic constitutive relations : application to the fiber reinforced sand

P. Jouve, Rabah Bouzidi, Yvon Riou

► To cite this version:

P. Jouve, Rabah Bouzidi, Yvon Riou. Resolution of elastoplastic constitutive relations: application to the fiber reinforced sand. Computers and Geotechnics, 1996, 17 (3), pp.327-347. 10.1016/0266-352X(95)99216-E. hal-01007255

HAL Id: hal-01007255

https://hal.science/hal-01007255

Submitted on 6 Dec 2016

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# RESOLUTION OF ELASTOPLASTIC CONSTITUTIVE RELATIONS
APPLICATION TO THE FIBER REINFORCED SAND

P. Jouve & R. Bouzidi

Laboratoire de Mécanique et Géomécanique

Université de Nantes

Y. Riou

Laboratoire de Génie Civil

Ecole Centrale de Nantes

GRECO - Géomatériaux : FRANCE

The aim of this paper is to introduce a new formulation for the resolution of elastoplastic constitutive relations. A variant of the initial stress method is proposed to reduce the number of increments and iterations and avoid a drift of results. A validation of this new algorithm is performed in two homogeneous tests : unreinforced sand and fiber reinforced sand in triaxial compression. For the sand, the basic elements used in the study are a double yield surface with isotropic hardening and a non-associated flow rule. Then, in order to check our computer program and the mechanical principle of superposition of the two continua (sand and fibers), a finite element simulation of the construction of a wall under gravity forces is presented.

## INTRODUCTION

In the past decade, major efforts have been made to improve the elastoplastic constitutive relations of soils. In particular, the models describing the behaviour of granular soils are more and more efficient. Nowadays, most structures can be studied by using computers. However, in spite of the capacity of these computers, CPU time is often prohibitive. Clearly, numerical methods are not really efficient in computational time, because the numerical integration of the constitutive laws, connecting stress and strain increments, is not straightforward. This integration requires some tedious incremental and iterative methods, and their convergence in structural simulations is also very critical. Consequently, the industrial use of geotechnical software is still limited.

The aim of this study is to present a more cost-effective numerical method, one that will reduce the computational effort, or make it possible to study some ticklish problems. More precisely, this paper is concerned with a new formulation of the basic elements of the Newton-Raphson's method. A new relation for initial stress is expressed. This relation is introduced in the particular instance of plane strain. The applicability of the model is demonstrated for two homogeneous tests : unreinforced sand and fiber reinforced sand in triaxial testing. A complete structure (a sand-reinforced wall) is then dealt with using the finite element method.
ELASTOPLASTIC CONSTITUTIVE LAWS

Following plasticity theory, the strains are divided into an elastic part and a plastic part, respectively represented by superscripts e and p :

$$
\{\epsilon\} = \{\epsilon^e\} + \{\epsilon^p\} \tag{1}
$$

We assume the following non-linear elastic behaviour :

$$
\{\epsilon^e\} = \{L(\sigma)\} \tag{2}
$$

The elastic strain rates are related to stress rates by a symmetrical matrix :

$$
\{ \dot{\epsilon}^e \} = [D^e]^{-1} \{ \dot{\sigma} \} \quad \text{with:} \quad [D^e] = ( [\nabla_{\sigma} L]^T )^{-1} \tag{3}
$$

[ D^e ] is the elastic tangential matrix ; this matrix is non-singular.

We consider a model characterized by two yield surfaces with isotropic hardening. Let r<sup>j</sup>
be the hardening parameter of mechanism j ; the elastic strains are contained within the current
yield surface in stress space if :

$$
F^j(\sigma, r^j) \le 0.
$$

When plastic yield occurs, the stresses are to be found on the following yield surface :

$$
F^j(\sigma, r^j) = 0 \tag{4}
$$

By differentiating this equation, we can write the consistency condition :

$$
\dot{F}^j = \{\nabla_{\sigma} F^j\}^T \{\dot{\sigma}\} + \frac{\partial F^j}{\partial r^j} \quad \dot{r}^j = 0 \quad (5)
$$

In this case, the description of the strain rate tensor is defined by means of the flow rule :

$$
\{\dot{\epsilon}^{pj}\} = \dot{\lambda}^{j}\{\nabla_{\sigma}G^{j}\} \quad \text{with:} \quad \dot{\lambda}^{j} \geq 0 \tag{6}
$$

where G^j : plastic potential of mechanism j,

$\dot{\lambda}^j$: plastic multiplier of mechanism j.

The classical formulation assumes that the hardening parameter is defined by the plastic
strain rate :

$$
\dot{r}^j = \| \dot{\epsilon}^{pj} \|
$$

The norm $\| \cdot \|$ will be determined later. Thus, we obtain :

$$
\dot{r}^j = N^j \dot{\lambda}^j \quad \text{with} \quad N^j = \| \nabla_{\sigma} G^j \| \tag{7}
$$

If mechanism j is active, relations (4) or (5), (6) and (7) stand for the corresponding plastic
behaviour of the material. If we generalize these relations for all mechanisms, we obtain :

$$
\{F(\sigma, r)\} = \{0\} \tag{8}
$$

$$
\text{or} \quad [\nabla_{\sigma} F]^T \{\dot{\sigma}\} + [\nabla_r F]^T \{\dot{r}\} = \{0\} \quad (9)
$$
The total plastic strain {$\epsilon^P$} = {$\epsilon^P1$} + ... + {$\epsilon^Pn$}, may be expressed as :

$$ \{\epsilon^P\} = [\nabla_{\sigma} G]\{\dot{\lambda}\} \quad \text{with:} \quad \{\dot{\lambda}\} \geq \{0\} \quad (10) $$

In the case of a multimechanism, the hardening parameter tensor can be expressed by the following formula :

$$ \{\dot{r}\} = [N]\{\dot{\lambda}\} \quad (11) $$

If all mechanisms are active, elastoplastic behaviour is entirely defined by relations (1), (2), (8), (10) and (11).

# PROPOSED METHOD OF RESOLUTION

## Basic relations

Let us consider a structure loaded by nodal forces {$P$}, for which plastic yielding has occured. For this state (noted A in figure 1) we consider the "principal" variables at any point : {$\epsilon$}, {$\epsilon^P$}, {$r$} and {$\sigma$}. Because {$\epsilon^e$} is obtained directly by means of {$\sigma$} and equation (2), we call {$\epsilon^e$} a "secondary" variable. The following relations have to be satisfied at each point :

$$ \{\epsilon\} - \{L(\sigma)\} - \{\epsilon^P\} = \{0\} \quad \text{and} \quad \{F(\sigma, r)\} = \{0\}. \quad (12) $$

For this state, the principal variables being obtained by computations, we assume that the previous relations are not entirely satisfied. So, we define the corresponding residuals {$\delta\epsilon$} and {$\delta F$} as follows :

$$ \{\delta\epsilon\} = \{\epsilon\} - \{L(\sigma) - \{\epsilon^P\}\} \neq \{0\} \quad (13) $$

$$ \{\delta F\} = \{F(\sigma, r)\} \neq \{0\} \quad (14) $$

From the initial state to another one (noted B in figure 1), there is a small increment of load {$\Delta P$}. We are going to determine the relation between change of stress and change of strain, and to take account of the residuals {$\delta\epsilon$} and {$\delta F$}.

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Proposed constitutive terms for the incremental flow.</figcaption></figure>
Because of the form of equations (10) and (11), an incremental procedure, with suitably small increments, must be used to solve the equations. Within a loading step, we can write:

$$ \{\Delta\epsilon^P\} = [\nabla_\sigma G]\{\Delta\lambda\} \quad (15) $$

$$ \{\Delta r\} = [N]\{\Delta\lambda\} \quad (16) $$

As we must have $\{\epsilon+\Delta\epsilon\} - \{L(\sigma+\Delta\sigma)\} - \{\epsilon^P+\Delta\epsilon^P\} = \{0\}$, a first order expansion gives:

$$ \{\epsilon\} + \{\Delta\epsilon\} - \{L(\sigma)\} - [D^e]^{-1}\{\Delta\sigma\} - \{\epsilon^P\} - [\nabla_\sigma G]\{\Delta\lambda\} = \{0\} $$

By using $\{\delta\epsilon\}$, expressed in equation (13), we obtain:

$$ \{\delta\epsilon\} + \{\Delta\epsilon\} = [D^e]^{-1}\{\Delta\sigma\} + [\nabla_\sigma G]\{\Delta\lambda\} \quad (17) $$

A similar expression is obtained for $\{F(\sigma+\Delta\sigma, r+\Delta r)\} = \{0\}$, which gives after expansion:

$$ \{F(\sigma, r)\} + [\nabla_\sigma F]^T \{\Delta\sigma\} + [\nabla_r F]^T \{\Delta r\} = \{0\} $$

By using $\{\delta F\}$, expressed in equation (14), it follows:

$$ \{\delta F\} + [\nabla_\sigma F]^T \{\Delta\sigma\} + [H]\{\Delta\lambda\} = \{0\} \quad (18) $$

$$ \text{with: } [H] = [\nabla_r F]^T [N]. $$

The elimination of $\{\Delta\sigma\}$ between the relations (17) and (18) gives the plastic multipliers:

$$ \{\Delta\lambda\} = [W]^{-1} \{[\nabla_{\sigma} F]^T [D^e] \{[\Delta\epsilon] + \{\delta\epsilon\}\} + \{\delta F\}\} \quad (19) $$

$$ \text{with: } [W] = [\nabla_{\sigma} F]^T [D^e] [\nabla_{\sigma} G] - [H]. $$

$$ \text{Let: } [D^{ep}] = [D^e] - [D^e] [\nabla_\sigma G] [W]^{-1} [\nabla_\sigma F]^T [D^e] \quad (20) $$

$$ \text{and } \{\delta\sigma\} = [D^e] [\nabla_\sigma G] [W]^{-1} \{\delta F\} - [D^{ep}] \{\delta\epsilon\}. \quad (21) $$

Thus, the relation between changes of stress and strain is of the form:

$$ \{\Delta\sigma\} + \{\delta\sigma\} = [D^{ep}] \{\Delta\epsilon\} \quad (22) $$

For a given state and for a change of stress $\{\Delta\sigma\}$, $\{\Delta\epsilon\}$ can be calculated by means of equation (22) (see figure 1).

The matrices $[D^{ep}]$ and $\{\delta\sigma\}$ are obtained by the matrices: $[D^e]$, $[\nabla_\sigma F]$, $[\nabla_r F]$ and $[\nabla_\sigma G]$. For each increment, there are two possible ways to calculate these matrices:

a - we update them at the beginning of the stress increment (at point A in figure 1),

b - we update them in the middle of the stress increment.

This last way is more accurate than the first, but it needs two computations for each increment.
The above analysis is more rigorous than the usual ones because :

- we take account of the approximation due to the expansion of {L(σ)},

- we take account of the approximation of a point out of the yield surfaces, without having to project it on to these surfaces.

- the matrix [D<sup>EP</sup>] is non-symmetrical for a non-associated material like a soil, and we use this form for the computation of {ΔU}.

In a structural example, where {P} and {U} stand respectively for nodal forces and nodal displacements, we have the well known relations :

$$ \{\varepsilon\} = [B]\{U\} \quad \text{and} \quad \{P\} = \sum_k \int_{V_k} [B]^T \{\sigma\} \, dV $$

Within each load increment {ΔP} in elastoplastic behaviour, the incremental constitutive relation can be written as :

$$ \{\Delta P\} = \sum_k \int_{V_k} [B]^T \{\Delta \sigma\} \, dV = \sum_k \int_{V_k} [B]^T \{[D^{ep}]\{\Delta \varepsilon\} - \{\delta \sigma\}\} \, dV \\ \text{or} \\ \{\Delta P\} + \{\delta P\} = [K^{ep}]\{\Delta U\} \qquad (23) $$

$$ \text{with } \{\delta P\} = \sum_k \int_{V_k} [B]^T \{\delta \sigma\} \, dV \quad \text{and} \quad [K^{ep}] = \sum_k \int_{V_k} [B]^T [D^{ep}] [B] \, dV \quad (24) $$

[K<sup>ep</sup>] is the tangential stiffness matrix of the structure. We notice that, if [D<sup>ep</sup>] is not symmetrical, [K<sup>ep</sup>] is not either.

## General procedure for solution of structural problems

In a structural example we use an incremental procedure characterized by increments of load {ΔP}. Inside each increment we consider an iterative method based on the equations (22) and (23). The term {δσ} is considered as unbalanced stresses and the term {δP} as an initial load. The principal variables are known at the beginning of the step {ΔP} (or at the end of the previous step) and they are denoted as : {U₀}, {ε<sub>0</sub><sup>P</sup>}, {r₀} and {σ₀}. It follows : {ε₀} = [B] {U₀}, {ε₀<sup>e</sup>} = {L(σ₀)}. We note that the matrices [∇<sub>σ</sub>G₀] and [N₀] are calculated from {σ₀} and {r₀}.

At the iteration k of this loading step, we use the following technique :

a - For each element and from {σ<sub>k-1</sub>} and {r<sub>k-1</sub>}, calculate the matrices : [D<sub>k-1</sub><sup>e</sup>], [∇<sub>σ</sub>F<sub>k-1</sub>], [∇<sub>r</sub>F<sub>k-1</sub>], [W<sub>k-1</sub>] and [D<sub>k-1</sub><sup>EP</sup>], and calculate the unbalanced stresses {δσ<sub>k-1</sub>} from relations (13), (14) and (21);
b - Calculation of the elastoplastic stiffness matrix $[K_{k-1}^{ep}]$ and calculation of the unbalanced nodal forces vector $\{\delta P_{k-1}\}$: relation (24);

c - Resolution of the linear equations system (non-symmetrical for a non-associated law):

for $k = 1$, the system is:

$$ \{\Delta P\} + \{\delta P_0\} = [K_0^{ep}] \{\Delta U_1\}, \text{ where } \{\Delta P\} \text{ is the load increment,} $$

for $k > 1$, the system is:

$$ \{\delta P_{k-1}\} = [K_{k-1}^{ep}] \{\Delta U_k\}, $$

updating of the nodal displacements $\{U_k\} = \{U_{k-1}\} + \{\Delta U_k\}$;

d - For each element, determination of the change of total strain: $\{\Delta\varepsilon_k\} = [B]\{\Delta U_k\}$, calculation of the change of plasticity multipliers $\{\Delta\lambda_k\}$, equation (19), calculation of the changes of plastic strain, hardening parameter and stress:

$$ \{\Delta\varepsilon_k^P\} = [\nabla_\sigma G_0]\{\Delta\lambda_k\}, \quad \{\Delta r_k\} = [N_0]\{\Delta\lambda_k\} \text{ and } \{\Delta\sigma_k\} = [D_{k-1}^{ep}]\{\Delta\varepsilon_k\}; $$

e - Updating of the plastic and total strains, hardening parameters and stresses:

$$ \{\varepsilon_k^P\} = \{\varepsilon_{k-1}^P\} + \{\Delta\varepsilon_k^P\}, \quad \{\varepsilon_k\} = \{\varepsilon_{k-1}\} + \{\Delta\varepsilon_k\}, \quad \{r_k\} = \{r_{k-1}\} + \{\Delta r_k\} \\ \text{and } \{\sigma_k\} = \{\sigma_{k-1}\} + \{\Delta\sigma_k\}; $$

f - Return to “a” for a new iteration when the convergence criterion is not satisfied:

when $\max_i \frac{\|\Delta U_i, k\|}{\|U_i, k\|} \le \mu$, with $\|U_i, k\|$ is the displacement length of the node i, and generally $\mu$ is equal to $10^{-4}$.

This technique is a Newton-Raphson procedure, which satisfies the relations (12) at each step when convergence occurs. The imprecision of the calculation is mainly due to the incremental form for calculating the change of plastic strains during the iterations, the matrix $[\nabla_\sigma G]$ being updated at the beginning of each step. For limiting this imprecision, it is necessary to use small load increments $\{\Delta P\}$.

For a more accurate procedure applied to the increment $\{\Delta P\}$, the best way is to use a second iterative method with the matrix $[\nabla_\sigma G]$ updated at the mean values between the initial stresses and the final ones (obtained at the end of the first iterative process): $(\{\sigma_0\} + \{\sigma_k\})/2$. This is a second order Runge-Kutta procedure which allows load increments greater than the first one with the same accuracy. This procedure is particularly interesting when there is some rotation of the principal directions of the stresses. For the examples presented here, it is not entirely the case. Thus, we only use the first procedure, $[\nabla_\sigma G]$ being updated at the beginning of each load increment.
# PRESENTATION OF THE CONSTITUTIVE MODEL

The incremental procedure for solution and its finite element formulation, described in the 3-dimensional case, are first used in a homogeneous test. We then consider a structure of sand reinforced by continuous fibers. The Vermeer's model is adopted in this paper, to describe the elastoplastic behaviour of sand [2]. This model has been deliberately chosen here because quite a number of successful calculations for granular soils have been reported. In order to specify the version of the model used in this study, a brief description of the formulation is given below.

Five parameters are used to characterize the double hardening behaviour of the sand :

$ε_0^e, ε_0^c, β, φ_p et φ_{cv}$.

## Elastic strain

In contrast to some variants of the model, we consider a non-linear elastic constitutive law, expressed, in the 3-dimensional case, by :

$$ \{\epsilon^e\} = \frac{1}{2} G_s (\sigma) \quad \text{where : } G_s = \frac{3 p_o}{2 \epsilon_o^e} \left(\frac{\sigma_n}{p_o}\right)^{(1-\beta)} \quad \text{and} \quad \sigma_n = \sqrt{\frac{1}{3} \sum_i \sigma_{ii}^2} $$

$p_o$: the reference isotropic compression,

$ε_o^e$: the elastic volumetric strain at isotropic compression $p_o$,

$β$: the basic material constant.

This formulation involves no Poisson's ratio at all.

## Plastic strain

In the Vermeer's model, the plastic strain is subdivided into two independent parts : a deviatoric part and a volumetric part. Each of these strains is defined by means of a yield function F and a plastic potential function G.

### Yield surfaces

Two surfaces act as the boundary of the elastic region in the stress space.

- The first one is cone-shaped around the space diagonal. The yield surface is defined by the equation :

$$ F^1(\sigma, r^1) = I_\sigma \Pi_\sigma + \Pi_\sigma A(\phi_m) = 0 $$

where $I_\sigma, \Pi_\sigma$ and $\Pi_\sigma$ are the three invariants of the stress tensor, and $\phi_m$ is the mobilized
friction related to $\phi_p$. Because of the isotropic hardening, the variable $\phi_m$ is also a function of stress level, and the "generalized" plastic distorsion $r^1$ is a scalar acting as a hardening parameter. The rate of this parameter is defined by:

$$ \dot{r}^1 = \frac{1}{\sqrt{2}} \sqrt{(\dot{\epsilon}_{xx}^{p1} - \dot{\epsilon}_{yy}^{p1})^2 + (\dot{\epsilon}_{yy}^{p1} - \dot{\epsilon}_{zz}^{p1})^2 + (\dot{\epsilon}_{xx}^{p1} - \dot{\epsilon}_{zz}^{p1})^2 + 6[(\dot{\epsilon}_{xy}^{p1})^2 + (\dot{\epsilon}_{yz}^{p1})^2 + (\dot{\epsilon}_{xz}^{p1})^2]} $$
where $\dot{\epsilon}_{ij}^{p1}$ is a plastic strain rate.

In the theory of isotropic hardening, the limit position of the expansion of this yield surface is the failure surface. The Vermeer's model uses the failure criterion of Matsuoka and Nakai expressed by:

$$ F_r^1(\sigma) = I_\sigma \Pi_\sigma + \Pi_\sigma A_r(\phi_p) = 0 $$

$$ \phi_p \text{ is the peak friction angle, and } A_r(\phi_p) = \frac{9 - \sin^2\phi_p}{1 - \sin^2\phi_p} $$

The second yield surface is used to introduce the plastic volumetric strain. This part is not fully explained by the first mechanism. We consider here the spherical yield surface expressed by:

$$ F^2(\sigma, r^2) = \varepsilon_0^c \left(\frac{\sigma_n}{p_0}\right)^\beta - r^2 = 0 $$

with: $p_0$: the reference isotropic compression,

$\varepsilon_0^c$: the plastic volumetric strain at reference isotropic compression $p_0$,

$r^2$: the second hardening parameter.

$$ \dot{r}^2 = \sqrt{3 \Sigma (\dot{\epsilon}_{ij}^{p2})^2} $$

### Plastic potential functions

In the Vermeer's model, the cone-flow rule is non-associated. Here, in order to obtain accurate results for dense and loose sands, and to compare our results with some analytical studies, two plastic potentials are considered:

The first one is formulated on the basis of the yield surface:

$$ G^1(\sigma, r^1) = I_\sigma^* \Pi_\sigma^* + \Pi_\sigma^* A(\psi_m) $$

where $I_\sigma^*$, $\Pi_\sigma^*$, $\Pi_\sigma^*$ are the three invariants of the stresses $\sigma_{ij}^* = \sigma_{ij} - a \delta_{ij}$,

with $\delta_{ij}$ Kronecker's symbol and $a$ a pseudo constant so that $G^1(\sigma, r^1) = 0$.

$$ \psi_m \text{ is the dilatance angle: } \sin \psi_m = \frac{\sin \phi_m - \sin \phi_{cv}}{1 - \sin \phi_m \sin \phi_{cv}}, $$

and $\phi_{cv}$ is the friction angle at constant volume.
The second one is a simple plastic potential, proposed by Vermeer, not involving the pseudo constant a :

$$G^1(\sigma, r^1) = \sqrt{\frac{2}{3} J_{2s}} - \frac{4}{3} p \sin \psi_m$$

with $J_{2s}$ the second invariant of the deviatoric stresses.

- Here we assume that the volumetric yield surface can be identified with the plastic potential: $G^2(\sigma, r^2) = F^2(\sigma, r^2)$.

## SIMULATION OF TRIAXIAL TEST ON UNREINFORCED SAND

To demonstrate the capability of the integration scheme described in the previous section, conventionally drained triaxial compression of unreinforced sand is simulated. This first example is studied with smooth contact between sample and platens. Consequently, homogeneous stress and strain states are generated in the whole sample. The aim is to test the procedure of resolution of elastoplastic equations, not the validity of the behaviour model. Three factors are examined and presented below:

- accuracy of solution when the law equations may be integrated by an analytical method,

- convergence expressed by a number of iterations,

- stability of the scheme with regard to the loading process.

In this study, three compression tests, for which numerical results obtained by Vermeer [2] and Mestat [3] are available, are performed:

- isotropic compression,

- constant stress-ratio path,

- conventional triaxial compression.

### Isotropic compression

According to the resolution procedure, the load is applied in several steps from initial state, defined by:

$$\sigma^0_1 = \sigma^0_2 = \sigma^0_3 = 200 \text{ kPa}$$

To solve this problem, a quantitative parameter study with four loading processes is carried out to verify the numerical performance of this algorithm: 2, 5, 10 and 20 increments. For a given simulation, the variation of stress is constant. In order to describe the sand behaviour, the selected values are:

$$p_0 = 200 \text{ kPa}, \varepsilon^0_0 = 0.004, \varepsilon^0_0 = 0.0035, \beta = 0.25, \phi_p = 43^\circ, \phi_{cv} = 32^\circ$$
These basic material constants have been proposed at the Grenoble Workshop which fit with dense Karlsruhe sand [2].

In this example, because of homogeneous stresses generated in the whole sample, the convergence criterion concerns the unbalanced stresses, and not the unbalanced nodal forces. As a general rule for all the simulations in this application, the convergence is achieved when the norm of this vector is less than 0.01 Pa :

$$|| \delta\sigma || = \sqrt{(\delta\sigma_1)^2 + (\delta\sigma_2)^2 + (\delta\sigma_3)^2} < 0.01 \text{ Pa}$$

with {$\delta\sigma_i$} : the unbalanced principal stress vector.

These results (stresses and strains) obtained during the loading process, can be verified. Indeed, for such a stress path, the strain rate equation may be integrated analytically :

$$\epsilon = \frac{\epsilon_0^e + \epsilon_0^c}{3} \left[ \left( \frac{\sigma}{\rho_0} \right)^\beta - 1 \right]$$

Table 1 and figure 2 show the numerical results for various loading process.

<table><caption>TABLE 1 : Isotropic compression test axial, 212 strain at 210 kPa for various loading processes.</caption><thead><tr><th>simulation</th><th>number of increments</th><th>number of iterations</th><th>axial strain at 210 kPa</th></tr></thead><tbody><tr><td>1</td><td>2</td><td>3</td><td>0.306803 10<sup>-4</sup></td></tr><tr><td>2</td><td>5</td><td>3</td><td>0.306803 10<sup>-4</sup></td></tr><tr><td>3</td><td>10</td><td>3</td><td>0.306806 10<sup>-4</sup></td></tr><tr><td>4</td><td>20</td><td>3</td><td>0.306806 10<sup>-4</sup></td></tr><tr><td>analytical solution</td><td></td><td></td><td>0.306806 10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 : Isotropic compression test, loading from 200 kPa to 210 kPa.</figcaption></figure>

## Constant stress-ratio path

In this test, we consider a constant stress-ratio path starting from the initial stress state :

$$\sigma_1^e = 200 \text{ kPa} \quad \sigma_2^e = \sigma_3^e = 100 \text{ kPa}.$$

Then proportional stressing is applied : $\sigma_2 = \sigma_3 = 0.5 \sigma_1$.

As we have already mentioned, this load is applied with different numbers of increments in order to test the reliability of the solution. The previous Vermeer's model parameters are also
used in this example. To evaluate the hardening parameters, we consider a stress state on the yield surfaces. The compression test then starts with plastic strains. For this proportional loading and the simple plastic potential function of the second mechanism, the flow rule equation can be integrated. This leads to the following relation [2], [3]:

$$ \epsilon_1 = \epsilon_1^e + \epsilon_1^{P1} + \epsilon_1^{P2} = S(\sigma_1^\beta - \sigma_1^{\beta0}) \quad \text{with} \quad S = 1.00132 \cdot 10^{-6} \quad \text{and} \quad \beta = 0.25, $$

for current values of the Vermeer's model parameters. The results of the numerical model are presented in table 2 and figure 3.

<table><caption>TABLE 2: Constant stress-ratio path, axial strain at 210 kPa for various loading processes.</caption><thead><tr><th>simulation</th><th>number of increments</th><th>number of iterations</th><th>Total axial strain at 210 kPa</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>3</td><td>4.6212 10<sup>-5</sup></td></tr><tr><td>2</td><td>2</td><td>3</td><td>4.6213 10<sup>-5</sup></td></tr><tr><td>3</td><td>5</td><td>3</td><td>4.6213 10<sup>-5</sup></td></tr><tr><td>4</td><td>10</td><td>3</td><td>4.6212 10<sup>-5</sup></td></tr><tr><td>analytical solution</td><td></td><td></td><td>4.6212 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIGURE 3: Constant stress-ratio path, loading from 200 kPa to 210 kPa.</figcaption></figure>

### Triaxial compression test : $\sigma_3$-constant test

To further evaluate the capability of the resolution procedure, the axial strain of the specimen subjected to constant lateral stress is predicted:

$$ \sigma_2^0 = \sigma_3^0 = 100 \text{ kPa}, \sigma_1^0 \text{ increasing from } 100 \text{ kPa to failure.} $$

The Vermeer's model parameters are the same as the previous test but are applicable for a reference stress $p_0$ of 100 kPa. Only in terms of axial strains, this stress path involves a numerical process of resolution. So, we have made a comparison with the results of a classical resolution using the simple potential function of the first mechanism [3]. The axial strain at 420 kPa obtained by the two methods is reported in table 3.

It was then decided to push the solution almost to the failure point. For four loading processes, the results of the numerical model are represented in the $\sigma_a-\epsilon_a$ and $\epsilon_v-\epsilon_a$ diagrams (see figure 4).
<table><caption>TABLE 3 : Triaxial compression test, axial strain at 420 kPa for various loading processes.</caption><thead><tr><th rowspan="2">simulation</th><th rowspan="2">number of increments</th><th rowspan="2">average number of iterations</th><th colspan="4">Axial strain at 420 kPa</th></tr><tr><th>elastic strain (10<sup>-3</sup>)</th><th>cone-plastic strain (10<sup>-2</sup>)</th><th>cap-plastic strain (10<sup>-4</sup>)</th><th>total (10<sup>-2</sup>)</th></tr></thead><tbody><tr><td>1</td><td>5</td><td>6</td><td>1.4348</td><td>1.2415</td><td>4.2851</td><td>1.4278</td></tr><tr><td>2</td><td>20</td><td>5</td><td>1.4348</td><td>1.2328</td><td>4.4247</td><td>1.4205</td></tr><tr><td>3</td><td>80</td><td>4.8</td><td>1.4348</td><td>1.2274</td><td>4.4632</td><td>1.4155</td></tr><tr><td>4</td><td>320</td><td>3.8</td><td>1.4348</td><td>1.2258</td><td>4.4730</td><td>1.4140</td></tr><tr><td>program "Lois"</td><td>32000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1.4136</td></tr></tbody></table>

For all these cases, we note that the obtained solution is very regular until a maximum distortion of 15 %. This distortion has been compared with the analytical solution :

$$ \gamma = B \frac{C \frac{q}{p}}{C - \frac{q}{p}} \quad \text{with} \quad B = \frac{p_0}{2 G_0} \left(\frac{p}{p_0}\right)^\beta, \quad C = \frac{6 \sin \phi_p}{3 - \sin \phi_p} $$

The percentage of error on total distortion at the last calculated load point is less than 1%.

This test is specific because the stress field is homogeneous, the loading is regular and the plastic flow direction is almost constant. So, it cannot stand as a representative application of the algorithms performance. It stands as a first check of the computer program. However, this agreement in different simulations indicates the applicability of the algorithm. Contrary to the initial stress method with a constant elastic matrix and a simple formulation of the unbalanced stress, we note here a limited number of iterations. Furthermore, although the plasticity relations are incremental, it is interesting to note that large load increments yield close results. This would increase efficiency and economy of analysis in structural problems.

After the new algorithm of numerical integration is tested and to further evaluate the capability of the program presented here, the behaviour of a reinforced sand is predicted. It must be kept in mind that our computer program was originally written to solve the 2-dimensional problems of micro-reinforced structures.
<figure><img src="image_5.png" /><figcaption>FIGURE 4: Triaxial Compression Test, Effect of loading process on Distortion.</figcaption></figure>

SIMULATION OF TRIAXIAL TEST ON REINFORCED SAND

The procedure of iteration and the constitutive model are now used to predict the stress-strain behaviour of sand reinforced by a network of flexible and continuous fibers. Like previously, the test consists in drained conventional triaxial compression with constant lateral stress. A similar analysis in uniform stress field has been performed experimentally and numerically by C. di Prisco and R. Nova in 1991 [4]. Their approach will not be discussed here as it involves a different procedure. However, in this section, their numerical and experimental results are plotted on figure 7 and compared with the simulations yielded by our program.

The formulation of the constitutive law of sand reinforced by continuous threads is given in reference [5]. Namely, the model proposed here may be considered as a continuum obtained by superimposing two continua. These continua are distinct but closely linked: solid medium and 3D network. It might be interesting to note that this assumption can reflect upon the simulation at failure state. Indeed, experimental testing indicates that the collapse of material can occur by disconnection of these continua or by failure of fibers in a shear band. These mechanisms are not considered here.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 : Sample of Reinforced Sand.</figcaption></figure>

The two continua obey the laws of solid mechanics, and undergo the same strain :

- The threads are very thin and they can only bear tensile forces. Their behaviour studied from tensile tests, is considered as elastic non linear.

<figure><img src="image_7.png" /><figcaption>FIGURE 6 : Constitutive Law of Threads.</figcaption></figure>

The volumetric ratio of threads is expressed by $\rho$. The distribution of the directions of the threads can be fixed by the users. This particularity accounts for the anisotropy of the material and can favour a bedding plane, depending on the way in which it is put in place : see figure 5.b. In order to compare our results with those of C. di Prisco and R. Nova, we assume that the reinforcement acts isotropically. The following thread properties are assumed (see figure 6) :

$E_0$ : Young's modulus $E_0$ = 10 300 MPa,

$\sigma_m$ : "threshold" stress $\sigma_m$ = 610 MPa,

$\rho$ : volumetric ratio $\rho$ = 0.1 %.
- To describe the sand behaviour, we use the elastoplastic Vermeer's model, previously presented. The five model constants are determined from the triaxial compression tests data of C. di Prisco and R. Nova, and by a least squares process. The best fit with the model test data, up to a distortion of 15%, are obtained with:

$$p_0 = 200 \text{ kPa}, \varepsilon_0^e = 0.004, \varepsilon_0^c = 0.0035, \beta = 0.45, \phi_p = 41^\circ, \phi_{cv} = 31^\circ$$

The comparison between simulation and experimental axial strain is shown in figure 7.

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Comparison between measured and predicted data, for a constant confining stress.</figcaption></figure>

As can be seen from figure 7, the numerical solution shows a reasonable agreement with the experimental data. The axial strain is correctly predicted up to 8-10%. Beyond this value, our previous assumptions concerning the constitutive law of the composite assembly, in particular the connection of sand and threads, are no longer valid. It is important to note that the discrepancy between the two predictions is due to an elastic perfectly plastic behaviour of materials assumed by C. di Prisco and R. Nova.

It may be noted that disagreements appear in the prediction of volumetric strain. C. di Prisco and R. Nova pointed out the same divergent results with model test data, implicating the simple model and (or) the reliability of the measurements of volumetric strain. So, we are faced with the same problem. A new triaxial test will be needed to explain the reasons for this difference.

Before dealing with structures, we give some explanations about the mechanical principle of reinforcement. During the isotropic compression, the threads are compressed and can be ignored. In deviatoric loading, the deformation in a given direction induces a reaction of
the threads in tension. To a lesser extent, this reaction acting in an angular sector increases with the dilatancy of the sand. Thus a fictitious confinement operates on the sand, allowing the deviatoric stress to increase : see figure 8 "Stress in sand of Texsol". The reinforcement can be explained by the Q-P diagram : according to previous remarks on modelling and model test data, we have considered here that failure occurs at an axial strain of 6,8 or 10 % depending on lateral pressure. Figure 9 indicates a "fictitious" cohesion (65 kPa), experimentally provided and used in reinforced sand structures design.

<figure><img src="image_9.png" /><figcaption>FIGURE 8 : Mechanical principle of reinforcement, P-Q diagram of drained triaxial test on unreinforced sand and reinforced sand : $\sigma_r$ = 100 kPa</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9 : Mechanical principle of reinforcement, Failure surface of unreinforced sand and reinforced sand.</figcaption></figure>
INITIAL STRESS DISTRIBUTION IN A REINFORCED SAND WALL

This problem involves the construction phase of a reinforced sand wall, assumed to be supported by an embedded foundation. The geometry of the selected wall is presented in figure 10. We study the stress distribution of the structure subjected to the force of its own weight, prior to lateral earth pressure, under plane strain conditions. In order to estimate the weight stresses, we propose a numerical procedure. According to the requirement of the resolution process, the calculation is performed in several increments, each increment corresponding to a percentage of the weight.

<figure><img src="image_11.png" /><figcaption>FIGURE 10: Reinforced sand wall - geometry, initial stress and loading.</figcaption></figure>

During the construction of the wall, compaction and weight have led to a prestress of the fibers. Moreover, this prestress is necessary to start numerical resolution of this problem. According to a previous simulation, we consider here an isotropic initial stress (9kPa) in the sand corresponding to a prestress in the fiber induced by compaction. This initial "guess" is necessary for ensuring convergence of the iterative process of resolution, and small enough not to affect the final distribution of stress. The loading consists in solving the gravity stress field.

In order to avoid a stress state exceeding the yield condition at the early increments, we consider an arithmetic progression of gravity forces:

$$ \{F_v\} = \sum_{i=1}^{N} a_i \{F_v\} \quad \text{with} \quad a_i = \frac{2i-1}{(N)^2} $$
N: number of increments, $\{F_v\}$: total gravity forces.

As mentioned earlier, the convergence criterion of the resolution process is expressed by the unbalanced stress:

$$ \|\delta\sigma\| = \sqrt{(\delta\sigma_1)^2 + (\delta\sigma_2)^2 + (\delta\sigma_3)^2} < 0.01 \text{ Pa} $$
Vermeer's model and the hyperbolic model are used to model the sand and the threads. We have adjusted the constant parameters to the components of the geotechnical material originally conceived by LCPC: Texsol [6]. Using a least squares method in regression computation, the sand properties are taken as:

$$p_0 = 200 \text{ kPa}, \varepsilon_0^e = 0.004, \varepsilon_0^c = 0.0035, \beta = 0.25, \phi_p = 42^\circ, \phi_{cv} = 32^\circ$$

The following thread properties are assumed:

$$E_0 = 12400 \text{ MPa}, \sigma_m = 550 \text{ MPa}, \rho = 0.2 \%$$

According to the observations on Texsol structures, it is assumed that the threads are deposited in a 3D angular sector of 30°. The bedding plane is horizontal. The mesh of the wall consists of 64 8-noded isoparametric elements with 3*3 integration (see figure 11).

<figure><img src="image_12.png" /><figcaption>FIGURE 11: Mesh of the reinforced sand wall.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12: Reinforced sand wall, variations of the settlement at the top.</figcaption></figure>
The results of the numerical evaluation are given in figure 12. This figure shows the settlement of the top of the wall for various loading processes. For all these cases, the obtained solution is very regular. Furthermore, a comparison of the loading considering 20, 30 and 50 increments, shows that the distributions are almost identical. For each increment in the last loading (50 increments), only five iterations are needed for convergence of the iterative procedure. It is important to note that a 10 increment loading leads to sizeable differences in the stress distribution. The behaviour of the wall subjected to lateral pressure could be affected by this different initial stress state: see figure 13.

<figure><img src="image_14.png" /><figcaption>FIGURE 13: Reinforced sand wall: contours of vertical stress for various loading (kPa).</figcaption></figure>

Figure 14 shows the horizontal stress distribution in the anisotropic medium equivalent to the network of fibers. This stress state resulting from gravity complements the initial stress due to the isotropic compression. Maximum value of horizontal component of tension in medium equivalent to fibers is then 17 kPa (9 kPa for initial stress and 8 kPa for weight). According to the density of fibers and spatial distribution, their maximum stress value is 8.5 MPa. This tension must be compared with the ultimate limit stress of the fibers: $\sigma_m$ = 550 MPa. This example demonstrates that the network of fibers is weakly activated.
<figure><img src="image_15.png" /><figcaption>FIGURE 14 : Reinforced sand wall : distribution of stress in medium of fibers (kPa).</figcaption></figure>

The aim of this simulation was to present a numerical initial stress state of the components of the reinforced sand, for a ground stress field. For the moment, we have no experimental data at our disposal, justifying this computation. However computational results are in agreement and indicate that the simulation using the previous resolution process is feasible. Furthermore, convergence and accuracy are assumed for a reasonable number of increments and iterations.

## CONCLUSION

A former approach to geotechnical problems with a conventional resolution process, was not really encouraging. It was found that the resolution required very small steps to ensure convergence stability. In this paper, a new approach of analysing reinforced sand walls has been presented. The formulation, the implementation and the validation of the new resolution process have made it possible to acquire thorough knowledge of the plastic mechanisms and a
full analysis of numerical problems. The good agreement between analytical and computed values for different loadings, demonstrates the validity of the new method which seems to be more economical than the conventional initial stress method.

Our program needs some improvements for the specific simulation of fiber reinforced sand structures. However, in this version, it can be used for parametric study. This kind of study will be needed in order to fully validate the assumptions of the mix of sand and fibers, and more particularly the interaction of these two components. It seems especially necessary to adjust the number and direction of the fibers acting in reinforcement, for the study of structures. This study requires a comparison between the simulations and measures carried out on the structures. This is our current aim taking as experimental references two retaining walls subjected to lateral pressure and fitted with measuring instruments by the Centre d'Expérimentation Routière de Rouen (France).

## ACKNOWLEDGEMENTS

The French program "GRECO Géomatériaux" is gratefully acknowledged for its financial support.

## REFERENCES

1. Zienkiewicz O.C.,
   *The finite element method*, 3rd ed., Mac-Graw Hill (1977).

2. Vermeer P.A.,
   A five-constant model unifying well-established concepts,
   *Results of the Internat. Workshop on Constitutive Relations for Soils*, Grenoble (1982).

3. Mestat P.,
   Validation numérique du procédé CESAR-LCPC pour les lois de comportement élastique non-linéaire plastique,
   *Rapport Intersociété LCPC*, 1.16.21.0 (1990).

4. Di Prisco C. and Nova R.,
   A constitutive model for soil reinforced by continuous threads,
   *Rapport scientifique de l'EURO-GRECO Géomatériaux*, Paris (1991).

5. Villard P., Jouve P. and Riou Y.,
   Modélisation du comportement mécanique du Texsol,
   *Bulletin de Liaison des Laboratoires des Ponts et Chaussées*, 168 (1990) 15-27.

6. Khay M. and Gigan J.P.,
   Texsol ouvrage de soutènement,
   *Guide technique du L.C.P.C.*, Paris (1990).