Zhen-Yu Yin

Pierre-Yves Hicher

Identifying parameters controlling soil delayed behaviour from laboratory and in situ pressuremeter testing

Keywords: viscoplasticity, consolidation, Modified Cam Clay, parameter identification, soft clay, pressuremeter test

The aim of this paper is to present a methodology for identifying the soil parameters controlling the delayed behaviour from laboratory and in situ pressuremeter tests by using an elasto-viscoplastic model (EVP-MCC) based on Perzyna's overstress theory and on the elasto-plastic Modified Cam Clay model. The influence of both the model parameters and the soil permeability was studied under the loading condition of pressuremeter tests by coupling the proposed model equations with Biot's consolidation theory. On the basis of the parametric study, a methodology for identifying model parameters and soil permeability by inverse analysis from three levels of constant strain rate pressuremeter tests was then proposed and applied on tests performed on natural Saint-Herblain clay. The methodology was validated by comparing the optimized values of soil parameters and the values of the same parameters obtained from laboratory test results, and also by using the identified parameters to simulate other tests on the same samples. The analysis of the drainage condition and the strain rate effect during a pressuremeter test demonstrated the coupled influence of consolidation and viscous effects on the test results. The numerical results also showed that the inverse analysis procedure could successfully determine the parameters controlling the time-dependent soil behaviour.

INTRODUCTION

The identification of soil parameters from in situ tests, such as pressuremeter tests consisting in studying the expansion of a cylindrical cavity within the soil, has been commonly used since the invention of the pressuremeter apparatus by Menard [1]. Zentar et al. [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF] divided the methods used for this purpose into three categories: (a) those based on empirical or semi-empirical relationships between soil properties and pressuremeter parameters; (b) those based on the theoretical development of cylindrical cavity expansion solutions; and (c) those based on numerical simulations of the test. It is noted that the proper calibration of soil parameters depends on the number of parameters to be optimized for a selected model and the number of observations, as shown by Zentar et al. [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF] during the optimization procedure and as indicated by Calvello and Finno [START_REF] Calvello | Selecting parameters to optimize in model calibration by inverse analysis[END_REF].

In the past few decades, several studies for identifying soil parameters from pressuremeter tests have been found in the literature, i.e. Anderson et al. [START_REF] Anderson | Rate effects in pressuremeter tests in clays[END_REF], Bahar et al. [START_REF] Bahar | Forecast of creep settlements of heavy structures using pressuremeter tests[END_REF], Pye [START_REF] Pye | Influence of constitutive models on self-boring pressuremeter interpretation in clay[END_REF], Hicher and Michali [START_REF] Hicher | Identifying soil parameters by means of laboratory and in situ testing[END_REF], Zentar et al. [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF], Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF], etc. However, few of these are concerned by the delayed behaviour of fine soils. However, the strain-rate dependency of soft soils has a significant influence on the selection of design parameters. A constitutive model taking into account the time-dependent behaviour should be recommended, although additional parameters of viscosity for this kind of model would result in difficulties during the selection of the soil parameters by optimization procedure. In order to study this kind of problem, we based our analysis on experimental results obtained with the equipment called pressio-triax [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF][START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF], which could simulate a self-boring pressuremeter (SBP) with pore water pressure measurement at the cavity wall in tightly controlled laboratory conditions.

In the first part of this paper, we briefly present the development of an elasto-viscoplastic model with its implementation in a finite element program CESAR LCPC. The motivation behind this initial work lied in the fact that we wished to have a viscoplastic version of the Modified Cam Clay model for which parameter identification procedures had already been developed [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF][START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF]. Then a back analysis method for identifying soil mechanical and hydraulic parameters from pressuremeter tests with three different strain rates is presented and applied on tests performed on natural Saint-Herblain clay. The methodology is validated by comparing the optimized values of soil parameters with the values of the same parameters obtained from laboratory test results, and also by using the identified parameters to simulate other tests on the same samples. The drainage condition during a pressuremeter test as a function of the strain rate and the permeability as well as the strain rate effect on the pressuremeter curve is analysed. Finally, an application to in situ pressuremeter tests is carried out.

ELASTO-VISCOPLASTIC CONSTITUTIVE MODEL

An elasto-viscoplastic model, named EVP-MCC, based on the framework of Perzyna's overstress theory and Modified Cam Clay elasto-plastic model was developed by Yin et al. [START_REF] Yin | An elasto-viscoplastic model for soft clay[END_REF]. The proposed model has two additional viscosity parameters besides the parameters of the Cam clay model and is an attempt to provide a simple but realistic enough approach in modelling time-dependent behaviour of clayey soils subjected to monotonic proportional loading paths. Some limitations of this model are discussed in the following sections.

According to the assumption of Perzyna's overstress theory [START_REF] Perzyna | The constitutive equations for work-hardening and rate sensitive plastic materials[END_REF][START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF], the viscous effects are negligible in the elastic region. In other words, the elastic strains are time independent, whereas the inelastic strains are time dependent. We define a static yield criterion f s , which represents a reference yield surface for the material, corresponding to the elastic limit with null viscoplastic flow rate. Its initial shape depends on the consolidation pressure p s c . The expansion of the static yield surface, which represents the hardening of the material, is expressed by the variation in the consolidation pressure due to the inelastic volumetric strain vp v :

d p s c = p s c • 1+e 0 - •d vp v = p s c * •d vp v (1)
where p s c is the static consolidation pressure, vp v the inelastic volumetric strain, and * the compressibility index.

A dynamic loading yield criterion f d is defined to represent the current state of stress and is expressed as follows:

f d = q 2 M 2 + p •( p -p d c ) = 0 (2)
where p d c is the dynamic consolidation pressure, M the slope of the critical state line, p representing the effective mean stress p = ( 11 + 22 + 33 )/3, and q the deviatoric stress, which is given by

q = 3J 2 = [( 11 -22 ) 2 +( 11 -33 ) 2 +( 22 -33 ) 2 +6( 2 12 + 2 13 + 2 23 )]/2
We then define a distance F between the dynamic loading and the static yield surfaces at a given time t as F = p d c / p s c . Based on the distance F, the scaling function (F), which controls the amplitude of the viscoplastic strain rate, is taken as an exponential form as follows:

• (F) = •(exp[N •(F -1)]-1) (3)
where and N are the viscosity parameters of the model. As expressed in Equation ( 3), the model cannot describe the volumetric relaxation process, unless the initial state is not considered at equilibrium. The flow rule for the viscoplastic strain rate, in a simple case of infinitesimal strain field, follows the form proposed by Perzyna [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF]:

˙ vp i j = (F) * f d * i j (4)
where the function of MacCauley is

F = 0 for F 0 F for F>0
The principles of the elasto-viscoplastic model EVP-MCC are illustrated by the effective stress path of an undrained triaxial test presented in Figure 1. The stress state 'A' represents an initially normally K 0 consolidation state. Along the loading stress path ' A-B-C', viscoplastic volumetric strains occur during loading and cause the static yield surface to expand in the stress space. As point C approaches C , corresponding to the critical state, the soil is subjected to a constant amount of overstress, which provokes an increase in the deviatoric strain at constant strain rate, without any volumetric strain.

As for the creep stress path 'B-D', the static yield surface expands with the time-dependent plastic volumetric strain as a function of the amount of overstress. If the static yield surface can reach the actual stress point, the equilibrium is obtained and the strain will be stabilized with time.

If not, the effective stress will continue to evolve until it reaches the critical state at point D where it will stop because no plastic volumetric strain will develop, but deviatoric strain will continue to increase.

Taking into account the elastic stress-strain relationships, the constitutive equations of the viscoplastic model for normally consolidated clays are derived as follows:

