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

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 fmite element simulation of the construction of a wall under gravity forces is presented.


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: We assume the following non-linear elastic behaviour: {£") = {L(a)} (2) The elastic strain rates are related to stress rates by a symmetrical matrix :  (3) [ oe] is the elastic tangential matrix ; this matrix is non-singular.
We consider a model characterized by two yield surfaces with isotropic hardening. Let ri be the hardening parameter of mechanism j ; the elastic strains are contained within the current yield surface in stress space if : Fi(a, d) ::;; 0.
When plastic yield occurs, the stresses are to be found on the following yield surface : By differentiating this equation, we can write the consistency condition : ·. ( ·)T(. ) ()pi · FJ = V 0 P a + -. rl = 0 drJ (5) In this case, the description of the strain rate tensor is defined by means of the flow rule : (f:Pi} = j_i (V 0 Gi) with: i_i ~ 0 (6) where Gi : plastic potential of mechanism j, ).i : plastic multiplier of mechanism j.
The classical formulation assumes that the hardening parameter is defmed by the plastic strain rate : The norm II II will be determined later. Thus, we obtain : ti = Ni j_i with Ni = IIV aGill (7) If mechanismj 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 :   [ Thus, the relation between changes of stress and strain is of the form : For a given state and for a change of stress ( Acr} , ( Ae} can be calculated by means of equation (22) (see figure 1). 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. 1be above analysis is more rigorous than the usual ones because : we take account of the approximation due to the expansion of { L( a)}, -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 [ l)t'P) is non-symmetrical for a non-associated material like a soil, and we use this form for the computation of {AU} .
In a structural example, where {P} and (U} stand respectively for nodal forces and nodal displacements, we have the well known relations : c -Resolution of the linear equations system (non-symmetrical for a non-associated law) : . fork= 1, the system is: . fork> 1, the system is: 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 [ V cP ] being updated at the beginning of each step. For limiting this imprecision, it is necessary to use small load increments { AP } .
For a more accurate procedure applied to the increment {AP}, the best way is to use a second iterative method with the matrix [ V cP ] updated at the mean values between the initial stresses and the fmal ones (obtained at the end of the frrst iterative process) : 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.

PRESENTATION OF THE CONSTITJJTIYE MODEL
1be 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 nwnber 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. This formulation involves no Poisson's ratio at all.

Plastic strajn
In the Venneer'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.
-1be first one is cone-shaped around the space diagonal. The yield surface is defined by the equation : friction related to IPp· Because of the isotropic hardening, the variable ~Pm is also a function of stress level, and the "generalized" plastic distorsion r 1 is a scalar acting as a hardening parameter. 1be rate of this parameter is defined by : 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 : !pp is the peak friction angle, and A r (!pp) = 9 -sin 2 1Pp 1 -sin 2~P p -The second yield surface is used to introduce the plastic volwnetric strain. This part is not fully explained by the first mechanism. We consider here the spherical yield surface expressed by : Po with : p 0 : the reference isotropic compression, ec 0 : the plastic volumetric strain at reference isotropic compression p 0 , r2 : the second hardening parameter. r 2 = ,.j 3 l: (Erff 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 : where I;,., ~. IT~ are the three invariants of the stresses cr; j = Oi j -a Bi j , with ~j Kronecker's symbol and a a pseudo constant so that G 1 (cr, r 1 ) = 0. -Here we assume that the volumetric yield surface can be identified with the plastic potential: G2(a, r2) = F2(a, r2) .

SIMULATION OF TRIAXIAL TEST ON JJNREINFORCED SAND
To demonstrate the capability of the integration scheme described in the previous section, conventionally drained triaxial compression of unreinforced sand is simulated. This frrst 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.

lsotropjc compressjop
According to the resolution procedure, the load is applied in several steps from initial state, defined by : 00 1 = 00 2 = 00 3 = 200 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. 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 : llooll = 'V'(oo.) 2 +(oo2) 2 +(oo3) 2 < o.o1 Pa with { &Ji} : 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 :  figure 2 show the numerical results for various loading process.  In this test, we consider a constant stress-ratio path starting from the initial stress state :

Irjaxjal compression test : ~-copstapt test
To further evaluate the capability of the resolution procedure, the axial strain of the specimen subjected to constant lateral stress is predicted : 00 2 = 00 3 = 100 kPa , a 1 increasing from 100 kPa to failure.
The Venneer' s model parameters are the same as the previous test but are applicable for a reference stress Po 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.   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: 2Go Po 3-sin~ 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.

SIMULATION OF TRIAXIAL TEST ON REINFORCED SAND
1lte procedure of iteration and the constitutive model are now used to predict the stressstrain 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 30 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. The volumetric ratio of threads is expressed by p. 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 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.

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 sttess 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 sttesses, 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 geometry isotropic initial stress : 9kPa loading: gravity force, 18 kN/m3 FIGURE 10: Reinforced sand wall-geometry, initial stress and loading.
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 fmal 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 : Vermeer's nxxlel and the hyperbolic nxxlel are used to model the sand and the threads. We have ajusted 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 : Po = 200 kPa , Eeo = 0.004, Eco = 0.0035 , ~ = 0.25 , ' P = 42°, «Pcv = 32~ The following thread properties are assumed : Eo= 12 400 MPa, Om= 550 MPa, p = 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). 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. 10 increments 20 increments 30 increments 50 increments FIGURE 13: Reinforced sand wall: contours of vertical stress for various loading (kPa). 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 : am = 550 MPa. This example demonstrates that the network of fibers is weakly activated. 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