˙ i j = ṡ i j 2G + ṗ 3K i j + (F)
3s i j M 2 +(2 p -p d c ) i j 3 (5)
where i j is the effective stress tensor and ˙ vp i j the inelastic strain rate tensor. The basic finite element scheme for the proposed model is similar to that presented by Oka et al. [START_REF] Oka | Two-dimensional consolidation analysis using an elasto viscoplastic constitutive equation[END_REF] and Hinchberger and Rowe [START_REF] Hinchberger | Modelling the rate-sensitive characteristics of the Gloucester foundation soil[END_REF]. For a coupled consolidation analysis based on Biot's theory, the relationship of the load increment is given by applying the principle of virtual work to the equilibrium equation as shown by Oka et al. [START_REF] Oka | Two-dimensional consolidation analysis using an elasto viscoplastic constitutive equation[END_REF]:

V [ ] T •[ ] dV = V [ ] T •[ ] dV + V [ ] T •[ u w ] dV = V [ d] T •[ F b ] dV + S [ d] T •[ T ] dS (6)
The equation of continuity for water flow can be expressed as follows:

vol

[b] 3 i=1 k i w * 2 u w *x 2 i + *v *t d(vol) = 0 (7)
where k i is the hydraulic conductivity in three coordinate directions, w is the unit weight of the pore water, and v is the seepage velocity. The coupled finite element equations were well documented by several researchers (e.g. [START_REF] Oka | Two-dimensional consolidation analysis using an elasto viscoplastic constitutive equation[END_REF][START_REF] Hinchberger | Modelling the rate-sensitive characteristics of the Gloucester foundation soil[END_REF][START_REF] Britto | Critical State Soil Mechanics via Finite Elements[END_REF]).

IDENTIFICATION OF SOIL PARAMETERS FROM PRESSUREMETER TESTS

Methodology of parameter identification by inverse analysis

The mechanical problem is traditionally resolved by calculating the response (R) of a mechanical system (S) including the constitutive model (M) and its parameters (P) subjected to the actions (A).

It corresponds to what is called the direct problems (Figure 2(a)) and is expressed in mathematical terms as follows:

R = F(S) (8)
However, in inverse analysis (Figure 2(b)), the parameters (P) of the given model (M) should be calibrated by using an optimization procedure that consists of iteratively changing the input values of the parameters until the calculated results match the observed data, i.e. laboratory or in situ tests results. The methodology of the inverse analysis is shown in Figure 3. During the optimization procedure, the difference between observation data and model prediction is expressed as follows:

L n (P) = 1 t 1 -t 0 R * (t)-R(S, t) dt (9)
where the notation . . . represents a norm in the space variable, t 1 -t 0 is the time of observation, and R * (t)-R(S, t) is the difference between experimental and numerical data. Owing to the measurements at discrete moments, the integral of L n (P) can be replaced by a sum and the length of observation by the number of measurements. The difference between measured and predicted data is then expressed as follows: L n (P) < ( is a given tolerance) can be used as a discriminate function to judge the convergence of the iterative process.

L n (P) = 1 M n M n i (R * i -R i ) 2 (10) A R S (M+P) (a) A R S (M+P) (b)

Experimental results

The pressuremeter tests were conducted by using a special pressuremeter apparatus called pressiotriax, which allows reproducing the pressuremeter test conditions in a triaxial cell. As described by Zentar et al. [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF] and Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF], the pressio-triax was developed in order to study the pressuremeter loading in tightly controlled laboratory conditions. A special aspect of this equipment is its ability to measure the pore water pressure evolution at the cavity wall during testing.

On the basis of the test results obtained with this apparatus, Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF] have developed a methodology for identifying mechanical and hydraulic parameters of clayey materials, using an elastoplastic model, the Modified Cam Clay model. In this study, we extended this methodology in order to take into account the determination of the viscosity parameters of the EVP-MCC. The tests were performed on Saint-Herblain clay, a river clayey alluvial deposit from the Loire Paleolithic period, characterized by a high plasticity index (I p = 42) and slightly or moderately organic. The consolidation stress state was defined by 1 = 85 kPa, 2 = 3 = 68 kPa, u = 56.5 kPa for samples taken at a depth of 5.5-6.5 m, and 1 = 110 kPa, 2 = 3 = 90 kPa, u = 68 kPa for samples at a depth of 6.5-7.5 m. The physical and mechanical characteristics for each test are summarized in Table I, which includes the strain rate values of the performed pressio-triax tests.

An oedometer test at a constant strain rate of 3.3×10 -6 s -1 was carried out on a sample of the Saint-Herblain clay at a depth of 6.9-6.95 m using a modified oedometer apparatus [START_REF] Rangeard | Identification des caractéristiques hydro-mécaniques d'une argile par analyse inverse des essais pressiométriques[END_REF]. Lateral stress could also be measured during the test. The parameters deduced from the one-dimensional consolidation curves are presented in Table II. The soil is slightly overconsolidated and K 0 = 0.55. ra,0 , radial preconsolidation pressure and v,0 , vertical preconsolidation pressure.

Numerical modeling of pressuremeter test

Numerical modelling of pressuremeter tests was performed by using the proposed elastoviscoplastic model coupling the Biot consolidation theory in order to predict the stress evolution as well as the pore pressure change caused by the expansion of the cavity wall. During this expansion test, the soil was considered saturated and the permeability was kept constant. The initial stress state was taken to be equal to the in situ state described in the previous section.

The real geometry of the specimen was taken into account by considering a plane strain condition in the vertical direction and an axisymmetric condition in the horizontal direction with a = 6.5 mm, b = 35 mm, as shown in Figure 4. The value of the height H was taken to be equal to 6.5 mm. The mesh was composed of 27 elements with 138 nodes. The elements for all tests were isoparametric with eight integration points. The proposed elements and the mesh were verified dense enough to achieve accurate numerical solutions.

As for the boundary conditions, vertical displacement was prevented, whereas horizontal displacements could develop freely. The displacement was applied in the mini-pressuremeter probe on the side AD, and the cell pressure was maintained constant on the side BC. No water flow (flux = 0) was allowed to pass through the four faces for the undrained boundary conditions.

Sensitivity to model parameters

The influence of the model parameters on the pressuremeter curve and the evolution of the pore pressure at cavity wall were examined by simulating a pressuremeter test at a strain rate of 6×10 -5 s -1 with an initial total stress 1 = 85 kPa, 2 = 3 = 68 kPa and initial pore pressure u 0 = 56.5 kPa. The reference values of the parameters for this test are given in Table III. The influence of the parameters common to the MCC and the EVP-MCC models {E, , * , M, p c0 } is in accordance with the results presented by Zentar et al. [START_REF] Zentar | Identification of soil parameters by inverse analysis[END_REF] and Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF]:

• the parameters E, M, p c0 have a significant effect on the stress-strain relationship as well as on the evolution of pore pressure; • the parameters and * have a very small effect on the evolution of the total stress and the pore pressure; • only Young's Modulus E can change the initial slope of the strain-stress curve.

Special attention was paid to the influence of the viscosity parameters (N , ). Their effects on the radial total stress and the pore pressure at the cavity wall are presented in Figure 5. Their evolutions are plotted as a function of the relative displacement at the cavity wall ra = ratio of the cavity wall displacement to the initial cavity radius. We can note that N and have an important effect on the computed stress, especially for 5<N <20 and <1×10 -6 , and a small effect on the pore pressure. The information can be used to determine these two parameters from the sole total stress curve. To obtain the same stress evolution for a given constant strain rate, N and log( ) are related by a linear relationship as shown in Figure 6(a). For other strain rates, the selected values of (N , ) produce much more changes on the stress-strain curve, as presented in Figure 6(b), where * ra represents the computed total radial stress for different values of the viscosity parameters normalized by the stress corresponding to (N = 10, = 1×10 -4 ). From this information, the values of (N , ) can be determined from pressuremeter tests at three different strain rates.

The influence of (N , ) on a strain holding stage (stress relaxation test) at a strain of 8% and on a stress holding stage (creep test) at a stress increment of 41 kPa after an initial constant strain rate loading was also studied and is presented in Figure 7. The results show the influence of (N , ) on the evolution of the radial stress during a strain holding stage or on the radial strain during a creep test. It is therefore possible to determine the values of (N , ) from either a relaxation test or a creep test. The soil permeability k has to be taken into account in pressuremeter tests. For tests in perfectly drained or undrained conditions, k has no influence on the stress and pore pressure responses, whereas for tests in partially drained condition, k has a significant effect. Figure 8 shows that the permeability has a significant effect on the pore pressure evolution for values between 10 -8 and 10 -12 m/s, and a very small effect for other values, which correspond to either perfectly drained or undrained conditions during the pressuremeter test. However, the permeability has a very small effect on the radial stress. We can point out that the influence of the viscosity parameters and of the soil permeability on the stress-strain relationship and on the pore pressure evolution is independent of each other. This information can be used to determine the soil permeability from the sole evolution of the pore pressure at cavity wall for given values of EVP-MCC parameters.

Identification procedure

Based on the previous parametric study, we firstly determined several parameters as follows:

• Poisson's ratio was taken to be equal to 0.3, which is a common value for clay. • Young's modulus, which controls the initial slope of the radial stress-radial strain curves in the elastic region, was determined as E = 3500 kPa. • The compressibility index was determined from oedometer tests: * = ( -)/(1+e 0 ) = (0.43-0.02)/(1+2.1) = 0.13.

The procedure of parameter identification by inverse analysis, as presented in Figure 3, was employed to determine the other parameters from pressio-triax tests with three different strain rates. The iterative process to update the parameters takes the following order: M → p c0 → N → m and k, as presented in Figure 9. The optimization loops (Table IV) were carried out as follows:

1. At first, an interval of acceptable values for the set of parameters {M, p c0 , N , , k} was selected, such as in the first two lines in Table IV. The initial value of each parameter was then given as the average of the upper and lower bound values. 2. The incremental loop for determining M was carried out by using, at each step, the average value of the defined interval. For each value of M, the incremental loop for the determination of p c0 was carried out by using the average value of the selected interval as for M. The interval of parameter values was re-defined by the updated parameters. In Table IV, 'UP' and 'LOW' mean selected parameter values defined as new upper and lower bounds, respectively, which result in predicted curves above and under, respectively, the experimental ones. 3. For each value of M and p c0 , the values of N , , and k could be determined by using the same process as for p c0 , and also with the re-definition of their bound values by the updated parameters. Noting the small influence of the viscosity parameters on the evolution of the pore pressure and the small effect of the permeability on the total stress evolution, the iterative process for the parameters and k was set in the same iterative loop. 4. The values of the parameters were assessed after the convergence of the calculation according to a criterion on the difference between the current step and the last step smaller than a given tolerance, as presented in Table IV for the iterating process concerning p c0 , N , , and k with M = 1. 5. For the iterating process concerning M, once the calculation converged, that is, the difference between experimental data and calculation results became sufficiently small and stable, according to a criterion that allowed us to calculate the overall error between experimental and numerical results, we then obtained the final values of the parameters. In Table IV, we give the final values of the parameters at the end of the iterative process.

Figure 10 shows the predicted results after the convergence for p c0 = 40.6 with M = 1.0. A good performance is achieved for the total stress curves, but not for the pore pressure curves. The iterative process goes to the next loop of M until a good general agreement is achieved, i.e. L n (P)< . In Figure 10 one can see the numerical simulations obtained with the final set of parameters (M = 1.25). The predicted results after the convergence of the whole iterative process are in good agreement with the whole set of experimental results. The final values of the parameters determined by the proposed procedure, in the last row of Table IV, are summarized in Table V.

VERIFICATION OF THE IDENTIFICATION PROCEDURE

The value of the slope of critical state line M = 1.25 corresponds to a friction angle equal to 31 • , which agrees well with the results of triaxial tests in compression performed on the same samples (Figure 11(a)) reported by Zentar [START_REF] Zentar | Analyse inverse des essais pressiometrique, application a l'argile de Saint-Herblain[END_REF]. The preconsolidation pressure p c0 = 30.1 kPa corresponds to an overconsolidated ratio OCR = 1.1 ( v = 33.3 kPa calculated at the depth of 6 m, K 0 = 0.55 as measured by Zentar [START_REF] Zentar | Analyse inverse des essais pressiometrique, application a l'argile de Saint-Herblain[END_REF] and Rangeard [START_REF] Rangeard | Identification des caractéristiques hydro-mécaniques d'une argile par analyse inverse des essais pressiométriques[END_REF]). Owing to the strain-rate-dependency of the soil strength, it appears reasonable that the optimized preconsolidation pressure was found smaller than the measured one ( p c0 = 37-50 kPa by Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF]). The value of the soil permeability k = 1.78×10 -11 m/s is smaller than the one measured from conventional oedometer tests (Figure 11(b)) by Yin [START_REF] Yin | Modelisation viscoplastique des argiles naturelles et application au calcul de remblais sur sols compressibles[END_REF]. It could be due to the fact that the permeability value was kept constant during the numerical simulation, not taking into account the change of void ratio during consolidation. Furthermore, in order to examine the coherence of the set of parameters obtained by the identification procedure, a pressuremeter test at a constant strain rate of 6×10 -7 s -1 (Pre3B) and a constant strain rate test (6×10 -5 s -1 ) with a relaxation stage (Pre3H) were simulated using the determined parameter values. The comparison between predicted and measured results, presented in Figures 12 and13, showed that the EVP-MCC is able to describe the time-dependent behaviour observed in pressuremeter tests and that the proposed method can be successfully applied for identifying the soil parameters. The time for dissipating the excess pore water pressure is identical for both numerical and experimental results during the relaxation stage, although the value of the predicted pore pressure is slightly smaller than the measured one. Moreover, the determined set of parameters was taken for simulating a pressuremeter creep test (Pre1A), except for the value of the preconsolidation pressure p c0 = 35 kPa due to a different initial state of consolidation for a 1 m deeper sample. A good agreement between predicted and measured results was achieved, as presented in Figure 14. These different cases show the ability of the proposed procedure to identify correctly the soil parameters by means of pressuremeter tests with different strain rates. The identification procedure of the EVP-MCC model parameters was also tested by modelling an undrained triaxial test with three levels of constant strain rates (Tri1B), an undrained triaxial creep test with two levels of deviatoric stress (Tri1F), an oedometer test with different strain rates (Oed1), and an oedometer test under constant loading (Oed2) performed on the same clay sample, as presented in Figures 15 and16. The value of Young's modulus (E = 5000 kPa) for the simulation of the triaxial and oedometer tests was taken higher than the one determined from pressuremeter tests due to a different rigidity between vertical and horizontal directions of this anisotropic soil, while the proposed model does not take into account the anisotropy feature of natural soils. Besides this change, all the other parameters were kept at the same values as the ones determined by the identification procedure from pressuremeter tests. We can see in Figures at failure, called tertiary creep, but predicts instead a constant creep rate when failure is reached. Note that the non-homogeneity of the tested samples and the initial anisotropy feature of natural clay can result in difficulties when several different tests on different specimens are concerned. The identified value of the preconsolidation pressure p c0 is found to agree well with the results obtained from modelling the oedometer test at a constant strain rate.

ANALYSIS OF DRAINAGE CONDITIONS

The phenomenon of local consolidation can take place during a pressuremeter test, depending on the strain or stress rate at the cavity wall as well as on the soil permeability. Based on the numerical tests at constant strain rates varying from 6×10 -2 to 6×10 -7 s -1 , an analysis of the evolution of the total stress ( ra -u 0 ) and of the effective stress ( ra -u) at a strain of 8% as a function of (d a /dt)/k was made using the EVP-MCC parameters determined in Table V. Four different types of behaviour during a pressuremeter test could be observed, as already described by Rangeard et al. [START_REF] Rangeard | Determining soil permeability from pressuremeter tests[END_REF] using the Modified Cam Clay model (Figure 17):

• Fully drained. No excess pore pressure during the pressuremeter test.

• Partially drained type A. The total pressure and the pore pressure are influenced by drainage condition. • Partially drained type B. The pore pressure is influenced by drainage condition but not the total pressure. • Fully undrained. The pore pressure evolution remains the same, whereas the total pressure increases with strain rate.

To analyse the drainage conditions for given values of permeability and strain rate, we can transform Figure 17 into the one presented in Figure 18, as described by Rangeard [START_REF] Rangeard | Identification des caractéristiques hydro-mécaniques d'une argile par analyse inverse des essais pressiométriques[END_REF]. The points representing each test allow us to determine their drainage conditions. From Figure 18, we observe that Test Pre3B, with the lowest strain rate, is close to fully drained condition. Test Pre3F, with an intermediate strain rate, is in partially drained type A condition. For Tests Pre3A and 3D, both in partially drained type B condition, the difference in total pressure for different strain rates is due to the viscous behaviour of the soil skeleton. The strain rate effects in pressuremeter tests on Saint-Herblain clay were studied based on the determined soil parameters. The total stress and pore pressure at the cavity wall, at a strain of 8%, as a function of the logarithm of the strain rate are presented in Figure 19. The prediction for the strain rate varying from 1.5×10 -4 to 6×10 -7 s -1 agrees with experimental results. The total stress and the pore pressure increase with the strain rate due to the viscous properties as well as due to the different drainage conditions. As for other strain rates, we note that:

• For a strain rate higher than 6×10 -4 s -1 , the consolidation effect is erased due to fully undrained condition. The total stress or the effective stress increases linearly with the logarithm of the strain rate due to the viscous properties only, as described by Vaid and Campanella [START_REF] Vaid | Time-dependent behavior of undisturbed clay[END_REF] for undrained triaxial tests. • For a strain rate smaller than 6×10 -8 s -1 , no pore pressure is generated due to fully drained condition. The total stress or the effective stress decreases down to a constant value with the decrease in the strain rate only because of the viscous properties. This effect decreases for low strain rates, and therefore strain rates smaller than 1×10 -8 s -1 on Saint-Herblain clay can erase either viscous or consolidation effect.

APPLICATION TO IN SITU PRESSUREMETER TESTS

SBP tests were performed in the Burswood Peninsula site at a depth of 5.25 m by Lee and Fahey [START_REF] Goh | Application of a 1-dimensional cavity expansion model to pressuremeter and piezocone tests in clay[END_REF]. The expansion phase was conducted at two different rates to around 10% cavity strain: 0.167 and 0.0185%/s. The parameters from oedometer and triaxial tests by Lee and Fahey [START_REF] Goh | Application of a 1-dimensional cavity expansion model to pressuremeter and piezocone tests in clay[END_REF] with viscosity parameters determined by the proposed procedure of parameter identification are summarized in Table VI. Figure 20 compares the predicted and measured results for SBP tests in terms of total pressure at the cavity wall versus cavity strain. The proposed methodology by inverse analysis is found efficient for identifying soil parameters from in situ testing.

CONCLUSION

The possibility of determining viscous properties from in situ testing has been discussed. An inverse analysis method for identifying soil parameters from constant strain rate pressuremeter tests has been proposed. This method of identification was employed to determine the parameters of an elasto-viscoplastic model, the EVP-MCC model, as well as the soil permeability. The procedure was applied to a soft clay, the Saint-Herblain clay. The procedure was then verified by comparing the optimized values of the soil parameters and the values of the same parameters obtained from laboratory test results. Additional validation was also performed by using the identified parameters to simulate other tests on the same samples, such as pressuremeter tests at a constant strain rate and with a relaxation or a creep stage, triaxial and oedometer tests. The values of the parameters obtained from pressuremeter tests agreed with those obtained from triaxial and oedometer tests. This showed that the soil parameters could be adequately identified by pressuremeter tests with three different strain rates.

A drainage condition analysis has then been performed by observing the evolution of the total and the effective pressure at the cavity wall as a function of the strain rate and the permeability. Four different types of drainage conditions have been defined by the value of the ratio between strain rate and permeability. The drainage conditions of the pressuremeter tests performed on Saint-Herblain clay were then determined. The results show that the strain rate effects are due to both the viscous properties and the consolidation condition during a pressuremeter test in field.

All the numerical simulations have shown that the proposed methodology is efficient to be used for identifying soil parameters, and the EVP-MCC model could adequately describe the timedependent behaviour of soft soils under different monotonic loading conditions. Further validation for the parameter identification is needed, such as the prediction of the behaviour of infrastructures built on soft clay by using the recorded deformations during the construction.

NOTATIONS

Figure 1 .

 1 Figure 1. Schematic behaviour of the EVP-MCC model during CAU triaxial compression and triaxial creep tests.

Figure 2 .Figure 3 .

 23 Figure 2. Resolution of the (a) direct problem and (b) inverse problem.

Figure 4 .

 4 Figure 4. Finite element geometry for pressuremeter tests.

Figure 5 .

 5 Figure 5. Effect of viscosity parameters on the total radial stress and the pore pressure for pressuremeter test.

Figure 6 .

 6 Figure 6. Coupling effect of viscosity parameters on the total radial stress for pressuremeter test.

Figure 7 .

 7 Figure 7. (a) A comparison of stress relaxation for different (N , ) during strain holding stage and (b) a comparison of strain evolution for different (N , ) during stress holding stage.

Figure 8 .

 8 Figure 8. Permeability effects on radial stress and pore pressure responses for pressuremeter tests.

Figure 9 .

 9 Figure 9. Iterative process to identify soil parameters in the order of M → p c0 → N → m and k.

Figure 10 .

 10 Figure 10. (a) Results for different p c0 with M = 1.0 and (b) results for different values of M.

Figure 11 .Figure 12 .

 1112 Figure 11. (a) Results in q-p space of the triaxial tests in compression and (b) soil permeability versus voids ratio by conventional oedometer tests (H, sample in the horizontal direction and V, sample in the vertical direction).

Figure 13 .

 13 Figure 13. A comparison between experiment and simulation of a pressuremeter strain rate test with a relaxation stage (Pre3H).

Figure 14 .

 14 Figure 14. A comparison between experiment and simulation of a pressuremeter creep test (Pre1A).

Figure 15 .

 15 Figure 15. A comparison between experiment and simulation of triaxial tests (Tri1B, Tri1F).

Figure 16 .

 16 Figure 16. A comparison between experiment and simulation of oedometer tests (Oed1, Oed2).

Figure 17 .

 17 Figure 17. Evolution of ( ra -u 0 ) and ( ra -u) as a function of (d a /dt)/k.

Figure 18 .

 18 Figure 18. Drainage conditions as a function of strain rate and permeability.

Figure 19 .

 19 Figure 19. Strain rate effects during pressuremeter tests on Saint-Herblain clay.

Figure 20 .

 20 Figure 20. A comparison between predicted and experimental results for SBP tests on Burswood clay.

Table I .

 I Physical and mechanical characteristics of Saint-Herblain clay samples.

	Specimen	Depth (m)	w (%)	e i	(kN/m 3 )	Strain rate (s -1 )
	Pre1A (CR)	6.5-7.5	113	2.69	14.91	3×10 -5
	Pre3A (SR)	5.5-6.5	97	2.23	13.96	3×10 -5
	Pre3B (SR)	5.5-6.5	96	2.48	14.21	6×10 -7
	Pre3D (SR)	5.5-6.5	85	2.16	14.22	1.5×10 -4
	Pre3F (SR)	5.5-6.5	90	2.51	14.41	6×10 -6
	Pre3H (RE)	5.5-6.5	93	2.29	14.35	6×10 -5
	Tri1B (SR)	5.5-6.5	89	2.32	14.76	2.8×10 -5 to -7
	Tri1F (CR)	5.5-6.5	86	2.84	14.87	-
	Oed1	6.5-7.5	87	2.26	14.85	3.3×10 -4 , 6.6×10 -5
	Oed2	5.5-6.5	93	2.61	14.88	-

CR, creep test; SR, strain rate test; and RE, relaxation test.

Table II .

 II Parameters deduced from oedometer test.

	Depth (m) w (%)	e i	(KN/m 3 )			e 0	ra,0 (kPa)	v,0 (kPa)
	6.9-6.95	87	2.26	14.85	0.43	0.02	2.1	29	52

Table III .

 III Reference values of EVP-MCC model and hydraulic parameters for natural Saint-Herblain clay.

				M	p c0 (kPa)	N	(s -1 kPa -1 )	k (m/s)
	3500	0.3	0.16	1.25	30	10	1×10 -9	1×10 -11

E (kPa) *

Table IV .

 IV Optimization loops for identifying soil parameters by inverse analysis.

	Loop of M	Loop of p c0	Loop of N		Loop of , k
	M		p c0 (kPa)	N		Log( )(s -1 kPa -1 )	Log(k) (m/s)
	0.5	Low	25	Low	15	Low	-6	L o w	-7	L o w
	1.5	Up	50	Up	5	Up	-15	Up	-15	Up
							-10.5	U p	-11	Low
							-8.25	Low	-13	Low
			37.5		10		-9.375	Low	-14	Low
							-9.938	Up	-14.5	L o w
							-9.656	Up	-14.75	Low
			37.5	Low	10	Converge -9.516 Converge -14.875 Converge
	1						-10.5 -8.25	U p Low	-11 -13	Low Low
							-9.375	Up	-14	Low
					10		-8.813	Low	-14.5	L o w
							-9.094	Up	-14.75	Low
							-8.953	Low	-14.875	Low
			43.75		10	Low	-9.023 -9.059 Converge -14.969 Converge Low -14.938 Low
							-10.5	U p	-11	Low
							-8.25	Up	-14	Low
					7.5		-7.125 -7.688	Low Low	-14.5 -14.75	L o w Low
							-7.969	Up	-14.875	Low
							-7.828	Up	-14.938	Low
			43.75	Up	7.5	Converge -7.758 Converge -14.969 Converge
							-10.5	U p	-11	Low
							-8.25	Low	-14	Low
			40.625		10		-9.375 -8.813	Up Low	-14.5 -14.75	L o w Low
							-9.094	Up	-14.875	Low
							-8.953	Low	-14.938	Low
	1	L o w	4 0 .625 Converge 10	Converge -9.023 Converge -14.969 Converge
			37.5	U p	7 .5	( F i n a l )	-7.828	(Final) -10.75	(Final)
			31.25	Up	9.375	(Final)	-8.813	(Final) -10.75	(Final)
	1.25		28.125	Low	11.25	(Final)	-9.375	(Final) -10.75	(Final)
			29.688	Low	10	(Final)	-9.094	(Final) -10.75	(Final)
			30.469	Up	10	(Final)	-9.023	(Final) -10.75	(Final)
	1.25	Converge 30.078 Converge 10	(Final)	-8.988	(Final) -10.75	(Final)
	up, upper bound value; low, lower bound value; and final, after convergence.

Table V .

 V Values of EVP-MCC model and hydraulic parameters of Saint-Herblain clay.

	E (kPa)		*	M	p c0 (kPa)	N	(s -1 kPa -1 )	k (m/s)
	3500	0.3	0.13	1.25	30.1	10	1.03×10 -9	1.78×10 -11

  [START_REF] Rangeard | Identification des caractéristiques hydro-mécaniques d'une argile par analyse inverse des essais pressiométriques[END_REF] and 16 that a satisfactory agreement was achieved between numerical and experimental results on multiple stages strain rate tests, whereas a global agreement was achieved for creep tests. For these last tests, one should precise that the model cannot reproduce a phase of accelerated creep rate

	100									100
		75	1%/h Experiment EVP model	10%/h	0.1%/h		80	Experiment EVP model
	q, δu (kPa)	50		0.1%/h					q (kPa)	40 60		M=1.25
		25
					Tri1B Initial stress state:		20
					σ'	v0 =87 kPa, σ' h0 =77 kPa
		0									0
		0	1	2	3		4	5	6	7	0	20	40	60	80	100
					ε a	(%)						p' (kPa)
		2									25
				EXP. q=34 kPa
				EXP. q=39 kPa
	1.5		EVP model					18.75
			Tri1F Initial stress state:
	ε a (%)	1			σ' v0 =64 kPa, σ' h0 =35 kPa			δu (kPa)	12.5
	0.5									6.25			EXP. q=34 kPa
														EXP. q=39 kPa
														EVP model
		0									0
		10 0	10 1		10 2		10 3	10 4	10 5		1 10	4	2 10	4	3 10	4
					t (min)						t (min)

Table VI .

 VI Values of EVP-MCC model and hydraulic parameters for Burswood clay.

	Sites	E (kPa)			*	M	p c0 (kPa)	N	(s -1 kPa -1 )	k (m/s)
	Burswood	4000	0.3	0.07	1.45	38		12	2×10 -7	3.3×10 -9
				160
				140
			(kPa)	120
			σ ra	100
								Exp. 0.167%/s
				80				Exp. 0.0185%/s EVP model
				60
				0	2	4	6	8	1 0	1 2
						δ	ra (%)

ACKNOWLEDGEMENTS

The research reported in this paper is funded by a grant from 'Conseil Regional des Pays de la Loire'. The authors thank Dr Christophe Dano and Dr Yvon Riou, associate professors at Research Institute in Civil and Mechanical Engineering of Ecole Centrale de Nantes, and professor Hong-Wei Huang in Tongji University for their support.


L Filice

I Alfaro

F Gagliardi

E Cueto

F Micari

F Chinesta

A preliminary comparison between finite element and meshless simulations of extrusion

Keywords:

In this paper the extrusion process of a cross-shaped profile was investigated. In particular, the study was focused on the distortion of extruding profiles when the workpiece and die axis are not aligned. The process was simulated using the finite element method (FEM) and the natural element method (NEM), both implemented in an updated-Lagrangian formulation, in order to avoid the burden associated with the description of free surfaces in ALE or Eulerian formulations. Furthermore, an experimental equipment was developed in order to obtain reliable data in terms of deformed entity, required process load and calculated pressure. At the end, a comparison between the numerical predictions and the experimentally measured data was carried out. The main results are presented in the paper.

Introduction

Extrusion is one of the most utilised bulk metal forming processes. It is particularly suitable when high reductions in area are required or when hollow components have to be manufactured. Extrusion analysis has been carried out, in the past, utilising some analytical techniques based on the slab method although such approach results applicable only if simple geometries are taken into account. The introduction of finite elements and, in particular, the development of efficient mesh management procedures, allowed the simulation of very complex processes characterised by large surfaces generation, high strains and strain rates.

Up to now, the most relevant limitation is associated with the numerical dissipation introduced with remeshing procedures. When the mesh becomes too distorted in a particular domain, a new mesh must be generated and the variables associated with material history must be mapped from the old mesh to the new. This introduces additional errors that meshless formulation can avoid. On the other hand, the simulation of hollow components extrusion introduces another critical aspect since material joining (after the previous separation in correspondence of the porthole) has to be properly modelled. According to the above considerations, the development of a new modelling technique able to introduce a sort of discontinuity with respect to the currently applied ones can be strategic.

Very few examples exist of the application of meshless methods to the simulation of forming processes, and none (up to our knowledge) to extrusion. Nevertheless, we can cite the pioneering work by [START_REF] Chen | Analysis of metal forming process based on meshless method[END_REF] as one of the earliest examples, together with [START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF]. They use, however, early implementation of meshless methods that still have some drawbacks. One of the most important is the lack of proper interpolation on the boundary. This produces some problems in the correct imposition of essential boundary conditions and, of course, in the appropriate simulation of contact and friction.

In this paper a well-established meshless technique known as natural element method (NEM, [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] is adapted to the simulation of extrusion processes. The NE technique presents some advantages over finite element simulations: no remeshing is necessary since a satisfactory accuracy of the approximation is achieved even with highly distorted triangles or tetrahedra. The price to pay is, obviously, a higher computational cost. The method is here employed in conjunction with the notion of ˛-shapes of the cloud of nodes in order to extract the geometry of the extruded part as it evolves. This state-of-the-art geometrical concept allows to avoid complex geometrical checks of self-contact on the boundary of the domain, as shown by Alfaro et al. (2006a) or Alfaro et al. (2006b).

In the study here addressed a simple 3D extrusion process was taken into account, characterised by an asymmetric die which determines a distortion of the extruded profile. The study of profile distortion is important for practical purposes of real-life extrusions, but also to test the geometrical algorithms employed to track the free surfaces in the simulation. As far as the proposed approach is concerned, home-made FE and NE codes were utilised for the process modelling. In order to analyse the origins of possible discrepancies between both techniques, the heat generation due to plastic deformation has been neglected, although it can be considered without any limitations to the technique, as in some of our previous works, see for instance Alfaro et al. (2006a) or Alfaro et al. (2006b). Finally, a proper experimental equipment was developed. The numerical results were compared to the experimental evidences taking into account some interesting process variables, namely the required process load, the pressure on the die and the geometry of the extruded profile.

2.

The case under study

Constitutive modelling

The process geometrical data are reported in Fig. 1. The chosen material for all the tests is UNS L51120 chemical lead, due to its property to yield at room temperature with moderate loads. Its behaviour was modelled through a power-law type constitutive equation, reported in

= 60ε 0.05 (MPa) (1)
More in detail, three sets of upsetting tests were carried out on cylindrical specimens, reducing as much as possible the friction between the basis and the press-plates. The tests were carried out utilising an Instron 8501 electronically controlled hydraulic testing machine and the average strain rate was fixed as 0.05 (hence, quasi-static), 50 and 100 s -1 in order to range into a wide area. Results are depicted in Fig. 2. The flow rule was derived by using an inverse approach based on the predicted data of a simple numerical simulation. To this end, a general viscoplastic Norton-Hoff law was adjusted with the experimental data.

The general Norton-Hoff viscoplastic model can be stated as

dev = 2 3 d 0 ( √ 3) n+1 e ˇT dn εh d (2)
where the superscript dev refers to the deviatoric part of the stress tensor, d represents the equivalent strain rate, and n, h and ˇare parameters of the model. We have neglected the dependence of the model with temperature and strain, assuming dependence on strain rate only. After particularising this model for the uniaxial case, we arrive to 0 = 33.70 and n = 0.05. We thus obtain a rigid viscoplastic model (similar to the flow formulation by [START_REF] Zienkiewicz | Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes[END_REF] in the form:

dev = 2 d, (3)
where

= 2 3 0 ( √ 3) n+1 dn-1 .
The requested process load and the shape of the obtained extruded profiles were measured. In particular, a configuration with an asymmetric die with 10 mm of misalignment (Figs. 1 and3) was investigated in order to induce the distortion of the extruded profile. The geometry of the profile was obtained after [START_REF] Zhou | 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach[END_REF].

Constitutive modelling of friction

In the simulations hare presented, a viscoplastic model for the friction between the lead and the walls of the die was considered. Following the notation established by [START_REF] Chenot | Numerical treatment of contact and friction in Fe simulation of forming processes[END_REF], the viscoplastic friction law is described by a non-linear constitutive equation, relating the shear stress vector and the relative velocity between the lead (in this case) and the walls of the die, v s :

= -˛f( n )K|v s | q-1 v s (4)
where n represents the normal stress at the interface. ˛f, the viscoplastic friction coefficient, can be established as a function of the normal stress n [START_REF] Chenot | Numerical treatment of contact and friction in Fe simulation of forming processes[END_REF] and the sensitivity parameter q is very often taken equal to m, and thus equal to 0.05.

Governing equations

We consider the balance of momentum equations, without inertia and mass terms:

∇ = 0, (5)
and the assumed incompressibility of a von Mises-like flow:

∇v = 0, ( 6
)
where v represents the velocity field. The stress-strain rate relationship is given by Eq. ( 2).

Together with these equations, appropriate boundary conditions are considered:

n = t in t (7) v = v in v ( 8
)
where t and v represent, respectively, the part of the boundary = ∂ where tractions and velocities are prescribed.

If we write the incremental variational equation at time t + t we arrive to

(t+ t) (-(p t + p)I + 2 (d t + d)(d t + d)) : d * d = 0. ( 9
)
Domain updating is done in an explicit procedure, given the last converged velocity field, but due to the non-linear character of the constitutive equations, an iterative approach has been applied to the conservation equations, using the Newton-Raphson scheme, thus leading to

(t+ t) - pI + 2 ∂ (d t+ t k ) ∂d : d d t+ t k +2 (d t+ t k ) d : d * d = - (t+ t) (-p t+ t k I + 2 (d t+ t k )d t+ t k ) : d * d (10)
where the subscript k indicates the iteration within a time increment. The incremental form of the incompressibility condition results:

(t+ t) ∇( v) p * d = - (t+ t) ∇(v t+ t k )p * d (11)
If we approximate the velocities and pressures, as well as their variations, by employing a finite-dimensional set of basis functions, we arrive to a discrete form of the previous equations (Bubnov-Galerkin method):

v h (x) = n I=1 I (x) v I (12) p h (x) = n I=1 I (x) p I . ( 13
)
where n represents the number of nodes considered in the approximation. The functions I (x) and I (x) in this work represent some form of finite element (piece-wise polynomials) or natural neighbour interpolation, which will be studied in the following section. This leads to the linear system:

K t+ t k G G T 0 v p = f t+ t k 0 (14)
being the updating given by

v t+ t k+1 = v t+ t k + v (15)
p t+ t k+1 = p t+ t k + p (16)
Obviously, Eqs. ( 10) and ( 11) will be solved iteratively until convergence. The convergence is assumed once the relative modulus of the residual gets under 10 -6 and the relative increment of the velocity norm is below 0.02.

Numerical modelling of the extrusion process

Finite element analysis

From the FE point of view, a linear velocity-constant pressure mixed formulation, based on four node tetrahedra, was employed. This simple formulation possesses obvious limitations (noteworthy some tendency to lock) but presents the same structure as the implemented NE formulation, and the same degrees of freedom per node. Thus it is in the opinion of the authors the most appropriate in order to make a comparison with the NEM results. Despite the lack of compliance with the LBB condition, no locking was observed during the simulations. It should be highlighted, in addition, that the FE formulation employed here is by no means standard in the sense that the FE nodes remain the same throughout the simulation. In this way, no projection of variables is made. This is not frequent in commercial FE codes. At each time step a Delaunay triangulation is performed in the FE nodes, even if it is well-known that the Delaunay triangulation is mainly used in two-dimensional applications.

Natural neighbour approximations

We review here with some more details about the basics of the natural element method, since it is a technique much less known than the finite element method. The vast majority of meshless methods are based on the employ of scattered data approximation techniques to construct the approximating spaces of the Galerkin method. These techniques must have, of course, low sensitivity to mesh distortion, as opposed to FE methods. Among these techniques, the natural element method employs any instance of natural neighbour interpolation [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF][START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF] to construct trial and test functions. Prior to the introduction of these interpolation techniques, it is necessary to define some basic concepts. The model will be constructed upon a cloud of points with no connectivity on it. We will call this cloud of points N = {n 1 , n 2 , . . . , n M } ⊂ R d , and there is an unique decomposition of the space into regions such that each point within these regions is closer to the node to which the region is associated than to any other in the cloud. This kind of space decomposition is called a Voronoi diagram of the cloud of points and each Voronoi cell is formally defined as (see Fig. 4):

T I = {x ∈ R d : d(x, x I ) < d(x, x J )∀J / = I}, (17)
where d(•, •) is the Euclidean distance function.

The dual structure of the Voronoi diagram is the Delaunay triangulation, obtained by connecting nodes that share a common (d -1)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there being some so-called degenerate cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in two-dimension, as depicted in Fig. 4, right). Another way to define the Delaunay triangulation of a set of nodes is by invoking the empty circumcircle property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called natural neighbours and hence the name of the technique.

In order to define the natural neighbour coordinates it is necessary to introduce some additional concepts. The secondorder Voronoi diagram of the cloud is defined as

T IJ = {x ∈ R d : d(x, x I ) < d(x, x J ) < d(x, x K )∀J / = I / = K}. ( 18
)
The simplest of the natural neighbour-based interpolants is the so-called Thiessen's interpolant [START_REF] Thiessen | Precipitation averages for large areas[END_REF]. Its inter- polating functions are defined as

I (x) = 1 ifx ∈ T I 0 elsewhere. ( 19
)
The Thiessen interpolant is a piece-wise constant function, defined over each Voronoi cell. It defines a method of interpolation often referred to as nearest neighbour interpolation, since a point is given a value defined by its nearest neighbour. Although it is obviously not valid for the solution of second-order partial differential equations, it can be used to interpolate the pressure in formulations arising from Hellinger-Reissner-like mixed variational principles such as the one here employed. [START_REF] Sibson | A vector identity for the Dirichlet tessellation[END_REF] defined the natural neighbour coordinates of a point x with respect to one of its neighbours I as the ratio of the cell T I that is transferred to T x when adding x to the initial cloud of points to the total volume of T x . In other words, if Ä(x) and Ä I (x) are the Lebesgue measures of T x and T xI , respectively, the natural neighbour coordinates of x with respect to the node I is defined as

I (x) = Ä I (x) Ä(x) . ( 20
)
In Fig. 5 the shape function associated to node 1 may be expressed as

1 (x) = A abfe A abcd . ( 21
)
It is straightforward to prove that NE shape functions (see Fig. 6) form a partition of unity [START_REF] Babuška | The partition of unity method[END_REF], as well as some other properties like positivity (i.e., 0 ≤ I (x) ≤ 1∀I, ∀x) and strict interpolation:

I (x J ) = ı IJ . ( 22
)
Recently, [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF] have generalised the form of natural neighbour interpolants. One different type of interpolation has attracted the interest of researchers, since it is slightly faster to compute, although gives less smooth interpolations. It has received the name of Laplace interpolant.

Consider the introduction of the point x in the cloud of nodes. Due to this introduction, the Voronoi diagram will be altered, affecting the Voronoi cells of the natural neighbours of x. The Laplace interpolant is defined by using geometrical entities of one dimension less than the original space under consideration. If we define the cell intersection t IJ = {x ∈ T I T J , J / = I} (note that t IJ may be an empty set) we can define the value:

˛J(x) = |t IJ | d(x, x J ) . ( 23
)
Thus, the point x shape function value with respect to node 4 in Fig. 7 is defined as

ns 4 (x) = ˛4(x) n J=1 ˛J(x) = s 4 (x)/h 4 (x) n J=1 [s J (x)/h J (x)] , (24)
where s J represent the length of the Voronoi segment associated to node J and n represents the number of natural neighbours of the point under consideration, x. Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see the work by [START_REF] Sukumar | Natural neighbor Galerkin methods[END_REF]. For the purposes of the work here presented, a mixed Thiessen-Laplace approximation for pressure and velocity, respectively, has been considered.

Laplace interpolants have some remarkable properties that help to construct the trial and test functional spaces of the Galerkin method (see [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF] for proofs of these properties). Besides properties like continuity and smoothness (everywhere except at the nodes for Sibson interpolants and at some other lines of zero measure for the Laplace interpolant), Sibson and Laplace interpolants posses linear completeness (i.e., exact reproduction of a linear field). A comparison of Sibson and Laplace shape functions for the same regular lattice of points is done in Fig. 8.

Sibson and Laplace interpolants can also reproduce linear functions exactly along convex boundaries. This is in sharp contrast to the vast majority of meshless methods. In addition, in [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled ˛-shapes[END_REF] and [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method[END_REF] distinct methods of imposing linear displacement fields along non-convex boundaries were developed.

If higher-order natural element approximations are needed, the approach presented in [START_REF] González | Higher-order natural element methods: towards an isogeometric meshless method[END_REF] allows for arbitrary orders of consistency and smoothness.

Functional approximation and discretisation of the problem

As mentioned before, in the approximation of Eqs. ( 12) and ( 13), linear approximation was chosen for velocities and constant pressure per elements, when FE technology was employed.

In the case of natural element approximation, velocities were approximated by making use of Laplace interpolation, whereas pressure was assumed constant per Voronoi cell by employing Thiessen interpolation. The choice of Laplace instead of Sibson interpolation is dictated by the compromise between accuracy and computational cost, as studied in Alfaro et al. (2007).

Tracking of the domain

Finally, in order to track the evolution of the free surface of the domain, shape constructors (in particular, ˛-shapes, see [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF] were employed. This allows us to proceed without the need of storing information about the boundary. At each time step, triangles belonging to the boundary are found by the method. ˛-shapes define a one-parameter family of shapes S ˛(being ˛the parameter), ranging from the "coarsest" to the "finest" Fig. 9 -Two snapshots of a two-dimensional simulation of an extrusion process. Equivalent plastic strain rate is depicted. level of detail. ˛can be seen, precisely, as a measure of this level of detail. This means that all details of size less than ˛will be ignored in the geometry of the domain. This is a common practice in many conversions from CAD models (very often geometrically "exact") to FEM models, in which very small features of the geometry are ignored from a mechanical point of view.

In order to clarify the before presented concepts, consider some examples of ˛-shapes computed from a cloud of points corresponding to the simulation of two-dimensional extrusion process. We restrict ourselves to geometrical concepts only.

Consider the extrusion example shown in Fig. 9, where the contour plot of equivalent plastic strain rate is depicted. The key idea of the method here proposed is to extract the shape of the domain at each time step by invoking the concept of ˛-shape of the cloud. The ˛parameter will be obtained by geometrical considerations. In this case the radius at the inlet of the die, for instance, seems to be the smallest level of detail up to which the domain (i.e., the billet) must be represented. In order to appropriately represent this value, the nodal distance h must be accordingly chosen.

In Fig. 10 some members of the family of ˛-shapes of the cloud of points in its final configuration (corresponding to Fig. 9b)) are depicted. In Fig. 10a) the member for ˛= 0, i.e., the cloud of points itself, is shown. Note how, as ˛is increased, the number and size of the simplexes (in this case, triangles) that belong to the shape is increasing. For ˛= 1.0 we obtain an appropriate shape for the cloud. Note, however, that this is not an exact value to be determined at each time step. Since the number of ˛-shapes is finite, there generally exists an interval of valid ˛values for a single shape. Finally, by increasing the ˛value, we arrive to the convex hull of the cloud of points (Fig. 10(f)).

Comparison between numerical and experimental predictions

Requested load during the process

The predicted punch load, as well as the experimental one, is reported in Fig. 11. Some underestimation of the punch load is obtained if friction is not taken into account, as expected. This is obvious, but still a good qualitative agreement is obtained. Improved results are obtained, however, if we consider in the simulation the viscoplastic friction model, detailed in Section 2.2. An overestimation of the punch load is obtained, however, at the initial time steps of the simulation. This has been repeatedly observed, since numerical simulations do not take into account imperfections in the material, voids, etc., that make the actual experiments to occur with lower loads at the initial stages of extrusion. After 3 mm of stroke displacement, FEM simulation (without considering friction) was not able to converge any more, due to the distortion of the mesh together with the incompressibility and the non-linear material model. Note that we are not allowing remeshing, neither in the FEM nor in the NEM simulations. This means that the usual practice of generating a new mesh, with improved quality tetrahedra, and projecting old variables to the new mesh, has not been accomplished. This is quite a fictitious procedure, but we tried to avoid the introduction of the well-known numerical diffusion associated to remeshing.

Shape of the extrudate

About 25 mm of extrudate were simulated with NEM. In order to verify the accuracy of the prediction, the deviation of the profile from the straight line is compared to that of the experiments and the FEM simulation. It is possible to verify the sufficiently good prediction of the distortion of the profile induced by the die asymmetry (see Fig. 12). Again, taking into account the friction, the obtained results fit very well to experimental data. In addition, comparing the results for 3 mm of stroke displacement (Fig. 13), FEM (without friction) and NEM results are very similar. Note however that at this displacement value FEM simulation could not achieve convergence due to mesh distortion. It is then concluded that the friction has an important effect on the determination of punch load, as expected, but not so much in the shape of the extrudate.

Pressure at the bottom of the die

Another interesting comparison was carried out with reference to the predicted pressure on the bottom die, along the symmetry plane, from point A to point B, according to the following Fig. 14. In this case, no experimental measurements were available, so only numerical predictions were compared. Fig. 15 shows the predicted pressure distribution. No friction was considered in this particular comparison in order to simplify the analysis. In this case the predicted pressure smoothness is better in the NE case. The FE results show some spurious oscillations, which are well-known in the FE community. The NEM, however, despite the simple formulation employed, does not show spurious oscillations. A more detailed discussion on the performance of mixed NEM approximations for incompressible materials can be found in [START_REF] González | Volumetric locking in natural neighbour Galerkin methods[END_REF]. Another important fact is that the FE shows some spurious pressure peaks due to the obvious distortion of the mesh. This is shown in Fig. 16.

The influence of the quality of the mesh on the results is highlighted in Fig. 17. In it, the simulation was ran by using NEM until a given time step. It was then stopped and one time step more was accomplished by means of both FEM and NEM. In this way, a direct comparison of the results and the influence of the technique on them can be achieved. Equivalent strain rate results are depicted. Note the spurious concentration of strain rate near a zone with distorted elements, which does not appear in the NEM results. In addition, overall FEM strain rate values are somewhat higher than those of the NEM simulation.

Conclusions

A first comparison between a FEM and a NEM formulation for the numerical simulation of extrusion processes was proposed in the paper, focusing the attention on the extrusion of a crossshaped component. At the end of this work some conclusions can be traced:

• Meshless methods, despite their relative lower state of development, show promising capabilities in the simulation of forming processes with large deformations. • NEM is today computationally more expensive because it needs more time in calculating the shape functions and, being a very novel approach, a relevant optimisation work remains to be done.

The computational cost of the NEM has been estimated in about four times that of the FEM in some of our previous works (Alfaro et al., 2007). Nevertheless, this cost is negligible, if compared to the cost of performing iterations in the Newton-Raphson loop. Very efficient algorithms have been developed in order to alleviate this computational cost (the natural neighbour search, for instance, see the work by Alfaro et al., 2007). At this moment, we believe that the NEM is a competitive technique in terms of CPU hours, if we deal with highly non-linear problems.

We believe that the use of meshless methods in the simulation of extrusion processes can be especially interesting for the extrusion of hollow profiles. Such processes involve a complex flow through the die, with melting fronts of material at the exit. Eulerian (fixed mesh) or ALE approaches can handle this type of flows with difficulties due to the complex geometry of the extrudate. This constitutes our current effort of research.

At this moment NEM formulation is not implemented in any commercial program even if it can supply a decisive answer to the need of simulating some specific problems, like the mentioned extrusion of hollow components. We believe that the NEM constitutes nowadays an appealing choice for the simulation of some complex phenomena, like those involving free surfaces, as mentioned before, or the need of an accurate description of evolving fronts.

Fig. 1 -

 1 Fig. 1 -Geometry of the die for the experiments (left) and the numerical simulation (right).

Fig. 2 -Fig. 3 -

 23 Fig. 2 -Evolution of punch load during the characterisation essays.

Fig. 4 -

 4 Fig. 4 -Delaunay triangulation and Voronoi diagram of a cloud of points.

Fig. 5 -

 5 Fig. 5 -Definition of the natural neighbour coordinates of a point x.

Fig. 6 -

 6 Fig. 6 -Typical Sibson function (x).

Fig. 7 -Fig. 8 -

 78 Fig. 7 -Definition of non-Sibsonian coordinates.

Fig. 10 -

 10 Fig. 10 -Some members of the family of ˛-shapes of the cloud of points used in the extrusion example. (a) S 0 (the cloud of points) (b) S 0.3 (c) S 0.5 (d) S 1.0 (e) S 1.5 and (f) S ∞ (the convex hull of the set).

Fig. 11 -

 11 Fig. 11 -Evolution of the punch load. Results with the FEM and NEM without friction, NEM with friction, and experiments. The best polynomial fit to the experimental results is also shown.

Fig. 12 -

 12 Fig. 12 -Deviation due to die misalignment. (a) Zone measured and (b) profile of the extrudate.

Fig. 13 -

 13 Fig. 13 -Exit velocity distribution: FEM (a) and NEM (b).

Fig. 14 -

 14 Fig. 14 -Pressure analysis section on the bottom die.

Fig. 15 -

 15 Fig. 15 -Pressure on the bottom die, along the A-B line.

Fig. 16 -

 16 Fig. 16 -Spurious pressure distribution obtained in the FE formulation during the simulation. FE results (a) and NE results (b).

Fig. 17 -

 17 Fig. 17 -Spurious equivalent strain rate distribution obtained in the FE formulation during the simulation. FE results (a) and NE results (b). The highlighted zone shows the anomalous values.

Acknowledgement

The authors would like to thank Mr. F. Pulice for his contribution in the development of the experimental equipment.