D Le Houédec

G Lefeuve-Mesgouez

B Picoux

SOIL-STRUCTURE INTERACTION FOR MOVING LOADS : APPLICATION TO RAILWAY TRAFFIC

Keywords:

The problem of determining the response of a soil medium under the action of moving loads has received considerable attention in the past. Work in this area has been motivated by the need to determine the vibratory motion on the ground surface and at depth caused by moving vehicles. Moreover, high-speed trains are becoming increasingly popular and freight trains increasingly heavier. Combined with this fact and the observation that Rayleigh wave speeds are slower in soft soils, we see that the study of moving loads is important for environmental and geotechnical engineers.

A lot of models are possible for the simulation of propagation phenomena. First, we can consider the case of a two-dimensional problem involving a moving harmonic strip load. The results derived by Fourier transform are valid for any frequency and load speed. For moving loads, to choose a suitable damping model, an original modified hysteretic damping is employed. Also, it is possible to investigate the transmission of vibrations due to a moving harmonic train of rectangular loads rigidly attached to an elastic multilayered ground. In this case, the previous damping model is used, combined with a dynamic stiffness matrix approach. In the spatial wavenumber domain, the dominant features encountered in a spectral decomposition corresponding to compression, shear and Rayleigh waves are illustrated. Also, results of the model at certain critical load speeds can be presented in dimensionless parameters.

For 3D problems, a three-dimensional semi-analytic model of ground vibration due to a moving harmonic rectangular load can be developed. The model consists of a half-space, an elastic layer over a rigid foundation, or a multilayered soil, excited by an uniform harmonic surface load moving with constant velocity. The displacements are obtained in a wave number domain after a double Fourier transform for the equations of motion. Actual displacements are then calculated using a FFT algorithm. Numerical results can be presented in the wave number and actual domains, showing the existence of different regimes (sub-Rayleigh and super-Rayleigh) depending on the load speed. Besides, for multilayered soils, we have to construct a global stiffness matrix using the same assembly process as in finite element analysis.

A more complete study needs to take into account the railway track. In this case, a railway track model lying on a layered ground and submitted to a moving train can be developed and the resolution method uses the formalism of Fourier transform for a semi-analytical resolution in the wave number domain. It includes all elements of the track (rails, pads, sleepers and ballast) and allows a parametric analysis of its different elements and evaluation of vertical displacements according to the speed, the weight and the composition of the trains. The handwriting of the stiffness matrix for a layered ground with the help of a fitted phase angle of Helmotz functions provides a fast numerical approach of the problem. Then Mach cones for the super-Rayleigh regime can be obtained near the track involving very high displacements on waves front. In the continuation of these numerical results, in situ measurements can be performed with a view to the validation of the model. These studies allow the construction of an available data basis for a possible development of models.

INTRODUCTION

The motivation for this presentation is to develop analytical, semi-analytical and numerical approaches to solve problems, involving ground vibrations due to moving loads in order to understand phenomena of wave propagation induced by rail and road traffic. In fact, during the last thirty years, damage due to ground vibrations caused by rail and road traffic has increased considerably, leading to environmental consequences (vibrations, noise, discomfort and environmental deterioration ...). It is thus important to study these phenomena theoretically in order to envisage ways of reducing their impact on the built environment. This is particularly true with regard to rail and road traffic which is becoming heavier and faster. For example, in the North West of France, more particularly in the region of the Somme Bay where grounds are constituded mainly by peat and soft clay, observations have revealed the presence of important displacements at the soil surface near railway tracks [START_REF] Picoux | Propagation dans le sol des vibrations émises par un véhicule ferroviaire se déplaçant à vitesse constante[END_REF]. Also noticed in Sweden for grounds constituted by soft clay [START_REF] Madshus | High speed railvay lines on soft ground : dynamic behaviour at critical train speed[END_REF], these displacements involve with upsetting amplitudes in the track during heavy and high speed train passages.

To remember briefly some previous studies, we can notice the two-dimensional problem of the steady-state motion of a static line load in an unbounded body which was first considered by Eason et al [START_REF] Eason | Generation of waves in an infinite elastic solid by variable body forces[END_REF]. Eason [4] also studied the two-dimensional steady-state problem for an uniform half-space, using Fourier transform methods. Cole and Huth [START_REF] Cole | Stresses produced in a half-plane by moving loads[END_REF] considered the same problem for a constant normal line load and obtained analytical forms for the solutions of displacements, in the subsonic, transonic and supersonic cases, using Helmholtz decomposition. These various cases are defined as :

• subsonic if the load speed is lower than the shear wave speed c 2 ,

• transonic if the load speed is included between the shear wave speed and the compressive wave speed c 1 , • supersonic if the load speed is higher than the compressive wave speed. Also, sometimes it is interesting to appeal to an other classification, such as :

• sub-Rayleigh if the load speed is lower than the Rayleigh wave speed c R ,

• super-Rayleigh if the load speed is higher than the Rayleigh wave speed.

These typical speeds are defined as :

• for the compressive wave :

ρ μ + λ = 2 c 2 1 (1)
• for the shear wave :

ρ μ = c 2 2
(2)

• for the Rayleigh wave (Viktorov formula) :

ν + ν + ≅ 1 12 . 1 87 . 0 c c 2 R (3)
where λ et μ are Lamé constants, ρ is the density and ν is the Poisson's ratio.

Georgiadis and Barber [START_REF] Georgiadis | Steady-state transonic motion of a line load over an elastic half-space : the corrected Cole-Huth solution[END_REF] corrected the Cole and Huth solution for the transonic regime. Fryba [START_REF] Fryba | Vibrations of solids and structures under moving loads[END_REF] also investigated the possibility of shocks in the ground due to a line load moving at constant speed at or beyond the Rayleigh wave speed. The method is based on a coupling using a change of variables (linked to the moving load) and a triple Fourier transform. This approach was also used by other authors, for example : Niwa and Kobayashi [START_REF] Niwa | Stresses produced in an elastic half-plane by moving loads along its surface[END_REF], Alabi [START_REF] Alabi | Ground vibration excited by surface trains[END_REF], Jones et al [START_REF] Jones | Ground vibration in the vicinity of a moving harmonic rectangular load on a half-space[END_REF]. Similarly, Keer [START_REF] Keer | Moving and simultaneously fluctuating loads on an elastic half-plane[END_REF] extended the study of the subsonic case to a moving harmonic line load with some numerical solutions. More recently, Gunaratne and Sanders [START_REF] Gunaratne | Response of a layered elastic medium to a moving strip load[END_REF] also studied the response of a layered elastic medium to a moving strip using a layer stiffness approach and viscous damping. Also, Lefeuve-Mesgouez et al [START_REF] Lefeuve-Mesgouez | Ground vibration due to a high -speed moving harmonic load[END_REF] have obtained numerical results for the vertical displacements in the case of a vertical harmonic strip load moving with speed beyond the Rayleigh wave speed of the ground. The transmission of vibrations over the surface of the ground, due a high-speed moving, vertical harmonic strip load, is also investigated by the same authors [START_REF] Lefeuve-Mesgouez | Ground vibration in the vicinity of a high-speed moving harmonic strip load[END_REF], in the case of a two-dimensional problem with the ground modelled as an elastic half-space. Numerical results for the displacements on the surface are obtained for loads moving with speeds up to and beyond the Rayleigh wave speed of the half-space. De Barros and Luco [START_REF] Barros | Stresses and displacements in a layered half-space for a moving line load[END_REF] investigated the response of a layered viscoelastic half-space to a moving constant line load, and studied the stresses and strains. The transient motion for a line load which is suddenly applied on the ground and then moves with constant speed has been studied by Payton [START_REF] Payton | Transient motion of an elastic half-space due to a moving surface line load[END_REF].

Eason et al [START_REF] Eason | Generation of waves in an infinite elastic solid by variable body forces[END_REF] also studied the three-dimensional problem of a point load moving in an unbounded solid. Eason [START_REF] Eason | The stress produced in a semi-infinite solid by moving surface force[END_REF] and Lansing [START_REF] Lansing | The displacements in an elastic half-space due to a moving concentrated normal load[END_REF] investigated the three-dimensional problem for an uniform half-space. Alabi [START_REF] Alabi | A model for the problem of ground vibration induced by the wheels of a moving train[END_REF], [START_REF] Alabi | A parametric study on some aspects of ground-borne vibrations due to rail traffic[END_REF] represented a train by a set of vertical moving point loads, and also studied the case of an oblique load for speeds up to half the Rayleigh wave speed. The three dimensional problem of a point load moving at constant speed on the ground demonstrating the possibility of shock formation in the ground has been studied by a lot of authors, for example Fryba [START_REF] Fryba | Vibrations of solids and structures under moving loads[END_REF], De Barros and Luco [START_REF] Barros | Response of a layered viscoelastic half-space to a moving point load[END_REF], Krylov [START_REF] Krylov | Generation of ground vibrations by superfast trains[END_REF], ... . Ground vibration in the vicinity of a moving harmonic rectangular load on a half-space can be also deduced (Jones et al [START_REF] Jones | Ground vibration in the vicinity of a moving harmonic rectangular load on a half-space[END_REF]). For the same loading, displacements can be found in the case of a multilayered soil or of an elastic layer over a rigid foundation (Lefeuve-Mesgouez et al [START_REF] Lefeuve-Mesgouez | Ground vibration due to a high-speed moving harmonic load[END_REF]. For the same characteristics of ground, it is also possible to obtain the surface vibration due to a sequence of high speed moving harmonic rectangular loads (Lefeuve-Mesgouez et al [START_REF] Lefeuve-Mesgouez | Surface vibration due to a sequence of high speed moving harmonic rectangular loads[END_REF]). This brief presentation don't give a complete state of art about the general problem concerning with the ground response for moving loads. Nevertheless, this point of view brings useful informations allowing a correct approach for the next situations :

• two or three dimensional analysis ;

• models of ground : half-space, elastic or viscoelastic layer over a rigid foundation, multilayered soil ; • models of loadings : line load, strip load or rectangular load, sequence of loads.

In the objective to present the more important conclusions deduced from these papers, in the next section we will investigate theoretically some problems for the ground vibration in the vicinity of a high-speed moving harmonic load.

GROUND RESPONSE FOR HIGH-SPEED MOVING LOADS

For this presentation, we consider only the case of a strip load (two-dimensional problem) and a rectangular load (three-dimensional problem) moving with a constant speed. The model considered is shown in Figure 1. The strip load has a width of 2b and is aligned with respect to the x 2 -axis. It rests on an homogeneous, isotropic, half-space, with material properties E (Young's modulus), ρ (density), ν (Poisson's ratio). We consider the soil exhibiting hysteretic damping (coefficient η) which yields complex body wave speeds as follows :

) i 1 ( 2 c 2 1 η + ρ μ + λ = ( 4
) ) i 1 ( c 2 2 η + ρ μ = (5)
But for load speeds above the Rayleigh wave speed, yielding a quasi-causal solution, we need to modify the definition of the loss factor to :

[ ] ) c ( sgn i 1 2 c 2 1 β - ω η + ρ μ + λ = (6) [ ] ) c ( sgn i 1 c 2 2 β - ω η + ρ μ = ( 7
)
where ω is the load pulsation, c the load speed and β the Fourier parameter relative to the moving space variable x.

An harmonic vertical load acts uniformly over the strip.

Fig. 1 -Geometry of model By introducing a Helmholtz decomposition in the Navier's elastodynamic equations and the stress-strain relatinship, we obtain the following wave equations for the two-dimensional model :

0 t c 1 x x 2 2 2 1 2 3 2 2 1 2 = ∂ Φ ∂ - Φ ∂ ∂ - ∂ ∂ (8) 0 t c 1 x x 2 2 2 2 2 3 2 2 1 2 = ∂ Ψ ∂ - Ψ ∂ ∂ - ∂ ∂ (9)
with :

1 3 3 3 1 1 x x u , x x u ∂ Ψ ∂ + ∂ Φ ∂ = ∂ Ψ ∂ - ∂ Φ ∂ = ( 10
)
where Φ and Ψ are the Helmholtz potentials, the co-ordinates (x 1 , x 3 ) are the fixed frame of reference, u 1 and u 3 are the horizontal and vertical displacements respectively. In terms of the scalar potentials, the boundary conditions for this problem are as follows :

• at x 3 = 0 < - - = ∂ ∂ Ψ ∂ μ + ∂ Φ ∂ μ + λ + ∂ Φ ∂ λ = σ ω elsewhere 0 b ct x if e b 2 p x x 2 x ) 2 ( x 1 t i 3 1 2 2 3 2 2 1 2 33 (11) 0 x x x x 2 2 3 2 2 1 2 3 1 2 13 = ∂ Ψ ∂ + ∂ Ψ ∂ + ∂ ∂ Φ ∂ μ = σ (12)
• at x 3 = , there are no reflections, which means that all wave motions are in the positive x 3 -direction.

We introduce a frame of reference moving with the load (Eq. 13), in order to eliminate the time variable. We then define a Fourier transform on the space variable relative to the motion (Eq. 14) :

x = x 1 -ct y = x 2 z = x 3 ( 13
) +∞ ∞ - β = β dx e ) z , x ( h ) z , ( h x i (14)
The Fourier transform is applied for Eqs (8), ( 9), [START_REF] Keer | Moving and simultaneously fluctuating loads on an elastic half-plane[END_REF] and [START_REF] Gunaratne | Response of a layered elastic medium to a moving strip load[END_REF]. Taking into account that reflexions from z = are disallowed, the Helmholtz potentials and then the longitudinal and vertical displacements can be obtained in the transform domain. The actual displacements in the (x,y) domain are deduced by means of the inverse Fourier transform. To perform this inverse Fourier transform, we use the well-known Fast Fourier Transform (F.F.T.) algorithm [START_REF] Brigham | The Fast Fourier Transform[END_REF]. The displacements on the surface of the ground are then given by the following integral expressions :

• for the longitudinal displacement :

( ) t i x i 2 1 2 2 2 e d e 2 ) ( F b sin b 2 iP U ω β ∞ + ∞ - β α α - α + β β β πμ - = (15)
• for the vertical displacement :

t i x i 2 2 2 1 e d e 1 k k ) ( F b b sin 2 P W ω β ∞ + ∞ - β - β β α β β πμ - = (16)
with

2 2 i 2 2 i k 1 k β - - β = α , i = 1,2 (17)
i i c k , c k ω = ω = (18)
( )

2 2 1 2 2 2 2 4 ) , ( F β α α - β + α = γ β (19)
We note that if c = 0 m/s, Eqs ( 15) and ( 16) reduce to the known expressions for a stationary harmonic vibrating strip load given for example by Le Houédec [START_REF] Houédec | Réduction et propagation dans le sol des vibrations dues au trafic routier urbain : cas particulier des chaussées sur fondation élastique[END_REF], and equation ( 19) reduces to the well known Rayleigh function ; thus F is the Rayleigh function adapted to a moving load.

For the numerical results, the data concerning with the soil characteristics are given in Table 1. Particularly, the soil A is chosen from a British Rail site at which the constants have been measured. For this soil, the values of the Rayleigh, shear and compression wave speeds are c R = 242 m/s, c 2 = 263 m/s and 459 m/s respectively. The next results describe the behaviour of the vertical response for different load speeds. Concerning with the soil A, five different frequencies are chosen for the moving load : 8 Hz, 16 Hz, 32 Hz, 48 Hz and 64 Hz respectively. For the first results showing the effect of load speed on the vertical motion of half-space, this one is non-dimensionalized by strip length. In Fig. 2, the curves are perfectly symmetric and we note that the amplitude of the vertical displacement increases as the excitation frequency decreases. In Fig. 4, all curves are fully modified, particularly behind the load. In front of the load, the responses vanish practically. Besides, increasing the load speed further (for example c = 500 m/s), we can observe that the maximum displacement is almost equivalent for all frequencies. Fig. 4 shows also the appearance of oscillations behind the load and the apparent wavelength of these oscillations depends on both frequency and speed of the load. The following parametric study investigates the factors which describe this phenomenon. For this, we define the Mach number M R as the ratio between the load speed and the Rayleigh wave speed :

Soil A Soil B Soil C Soil D E(Pa) ρ(kg.m -3 ) ν η c R (m/s) c 2 (m/s) c 1 (m/s)
With the data given in Table 1 and for a frequency equal to 64 Hz, Fig. 5 shows that the maximum amplitude, non-dimensionalized against appropriate Rayleigh wavelength λ R , varies with different soil parameters. However, the apparent « wavelength » of oscillation is the same for all soil types, which will only be true for displacements over a half-space due to the dominant nature of the Rayleigh wave, Ewing et al [START_REF] Ewing | Elastic waves in layered media[END_REF]. In Fig. 6, for high Mach numbers, M R > 1, the non-dimensional vertical displacements are plotted against ( )

1 M x 2 R R - λ
; which shows that the « wavelength » of the vibrations of displacements behind the load is directly proportional to this factor similar to results obtained in acoustics, Morse and Ingard [START_REF] Morse | Theoretical Acoustics[END_REF]. The only material parameter which affects the apparent « oscillation » of the resultant wave is the Mach number M R . Thus, if we denote by λ* the « period » of non-dimensional oscillations and by λ ˆ the wavelength of the actual vertical displacements, we observe :

( ) 1 M 2 1 * 2 R R - = λ λ = λ
, and hence ( )

1 2 R 2 R 2 2 c c c ˆβ - β π = ω - π = λ (21) with R 1 c c + ω = β and R 2 c c + ω = β (22)
So, it is seen that the wavelength of the resultant wave is a combination of the two Rayleigh wavelengths. Finally, the wavelength of displacements, λ ˆ, depends on three factors : the Rayleigh wave speed c R , the load speed c and its frequency ω ; λ* depends on one factor, the Mach number M R . This justifies the interest of the factor ( )

1 M 1 2 R R - λ
in the supersonic region.

Strip load moving on an elastic layer overlaying a rigid foundation

Now, an other model is considered : the strip load rests on an homogeneous, isotropic, viscoelastic layer, with material properties E (Young's modulus), ρ (density), v (Poisson's ratio) and « modified » hysteretic damping as presented earlier ; this layer is overlaying a rigid foundation. The two-dimensional model is presented in Fig. 7. In this case, the transformed potentials are written as follows :

z z 1 1 e ) ( C e ) ( A α α - β + β = Φ (23) z z 2 2 e ) ( D e ) ( B α α - β + β = Ψ ( 24
)
where the functions -at the surface of ground, we obtain the same conditions as previously, -the displacements are equal to zero at the interface with the rigid foundation.

) ( A β , ) ( B β , ) ( C β
Combining the dynamic stiffness matrix analysis with the earlier approach (change of variable and Fourier transform), we obtain the following system :

[T]{ } { } τ = u ( 25
)
with

{ } { } t ) h ( u ), h ( w i ), 0 ( u ), 0 ( w i u = (26)
and

{ } { } t zx zz zx zz ) h ( ), h ( i ), 0 ( ), 0 ( i σ σ σ - σ - = τ (27)
and [T] is the symmetric 4 × 4 stiffness matrix whose terms are depending on λ and μ (Lamé constants), β, α i (Eq. 17), k and k i (Eq. 18) and h (layer depth). Because of computational difficulties due to the exponential terms in the matrix [T], it is necessary to divide the layer into several sub-layers. Morever, these sub-layers can have different material properties adapting for each sub-layer the matrix [T], and then assembling all the different matrices as in finite element analysis. The simplest case is of just two sub-layers with 2d = h : [START_REF] Jones | Ground vibration in the vicinity of a strip load : an elastic layer on a rigid foundation[END_REF] For c = 0 m/s, equation [START_REF] Jones | Ground vibration in the vicinity of a strip load : an elastic layer on a rigid foundation[END_REF] reduces to the expressions given by Jones and Petyt [START_REF] Jones | Ground vibration in the vicinity of a strip load : an elastic layer on a rigid foundation[END_REF] for a stationary vibrating strip load.

σ σ σ - σ - = × × × × × × × × × × × × × × × × × × × × × × × × × × × × ) d 2 ( ) d 2 ( i 0 0 ) 0 ( i ) 0 ( i 0 0 ) d ( u ) d ( w i ) 0 ( u ) 0 ( w i 0 0 0 0 0 0 0 0 zx zz zx zz
The influence of the layer's depth was studied by several authors, for example Waas [START_REF] Waas | Linear two-dimensional analysis of soil dynamics problems in semi-infinite layered media[END_REF] or Laghrouche [START_REF] Laghrouche | Simulation numérique de propagtion d'ondes dans les sols -Application à l'isolation vibratoire[END_REF]. For M R = 0 and for vertical vibrations, the critical depths are given by : .

... 3 , 2 , 1 n , f 4 c ) 1 n 2 ( h i n , i = - = (29)
where i = 1, 2 or R according to the wave type, and f is the load frequency. We can nondimensionalize the results by noting

R h * h λ =
, and in this case, * n , i h is depending on only (Poisson's ratio). For the soil A, the first adimensional critical depth is obtained for * 1 h = 0.474, very closed to 0.5. This result is verified in Fig. 8 where we notice that the vertical displacements are really higher for h* = 0.5, than for other values of h*. When the load speed increases, the responses of the soil surface are modified, as shown in Fig. 9. For the sub-Rayleigh case (M R = 0.5 < 1), the higher displacement is always obtained for h* = 0.5 behind the load, but now the curves are not any more symmetric and in front of the load conclusions are not distinct. For the super-Rayleigh case (M R = 2 > 1), in front of the load displacements are practically cancelled, and behind the load we notice the classical oscillations with on average a more important amplitude for h* = 0.5. For h* = 0.5, Fig. 10 shows a comparison of the vertical surface displacements for three values of M R : M R = 0 (stationary case), M R = 0.5 (sub-Rayleigh case) and M R = 2 (super-Rayleigh case). For the stationary case, we obtain more important adimensional displacements that for moving loads : 8.8 × 10 -9 for M = 0, 5.0 × 10 -9 for M R = 0.5 and 2.5 × 10 -9 for M R = 2. These values are obtained directly under the load, but beyond some wavelengths behind the load, we notice on the contrary opposite conclusions. In this case, equations are deduced taking into account the two previous approaches. First, we use a similar matrix system, as given in [START_REF] Jones | Ground vibration in the vicinity of a strip load : an elastic layer on a rigid foundation[END_REF], where displacements at the interface with the half-space are now unknown : [START_REF] Laghrouche | Simulation numérique de propagtion d'ondes dans les sols -Application à l'isolation vibratoire[END_REF] Besides, at the interface between the layer bottom and the half-space, we write the equality for stresses and for displacements. After the results obtained for a half-space, we can relate the displacements on the half-space surface according to the stresses, by means of a matrix equation : [START_REF] Jones | Ground vibration in the vicinity of a rectangular load on a half-space[END_REF] where [R] is a 2 × 2 stiffness matrix whose terms are depending on same data as [T].

σ σ σ - σ - = × × × × × × × × × × × × × × × × × × × × × × × × × × × × ) d 2 ( ) d 2 ( i 0 0 ) 0 ( ) 0 ( i ) d 2 ( u ) d 2 ( w i ) d ( u ) d ( w i ) 0 ( u ) 0 ( w i 0 0 0 0 0 0 0 0 zx zz zx zz
= σ σ ) 0 ( u ) 0 ( w i ] R [ ) 0 ( ) 0 ( i 1 1 1 zx 1 zz
Briefly, to simplify Eq. [START_REF] Laghrouche | Simulation numérique de propagtion d'ondes dans les sols -Application à l'isolation vibratoire[END_REF] in the objective to obtain six equations with six unknowns, we perform the two next operations :

first, we substract the matrix [R] from the last matrix group 2 × 2 of the global stiffness matrix [T], allowing to eliminate the two last terms of the stress vector ; -secondly, we introduce the boundary conditions for the moving load.

For the numerical results, we consider the case of a soft layer overlaying a more rigid half-space (Table 2 For these data, in Fig. 12, we obtain similar results as for Figs 8 and9, except in the case of M R = 0.5 where the maximum displacement is not obtained for h* = 0.5 but for h* = 0.25. So, for a moving load, the reference linked to the critical depth disappears, even for the subsonic regime. Fig. 12 -Adimensional displacements for an elastic layer overlaying a half-space.

Three-dimensional problem -Case of a rectangular load 2.2.1. Rectangular load moving on an elastic half-space

The model considered is shown in Fig. 13. Fig. 13 -Geometry of the model.

In this case, all data and assumptions written in the paragraph 2.1.1 are available. We only introduce :

-the dilatation

∂ + ∂ ∂ + ∂ ∂ = Θ (32)
-the double Fourier transform on the space variables x and y :

dy dx e ) z , y , x ( f 2 1 ) z , , ( f ) y x
( i γ + β - ∞ + ∞ - ∞ + ∞ - π = γ β (33)
with the above definition, the changes of variables and functions and introducing :

2 2 i 2 2 2 i c 1 k ω β - - γ + β = α , i = 1,2 (34)
with i i c k ω = , we finally obtain :

0 dz d 2 1 2 2 = θ α - (35)
Taking into account the condition to cancel the existence of waves propagating in the negative z-direction, the elastodynamic equation of half-space leads to the following system :

α - γ β μ μ + λ - = α - α - 1 z 2 2 2 2 i i Ae w v u dz d 1 (36)
The solution of Eq.( 36) can be obtained by a combination of the solution of the homogeneous form of (36) and a particular integral in this equation :

+ α γ - β - ω β - = α - α - D C B e i i c 1 k Ae w v u z 1 2 2 1 z 2 1 ( 37
)
If P is the total force acting on the rectangle and is uniformly distributed over it, the Fourier transform gives :

( ) ( ) 2 2 2 1 2 2 2 2 2 4 ) , ( F γ + β α α - γ + β + α = γ β ( 42
)
Transformed displacements are given by :

α γ - β - μ σ γ β γ + β + α = α - 1 z zz 2 2 2 2 i i e ) , ( F w v u 1 ( ) α γ + β γ β μ σ γ β α α + α - 2 2 2 z zz 1 i e ) , ( F 2 i 2 2 (43)
At the surface of the ground, the vertical displacement is written as follows :

( ) γ β γ β πμ ω β - α σ = γ + β ∞ + ∞ - ∞ + ∞ - d d e ) , ( F 2 / c 1 k ) 0 ( w ) y x ( i 2 2 2 1 zz (44)
that is to say in the case of a rectangular load :

( ) γ β γ β ω β - α γ γ β β μ π = γ + β ∞ + ∞ - ∞ + ∞ - d d e ) , ( F / c 1 k b b sin a a sin P 4 1 ) 0 ( w ) y x ( i 2 2 2 1 2 (45)
For a speed egal to zero, we can find results obtained by Jones and Petyt [START_REF] Jones | Ground vibration in the vicinity of a rectangular load on a half-space[END_REF], and for a stationary harmonic load we obtain the same results as Thomson and Kobori [START_REF] Thomson | Dynamical compliance of rectangular foundations on an elastic half-space[END_REF].

For the numerical results, we use the characteristics data given for the soil A (Table 1). In the transform domain, for non-zero damping, the real part of the transformed displacements over the surface of the ground is plotted in Fig. 14, showing a sub-Rayleigh case and a super-Rayleigh case. According to the hypothesis, we obtain the visualisation of an ellipse or an hyperbola. The comparison between the theoretical equation of the conic and the numerical results shows a good agreement. In this figure, only the contribution of the Rayleigh wave can be seen because it is dominant at the surface. When M R increases, the maximum of non dimensionalized displacements amplitude increases a little, from 1.63 × 10 -9 for M R = 0 to 1.74 × 10 -9 for M R = 1, but then it decreases continuously to 1.36 × 10 -9 for M R = 2. Morever, for low values of M R , the maximum is situated at the middle of the rectangle, whereas for greater Mach numbers, the maximum moves to the back edge of the rectangle.

When M R is low, the distribution of displacements is focused around the load in the near field, whereas with increasing M R , the distribution spreads behind the load and on the sides (Fig. 16). For M R > 1, the displacements are contained within the two Mach lines, and are almost equal to zero in front of the load due to a Rayleigh wave speed lower than the load speed.

Rectangular load moving on an elastic layer over a rigid foundation

The ground is modelled as an elastic layer of depth h overlaying a rigid foundation. Now, the boundary conditions are changed :

• for z = 0 = σ = σ < < ω - = σ 0 yz xz a y , b x t i e ab 4 P zz • for z = h : u 1 = u 2 = u 3 = 0
We use the same method to solve the problem i.e. changes of variables and functions, and double Fourier transform on the space variables x and y. Now, waves can propagate in the negative z-direction due to reflections at the boundary with the bedrock, and the resolution of [START_REF] Sheng | Ground vibration generated by a load moving along a railway track[END_REF] gives in this case :

+ α γ - β - ω β - = α - α - G D C e i i c 1 k Ae w v u z 1 2 2 1 z 2 1 + α - γ - β - ω β - + α α H F E e i i c 1 k Be z 1 2 2 1 z 2 1 (46)
where A, B, C, D, E, F, G, H are depending on β and γ. The constants G and H are calculated with the aid of the dilatation expression which gives :

γ + β α - = γ + β α = ) F E ( i H ) D C ( i G 2 2 (47)
and then the transformed displacements can be written with matrix expression as follows : In the first case, as for the 2D study (Figs 8 and9), for a critical depth closed to 0.5, the displacements are clearly higher, but mainly behing the load. Besides, except for h* = 0.1, under the load (x = 0), all maximum displacements are practically equal. Also, all curves are almost symmetric.

In the second case, we can see as usually the existence of disturbance and higher displacements behing the load. Then the condition h* = 0.5 don't give the higher responses and consequently for the super-Rayleigh case it is difficult to take into account the critical depth.

Multilayered soils

Now, the ground has several horizontal layers overlaying an elastic half-space. To obtain a solution of such a problem, we just have to combine the results of the two previous paragraphs (as for the 2D approach). We construct a global stiffness matrix using the same assembly process as in finite element analysis, (continuity of displacements and strains at each interface). Finally, using the continuity of displacements and strains at the interface between the last layer and the half-space, we obtain the following matrix system : yielding the transformed displacements at each interface. For the numerical analysis, for example two layers over an elastic foundation, the results can be deduced after a very long calculation, but are above all too complicated to explain the influence of characteristics concerning with the layers (material properties, depth, ...) and the load (frequency, speed particularly).

APPLICATION TO RAILWAY TRAFFIC

Previous studies are concerned with the dynamic behaviour of soil media for a moving harmonic loading. For the ground, a lot of models can be used : two-or threedimensional, half -space or multilayered soils. For a real excitation relative to the railway traffic, it is necessary to choose a more suitable model than a simple harmonic load.

Description of the track

Several track models have been already described in the literature. Some track models are based on Green's functions (Krylov [33]), others on a description in finite elements (Bode and Savidis [START_REF] Bode | Three dimensional time analysis of moving loads on railway tracks on layered soils[END_REF]). Others models using the formalism of Fourier transform consider the discrete nature of the track due to sleepers spacing (Metrikine and Popp [START_REF] Metrikine | Vibration of a periodically supported beam on an elastic half-space[END_REF]) or assume the homogeneity and the continuity of the track lying on the ground (Sheng et al. [START_REF] Sheng | Ground vibration generated by a load moving along a railway track[END_REF]).

Our numerical model is built with a semi-analytical method for rail, pads and sleepers. In a range of low frequencies (below 80 Hz), the track can be seen as a simple massspring system. So the continuous track model used is shown in Figure 19. For the rail, the Euler beam theory is used. In front of the foreseen wavelengths, the sleepers are described as a continuous mass system. Pads and ballast have a visco-elastic damping with an arbitrary frequency dependence. Consequently, in this mass-spring system, waves generated in the rail transmit energy to the ground through the sleepers. Note that the reaction of the ground to sleepers is a dependent function of the frequency and interact with the track.

x=x 1 -ct c F S O ) t , x ( ENG σ ) t , x ( WAG σ F B (x,t) Rail (EI,m R ) Pads (k P ,η P ) Sleepers (m S ) Ballast (E B ,ν B ,ρ B ,η B ) Layered ground (E n ,ν n ,ρ n ,η n ,h n ) x 3 =z W R W S Fig. 19 -Track-Soil interaction model
In a two-dimensional or three-dimensional model the ballast can be described as the first layer of the soil. Also railpads and ballast can prove non linear properties in some ranges of frequencies. Several models ( [START_REF] Suiker | Dynamic behaviour of a layer of a discrete particles -Part 2 : Response to an uniformly moving, harmonically vibrating load[END_REF], [START_REF] Nguyen | Comportement dynamique de structures non-linéaires soumises à des charges mobiles[END_REF]) tend to approach these effects. To describe the rail [START_REF] Jones | Simulation of ground vibration from a moving harmonic load on a railway track[END_REF], a mathematical model based on the Euler beam gives the following equation :

+ ∂ ∂ + ∂ ∂ 2 1 R 2 R 4 1 1 R 4 t ) t , x ( w m x ) t , x ( w EI [ ] = - ) t , x ( w ) t , x ( w k 1 S 1 R P < - σ lse e 0 b ct x if ) t , x ( 1 AIN TR (54)
Besides, equation for sleepers is written :

+ ∂ ∂ 2 1 S 2 S t ) t , x ( w m [ ] ) t , x ( F ) t , x ( w ) t , x ( w k 1 S 1 R 1 S P - = - (55)
where w R (x 1 ,t), w S (x 1 ,t) are respectively the vertical displacements of the rail beam and the sleeper mass. F S (x 1 ,t) is the force per unit length of track between the sleeper and the ballast.

Two -dimensional problem (multilayered soil)

In this case, vertical and longitudinal transformed stresses are given by :

ψ β - ∂ ψ ∂ - ∂ φ ∂ β μ = σ 2 2 2 xz z z i 2 and z i 2 z ) 2 ( 2 2 2 zz ∂ ψ ∂ μβ + φ λβ - ∂ φ ∂ μ + λ = σ (56)
where φ and ψ are Helmoltz's potentials.

The resolution method is described in figure 20 where the fitted phase angle is defined relating to Helmoltz's potentials such as : The application of this phase translation allows a good conditioning of the stiffness matrix of ground, then subsection of layers into sub -layers isn't necessary. A change of variables (translation in the mobile reference linked to the load) and functions (evaluation of the steady state solution) as well as the spatial Fourier transform allows the calculation of vertical displacement in the wave number domain.

+ = ψ + = φ α - - α α - - α z ) h z ( z ) h z (
Thus, for the ground, we can write a relationship between the displacement and stress for each layer interface (including the ballast) : Thus, vertical displacements on the surface of the ballast and the ground can be written, respectively :

+ + + + + + + + ) h ( * u ) h ( * w i ) h ( * u ) h ( * w i ) 0 ( * u ) 0 ( * w i R T R T T T 0 0 R T R T T T 0 0 T T T T T T
) ( F ] T [ ) ( w S 1 11 0 z β = β * - * = and ) ( F ] T [ ) ( w S 1 31 h z β = β * - * = (59)
Assuming that the model for sleepers is continuous in front of the wavelength then we can obtain : [ ]

β = β β - = β Ω - + β - β σ = β - β + Ω - β * - * * ) ( F ] T [ ) ( w ) ( F ) ( w ) m k ( ) ( w k ) ( )] ( w k ) ( w ) k m EI (
) k m EI ( ] T [ k ) m k )( k m EI ( ) ( k ] T [ ) ( w P 2 0 R 4 1 11 2 P 2 0 S P P 2 0 R 4 AIN TR P 1 31 h z + Ω - β + - Ω - + Ω - β β σ = β - - * = (61) ) ( w h z β * =
represents the vertical displacements in the wave number domain below the ballast (ground surface).

An analytic solution exists when the soil is considered as infinite in the positive zdirection. In this case, the Rayleigh's function leading to an infinite vertical displacement is determined by the cancellation of the denominator. It corresponds exactly to the speed of Rayleigh waves in the medium. In the case of a layered soil, displacements in the real numbers domain are calculated numerically by means of an inverse Fourier transform.

The vertical force exerted on the rail by a moving train is a complicated function of time. The spectrum of this function is wide but a significant part of the energy is concentrated between 0 and 20 Hz. Besides, two components characterize the source of vibrations on the ground surface : a static component coming from the weight of axles of the train acting directly on the track and a dynamic component due to flaws of the track and wheels. For a simple wavelength of the rail, the train vibrates at the frequency 2πc/λ 0 and each wheel exerts a vertical force t i 0 0 e P Ω . The rail is then constrained to its surface by a harmonic and moving function representing the train stress. This function depends on the weight of each axle load, the excitation frequency Ω 0 and the length parameters coming from bogies and axles. The railway train (engine and carriages) and its stress function are plotted in Figure 21. In the wave number, main peaks translate a frequency passage of a bogie or wheel. Finally this formulation is introduced in the model (at rail surface). For the numerical results, data for the ground are given in Table 3 In these conditions, the critical speed is close to 100 metres per second. Rails Young's modulus is equal to 2,11.10 11 N/m², the moment of inertia to 3055 cm 4 and the mass per unit of track length to 120 kg/m. The stiffness of the pads between the rail and sleepers is equal to 3,5.10 8 N/m with a damping coefficient equal to 0,15. The mass of wooden sleepers per unit of track length is equal to 490 kg/m. The ballast is a stiff layer compared to the soft peat with a Young's modulus equal to 400.10 6 N.m -2 , a Poisson ratio equal to 0.3 and a density equal to 1800 kg/m 3 . Damping is assumed to be hysteretic with a factor equal to 0.1.

Figure 22 shows a comparison between the displacements in the real domain for each speed regime : static, sub-Rayleigh and super-Rayleigh passing through the critical speed. The excitation frequency is chosen equal to 10 Hz according to real train spectra. The Sub-Rayleigh speed is equal to 0.5×c R and the super-Rayleigh speed to 1.5×c R where c R is the speed of the Rayleigh wave. All results are divided by the Rayleigh wavelength λ R for non-dimensional study. For non moving load, displacements are distributed symmetrically in the front and behind the load. Increase of speed involves increasing of vertical displacements amplitude and then the increase of the wavelength and oscillation behind the load. Increase of vibration can be noted for super-Rayleigh regime.

Besides, in super-Rayleigh regime, when the train speed on track exceeds the Rayleigh wave speed at the soil surface, Mach cones appear on the soil for each bogie load. So, Mach cones due to P and S-waves can be found below the surface and Mach cone due to Rwaves at the soil surface. Also, the mutual interaction between the engine, the carriage and the track was studied Hung and Yang [START_REF] Hung | Elastic waves in visco-elastic half space generated by various vehicle loads[END_REF], Wu and Shih [START_REF] Wu | Dynamic responses of railway and carriages under the high speed moving loads[END_REF] for example.

Three-dimensional problem

Then, the ground is modelled by either a half-space or a multilayered soil including or not the ballast. For the resolution and calculation of surface displacements, we combine the previous approaches taking into account the train, track and ground models, and the classical methods (change of variables, Fourier transform, choice of « fitted phase argument », numerical inverse Fourier transform). Here, we limit the presentation to some numerical results. For example, Fig. 23 gives the vertical surface displacements obtained for an unit load (frequency equal to 40 Hz) moving on a track connected with a ground constituted by an elastic sandy layer and a half-space. Similarly, Fig. 24 shows the vertical surface displacements deduced for a Corail train. In this case, the frequency is equal to 10 Hz, the ballast is modelled as the first layer of soil, and the ground is constituted by a peaty material. We confirm the previous results obtained at the ground surface for various speeds. Particularly, increase of speed implies more important vertical displacements and oscillations behind the train. Besides, we can see the Mach cones for each bogie group.

2. 1 .

 1 Two-dimensional problem -Case of a strip load 2.1.1. Strip load moving on an elastic half-space

Figs 2 ,

 2 Figs 2, 3 and give the results for c = 0 m/s (stationary load), c = 200 m/s (subsonic case) and c = 400 m/s (transonic and super-Rayleigh case).

Fig

  Fig. -Non-dimensional vertical motion for c = 0 m/s : solid line f = 8 Hz, dash -dot f = 16 Hz, dashed f = 32 Hz, light-dot f = 48 Hz, heavy -dot f = 64 Hz.

FigFig. 4 -

 4 Fig. -Non-dimensional vertical motion for c = 200 m/s : solid line f = 8 Hz,

Fig. 5 -

 5 Fig. 5 -Non-dimensional vertical motion for four soil material parameters at load speed c = 2 c R : solid line soil A, dash-dot soil B, dashed soil C, light-dot soil D.

Fig. 6 -

 6 Fig. 6 -Non-dimensional vertical motion at frequency f = 64 Hz for four load speeds :

Fig. 7 -

 7 Fig. 7 -Geometry of the two-dimensional model.

  and ) ( D β are determined taking into account the new boundary conditions :

Fig. 8 -

 8 Fig. 8 -Critical depth of the elastic layer (M R = 0).

Fig. 9 -

 9 Fig. 9 -Critical depth of the elastic layer (M R ≠ 0).

Fig. 10 -

 10 Fig. 10 -Comparison of displacements for various values of M R .

=

  Then, we can deduce expressions for A, B, C, D : for a moving load in the 3D case :

Fig. 14 -Fig. 15 -

 1415 Fig. 14 -Contour levels of the real part of the transformed vertical displacements for M R = 0.5 and M R = 2 .

Fig. 16 -

 16 Fig. 16 -Contour levels of the real part of vertical displacements for M R = 0.5 and M R = 2.

Fig. 17 -

 17 Fig. 17 -Vertical displacements for M R = 0.5.

Fig. 18 -

 18 Fig. 18 -Vertical displacements for M R = 2.

Fig. 20 -

 20 Fig. 20 -Resolution method

  BAL ] is the stiffness matrix of the ballast, [T n ] is the stiffness matrix taking into account the n layers of ground and [R] is the stiffness matrix of the halfspace. The stiffness matrix [T] is given by the product of the inverse of the displacement matrix [Q] by the stress matrix [S].

  of which is deduced the displacement on the soil surface (below the ballast) :

Train

  Fig 21-Railway excitation.

Fig. 22 -

 22 Fig. 22 -Maximum vertical displacements at the soil surface near the track for current train (c = 0 : dashed line, c=0.5×c R : dash-dotted line, c=1.5×c R : solid line).

Fig. 23 -

 23 Fig. 23 -Vertical surface displacements for a layered soil (static, sub-Rayleigh and super-Rayleigh regimes respectively).

Fig. 24 -

 24 Fig. 24 -Vertical surface displacements for a Corail train (on the left : sub-Rayleigh regime, on the right : super-Rayleigh regime).

CONCLUSIONS • 9

 9 th International Conference on Civil and Structural Engineering Computing, Egmond ann Zee (Netherlands), 2-4 September 2003. * Special lecture (invited paper) : D. LE HOUÉDEC, B. PICOUX « Soil-structure interaction for moving loads : application to railway traffic » * Communication orale : B. PICOUX, D. LE HOUÉDEC « Numerical and experimental comparison of D-model for the study of railway vibrations » * Organisation de la session spéciale « Dynamics of structures for moving loads » (une douzaine de présentations).

Table 1 -

 1 Material properties

	2.69 × 10 8	5.38 × 10 8	10.76 × 10 8	2.0 × 10 8
	1550	3100	2000	1250
	0.257	0.257	0.257	0.257
	0.1	0.1	0.1	0.1
	242	242	426	232
	263	263	463	252
	459	459	809	441

Table 2 -

 2 ). Characteristics of the layer and the half-space.

	Soft layer	« Rigid » half-space

  .

		h (m)	E (N/m²)	ν	ρ (kg/m 3 )	η	c P (m/s)	c S (m/s)
	Ballast	0,8	300×10 6	0.3	1800	0.1	547	292
	Layer	7	269×10 6	0.257	1550	0.1	459	263
	Substratum	∞	2040×10 6	0.179	2450	0.1	950	594

Table 3 -

 3 Parameters for numerical investigation

with

Fourier transforming the stress-strain relations, the stresses in the wave number domain can be written in terms of the constants A, B, C, D, E and F :

Combining ( 48) and (50), we obtain the following system :

With the different choices made for the expressions of the vectors of displacements and strains,

For the numerical results, Figs 17 and 18 show the amplitude of vertical displacements along the line y = 0 for different layer's depths. Fig. 17


Chady Ghnatios
email: chady.ghnatios@ec-nantes.fr

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Elías Cueto
email: ecueto@unizar.es

Adrien Leygue
email: adrien.leygue@ec-nantes.fr

Arnaud Poitou
email: arnaud.poitou@ec-nantes.fr

Optimizing Composites Forming Processes by Applying the Proper Generalized Decomposition

Keywords: PGD, parametric model, model reduction, optimisation

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Despite the spectacular progress in the context of computational mechanics in the last few decades, the efficient treatment of multidimensional models needs further developments [START_REF] Chinesta | On the use of model reduction in material and processes simulation: A way for realizing simulation dreams[END_REF]. In fact, this type of problems can't be solved by a classical technique like the finite element method, because of the time cost involved. Moreover, nowadays the fine description of the thermodynamics and structure of materials introduces some specific challenges related to the impressive number of degrees of freedom required for solving problems, which generates such meshes of multidimensional models. Furthermore, usual models encountered in computational physics are solved for specific conditions, but can also be transformed to multidimensional models allowing the computation of a very general solution. For example, in the pultrusion process, we can solve the equation of heat for all possible boundary conditions. Thus, we can optimize the system to a global optimum. Solving this equation is now possible due to the new numerical method so called the "Proper General Decomposition" or PGD, as we will describe later. In fact, when using PGD, the complexity of the problem scales linearly with its dimension, enabling the circumvention of the curse of dimensionality [START_REF] Chinesta | Archives of Computational Methods in Engineering[END_REF]. This technique has been successfully applied in a variety of linear, non linear, stationary and transient problems [3] [4].

In this paper, we will show the interest of multidimensional models by solving problems with parametric boundary conditions. Obviously, usual finite element procedures, for each choice of the boundary condition prescribed on a region of the boundary, a thermal problem must be solved. Thus, optimization becomes too expensive from the computation time point of view. The PGD allows solving only once a general model for any boundary condition prescribed on any chosen region of the mesh. Optimization as well as inverse identification becomes an easy matter.

PARAMETRIC BOUNDARY CONDITIONS

Often optimizing a process can be modelled as the optimal choice of a set of parameters, including sometimes some boundary conditions. This is, for example, the case of pultrusion, one of several manufacturing processes for polymer composites [5]. There is also too many other interests in computing parametric boundary conditions, which are currently being tested.

Illustrating the construction of a solution using PGD

In what follows, we show the construction of the solution, using PGD, of the following simple parametric heat transfer equation as an example:

∂ u ∂t -k∆u -f = 0 (1)
where (x,t, k) ∈ Ω × I × ℑ, and for the sake of simplicity the source term is assumed to be constant, i.e. f = cte. Because the conductivity is considered unknown, it is assumed as a new coordinate defined in the interval ℑ. So, in this case, the solution of Eq. ( 1) is searched under the form:

u (x,t, k) ≈ i=N ∑ i=1 X i (x) • T i (t) • K i (k) (2)
To start the algorithm, we assume that the approximation at iteration n is already known. At present iteration we look for enrichment, thus u writes:

u n (x,t, k) = i=n ∑ i=1 X i (x) • T i (t) • K i (k) + R (x) • S (t) •W (k) (3)
where

R (x) • S (t) • W (k) are respectively X n+1 (x) • T n+1 (t) • K n+1 (k). Before solving
the resulting non linear model related to the calculation of these three functions, a model linearization is performed. The simplest choice consists in using an alternating directions fixed point algorithm. First of all, we proceed by assuming S (t) and W (k) given at the previous iteration of the non-linear solver and then computing R (x). From the just updated R (x) and previously used W (k) we can update S (t). Finally from the just computed R (x) and S (t) we update W (k). The procedure continues until reaching convergence. The converged functions R (x), S (t) and W (k) allow defining all needed functions :

X n+1 (x) = R (x), T n+1 (t) = S (t) and K n+1 (k) = W (k) [1].
Remark 1. The construction of each term in the sum (2) needs a certain number of iterations because of the non linearity of the problem related with the approximation (3). Let m = ∑ i=N i=1 m i be the total number of iterations involved in the construction of the separated approximation (2), which rarely exceeds ten. On the other hand the number of sums N needed to approximate the solution of a given problem depends on the solution regularity itself, but all the experiments carried out reveal that this number ranges from few tens to few hundreds. Thus, we can conclude that the complexity of the solution procedure is of some hundreds of 3D solutions (the cost related to the one dimensional problems being negligible with respect to the 3D problems [START_REF] Chinesta | Archives of Computational Methods in Engineering[END_REF]). Now, if we assume a classical approach, we should solve a 3D problem at each time step and for each value of the parameter k. In [6] we proved that the CPU time savings by applying the PGD can be of several orders of magnitude. Remark 2. This section shortly described is the most basic algorithm for constructing the sum [START_REF] Chinesta | On the use of model reduction in material and processes simulation: A way for realizing simulation dreams[END_REF]. There are other ways to construct solutions using PGD, with faster convergence [START_REF] Pruliere | Mathematics And Computer Simualtion[END_REF].

Basic example using parametric boundary conditions

To fixe the ideas, let's start by considering the basic one dimension differential equation that follows:

- ∂ 2 T ∂ x 2 = 1 (4)
Let's consider that we need to solve this equation for x ∈ Ω where Ω = [0; L] and for the following boundary conditions:

T (0) = 0 T (L) = y ( 5
)
where y is a variable parameter inside a domain, noted α. Thus, in the classical methods, we will have to solve this equation for every single value of y. In the PGD method, we can include y as an extra coordinate in the system, since this doesn't involve dramatically the complexity of the problem. Thus, we first write a function, f(x,y) that satisfies all boundary conditions. For our example, we can use:

f (x, y) = x.y L (6)
Since, as noted on the page before, T can be written as a sum of N products of functions X(x) and Y(y), and introducing the values of our boundary condition f(x,y) as first mode of the solution, T(x,y) writes now:

T (x, y) = X 1 (x) ×Y 1 (y) + i=N ∑ i=2 X i (x).Y i (y) (7)
where :

X 1 (x) ×Y 1 (y) = x.y L (8)
Note that X 1 (x) and Y 1 (y) can be found using a singular value decomposition. Now we will go back to write the weak form of the problem, the trial function and test function write respectively as:

T (x, y) = X 1 (x)Y 1 (y) + i=n ∑ i=2 X i (x)Y i (y) + R(x)S(y) (9) T * (x, y) = R * (x)S(y) + R(x)S * (y) (10)
• By replacing equations ( 9) and (10) into equation ( 4) and considering S(y) known, so S * = 0 we have the first weak form to solve:

α Ω RY 1 ∂ 2 X 1 ∂ x 2 + R i=n ∑ i=2 ∂ 2 X i ∂ x 2 Y i + ∂ 2 R ∂ x 2 RS -R = 0 (11)
• Now repeating the same procedure considering R(x) known, thus R * = 0 and the second weak form:

Ω α SY 1 ∂ 2 X 1 ∂ x 2 + S i=n ∑ i=2 ∂ 2 X i ∂ x 2 Y i + ∂ 2 R ∂ x 2 S 2 -S = 0 (12)
Now, we should iterate equations ( 11) and (12) until convergence. Solutions are found for homogenous Dirichlet boundary conditions. This is a must, because our boundary conditions are already satisfied by first mode. The new modes will converge to the good solution taking in consideration the boundary conditions. Results of this problem are shown in figure 1. 4) for parametric boundary condition cases [START_REF] Martin | [END_REF]. In fact, each have chosen his own boundary conditions to solve a two dimensional heat transfer equation based on Fourier's law. This axisymmetric equation can be written as:

ρc ∂ T ∂t + u ∂ T ∂ z = 1 r ∂ ∂ r k r ∂ T ∂ r + ∂ ∂ z k z ∂ T ∂ z + q (13)
where k r and k z are the thermal conductivities in radial and axial directions respectively, q is the internal heat generated by the resin curing reaction, ρ is the density, c is the specific heat and u the pull speed. We will solve this problem in a die with three heating stages, each one is considered as a Dirichlet boundary condition. Obviously, each one will be considered as a parametric boundary condition, and a new coordinate. The pull speed u is also a parameter of the problem, and a coordinate is assigned to it. Thus, the solution T writes as :

T (x, y, θ 1 , θ 2 , θ 3 , u) = i=N ∑ i=1 F i (x, y) × Θ 1i (θ 1 ) × Θ 2i (θ 2 ) × Θ 3i (θ 3 ) ×U i (u) (14)
while the boundary conditions write as :

         T (0, y, θ 1 , θ 2 , θ 3 , u) = T 0 T (x, y, θ 1 , θ 2 , θ 3 , u) = θ 1 f or x ∈ [ε; L 1 -ε] and (y = 0 or y = h) T (x, y, θ 1 , θ 2 , θ 3 , u) = θ 2 f or x ∈ [L 1 + ε; L 2 -ε] and (y = 0 or y = h) T (x, y, θ 1 , θ 2 , θ 3 , u) = θ 3 f or x ∈ [L 2 + ε; L 3 -ε] and (y = 0 or y = h) ∇T (L 3 , y, θ 1 , θ 2 , θ 3 , u) = 0 (15)
where L 1 and L 2 are the values of x coordinate at the entrance of heater II and III respectively. L 3 is the value of the x coordinate at the exit of the die. h is the width of the die. To solve this problem, we will use a set of four initial modes to define (15) and use the same technique as described on the page before. Solution for one particular case is shown in figure 2.

We compared the solution represented in figure 2 to the finite element one for same boundary condition, and we have found a difference less than 10 -5 for N going up to 30 modes.

CONCLUSION

In this paper, we summarized a way to construct the PGD, and we proposed a new approach for solving problems involving parametric boundary conditions. This new view of problems promises very important changes in optimization. As we have seen, to many mechanical optimization problems consist of choosing a set of parameters. However, computing solution for all boundary conditions is not enough to optimize the process. Further treatment needs to be performed in a huge multidimensional space, where human imagination is limited. Some algorithms are currently being tested, trying to reach a global optimum of each process.

Figure 1 .

 1 Figure 1. Solution of equation (4) for parametric boundary condition

Figure 2 .

 2 Figure 2. Solution for U 0 =100K, heater temperatures are respectively 400K, 300K and 200K for heater I, II and III. Pull speed is 0.1 m/s

APPLYING PGD FOR OPTIMIZATION

In this section, we will treat the problem of axisymmetric pultrusion. This is a continuous process to produce constant cross-sectional profile composites. For decades, engineers relied on experience to define optimal parameters for pultrusion process [START_REF] Martin | [END_REF]. However, with the fast growing use of this technique [9], more effective methods must be used for the design of this process. Many numerical models have been presented, each with his own advantages, but most of the simulation procedures were designed for particular


Constantin-Irinel Morȃrescu

Bernard Brogliato

Passivity-based tracking control of multiconstraint complementarity Lagrangian systems

In this study one considers the tracking control problem of a class of nonsmooth fully actuated Lagrangian systems subject to frictionless unilateral constraints. A passivitybased switching controller that guarantees some stability properties of the closed-loop system is designed. A particular attention is paid to transition (impacting) and detachment phases of motion. This work extends previous works on the topic as it considers multiconstraint n-degree-of-freedom systems.

I. INTRODUCTION

This paper focuses on the problem of tracking control of complementarity Lagrangian systems [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] subject to frictionless unilateral constraints whose dynamics may be expressed as:

   M (X) Ẍ + C(X, Ẋ) Ẋ + G(X) = U + ∇F (X)λ X 0 ≤ λ X ⊥ F (X) ≥ 0, Collision rule
(1) where X ∈ R n is the vector of generalized coordinates, M (X) = M T (X) ∈ R n×n is the positive definite inertia matrix, F (X) ∈ R m represents the distance to the constraints, C(X, Ẋ) is the matrix containing Coriolis and centripetal forces, G(X) contains conservative forces, λ X ∈ R m is the vector of the Lagrangian multipliers associated to the constraints and U ∈ R n is the vector of generalized torque inputs. For the sake of completeness we precise that ∇ denotes the Euclidean gradient ∇F (X) = (∇F 1 (X), . . . , ∇F m (X)) ∈ R n×m where ∇F i (X) ∈ R n represents the vector of partial derivatives of F i (•) with respect to the components of X. We assume that the functions F i (•) are continuously differentiable and that ∇F i (X) = 0 for X with F i (X) = 0. It is worth to precise here that for a given function f (•) its derivative with respect to the time t will be denoted by ḟ (•). For any function f (•) the limit to the right at the instant t will be denoted by f (t + ) and the limit to the left will be denoted by f (t -). A simple jump of the function f (•) at the moment t = t ℓ is denoted σ f (t ℓ ) = f (t + ℓ )f (t - ℓ ). The admissible domain associated to the system (1) is the closed set Φ where the system can evolve and it is described as follows:

Φ = {X | F (X) ≥ 0} = 1≤i≤m Φ i , C.-I.
Morȃrescu is with INRIA, BIPOP team-project, Inovallée, 655 avenue de l'Europe, 38330, Montbonnot, France Constantin.Morarescu@inrialpes.fr B. Brogliato is with INRIA, BIPOP team-project, Inovallée, 655 avenue de l'Europe, 38330, Montbonnot, France Bernard.Brogliato@inrialpes.fr where Φ i = {X | F i (X) ≥ 0} considering that a vector is non-negative if and only if all its components are nonnegative. In order to have a well-posed problem with a physical meaning we consider that Φ contains at least a closed ball of positive radius.

Definition 1: A singularity of the boundary ∂Φ of Φ is the intersection of two or more codimension-one surfaces

Σ i = {X | F i (X) = 0}.
The presence of ∂Φ can induce some impacts that must be included in the dynamics of the system. It is obvious that m > 1 allows both simple impacts (when one constraint is involved) and multiple impacts (when singularities or surfaces of codimension larger than 1 are involved). The collision (or restitution) rule in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], is a relation between the post-impact velocity and the pre-impact velocity. Among the various models of collision rules, Moreau's rule is an extension of Newton's law which is energetically consistent [START_REF] Glocker | Concepts for Modeling Impacts without Friction[END_REF] and is numerically tractable [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. For these reasons throughout this paper the collision rule will be defined by Moreau's relation [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]:

Ẋ(t + ℓ ) = (1 + e n ) arg min z∈TΦ(X(t ℓ )) 1 2 [z -Ẋ(t - ℓ )] T × M (X(t ℓ ))[z -Ẋ(t - ℓ )] -e n Ẋ(t - ℓ ) (2)
where Ẋ(t + ℓ ) is the post-impact velocity, Ẋ(t - ℓ ) is the preimpact velocity, e n ∈ [0, 1] is the restitution coefficient and T Φ (X) is the tangent cone to Φ at X [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF]. Denoting by T the kinetic energy of the system, we can compute the kinetic energy loss at the impact t ℓ as [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF]:

T L (t ℓ ) = - 1 -e n 2(1 + e n ) σ ⊤ Ẋ (t ℓ )M (X(t ℓ ))σ Ẋ (t ℓ ) ≤ 0 (3)
The structure of the paper is as follows: in Section 2 one presents some basic concepts and prerequisites necessary for the further developments. Section 3 is devoted to the controller design. In Section 4 one defines the desired (or "exogenous") trajectories entering the dynamics. The desired contact-force that must occur on the phases where the motion is constrained, is explicitly defined in Section 5. Section 6 focuses on the strategy for take-off at the end of constraint phases. The main results related to the closed-loop stability analysis are presented in Section 7. One example and concluding remarks end the paper.

The tracking control problem under consideration was studied in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] mainly in the 1-dof (degree-of-freedom) case and in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] in the n-dof case. Both of these papers consider systems with only one unilateral frictionless constraint. Here we not only consider the multiconstraint case but the results in Section 7 relax some very hard to verify condition imposed in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF].

The following standard notations will be adopted: || • || is the Euclidean norm, b p ∈ R p and b n-p ∈ R n-p are the vectors formed with the first p and the last n-p components of b ∈ R n , respectively. N Φ (X p = 0) is the normal cone N Φ (X) to Φ at X (see [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF]) when X satisfies X p = 0, λ min (•) and λ max (•) represent the smallest and the largest eigenvalues, respectively. We also note that an LCP is a system as 0 ≤ λ ⊥ Aλ + b ≥ 0, which has a unique solution for all b ∈ R n if and only if A ∈ R n×n is a P-matrix.

II. BASIC CONCEPTS

A. Typical task

In the case m = 1 (only one unilateral constraint) the system dynamics alternates free-motion phases Ω 2k when the constraint is not active(F (X) > 0), and constraint-motion phases Ω 2k+1 when the constraint is active (F (X) = 0). Between free and constraint phases the dynamical system always passes through a transition phase I k containing some impacts (more details can be found in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF]). In the case m ≥ 2 (multiple constraints) things complicate since the number of typical phases increases due to the singularities of ∂Φ that must be taken into account. Explicitly, the constraint-motion phases need to be decomposed in sub-phases where some specific constraints are active. As we shall see later the tracking during this phases does not present particular difficulties for different number of active constraints. Thus the goal is to control the system during a generic phase (constraint or not) and the passage between phases when the number of active constraints increases or decreases. We stress out that an impacting transition occurs when the number of active constraints increases but there is no impact (and no change of the dynamics) when the number of active constraints decreases. Therefore, without any loss of generality we study the following typical task:

R + = k≥0 Ω J k 2k ∪ I J k k ∪ Ω J ′ k 2k+1 , J k , J k+1 ⊂ J ′ k ( 4
)
where the superscript J k represents the set of active constraints during the corresponding motion phase. Throughout the paper, the sequence

Ω J k 2k ∪ I J k k ∪ Ω J ′ k
2k+1 will be referred to as the cycle k of the system's evolution. Furthermore for robustness reasons, during impacting transition phases I J k k we impose a closed-loop dynamics (with impacts) that mimics somehow the bouncing-ball dynamics (see e.g. [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]).

B. Stability analysis criteria

The system (1) is a complex nonsmooth and nonlinear dynamical system which involves continuous and discrete time phases. A stability framework for this type of systems has been proposed in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] and extended in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF]. This is an extension of the Lyapunov second method adapted to closedloop mechanical systems with unilateral constraints. In the sequel we introduce some basic concepts in order to clarify the framework. The trajectories playing a role in the dynamics and the design of the controller are: • X nc (•) -the desired trajectory of the unconstrained system (i.e. the trajectory that the system should track if there were no constraints). We suppose that F (X nc (t)) < 0 for some t, otherwise the problem reduces to the tracking control of a system with no constraints. • X * d (•) -the signal entering the control input and playing the role of the desired trajectory during some parts of the motion.

• X d (•) -the signal entering the Lyapunov function. This function is set on the boundary ∂Φ after the first impact of each cycle. These signals may coincide on some time intervals as we shall see later.

Throughout the paper Ω denotes the complement of

I = k≥0 I J k k . The Lebesgue measure of Ω, denoted λ[Ω]
, is assumed infinite. Consider x(•) the state of the closed-loop system in (1) with some feedback controller [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF]): The closed loop system is called weakly stable if for each ǫ > 0 there exists

U (X, Ẋ, X * d , Ẋ * d , Ẍ * d ). Definition 2 (Weakly Stable System
δ(ǫ) > 0 such that ||x(0)|| ≤ δ(ǫ) ⇒ ||x(t)|| ≤ ǫ for all t ≥ 0, t ∈ Ω.
The system is asymptotically weakly stable if it is weakly stable and lim t∈Ω, t→∞

x(t) = 0. Finally, the practical weak stability holds if there exists

0 < R < +∞ and t * < +∞ such that ||x(t)|| < R for all t > t * , t ∈ Ω. Consider I J k k ∆ = [τ k 0 , t k f ] and V (•)
such that there exists strictly increasing functions α(•) and β(•) satisfying the conditions: α(0) = 0, β(0) = 0 and α(||x||) ≤ V (x, t) ≤ β(||x||).

In the sequel, we consider that for each cycle the sequence of impact instants t k ℓ has an accumulation point t k ∞ . We note that a finite accumulation period (i.e.

t k ∞ < +∞) implies that e n < 1 (e n = 1 ⇒ t k ∞ = +∞ see [3]
). The following is inspired from [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF], and will be used to study the stability of the system (1).

Proposition 1 (Weak Stability): Assume that the task admits the representation (4) and that

a) λ[I J k k ] < +∞, ∀k ∈ N, b) outside the impact accumulation phases [t k 0 , t k ∞ ] one has V (x(t), t) ≤ -γV (x(t), t) for some constant γ > 0, c) ℓ≥0 V (t k- ℓ+1 ) -V (t k+ ℓ ) ≤ K 1 V p1 (τ k 0 ), ∀ℓ ≥ 0 for some p 1 ≥ 0, K 1 ≥ 0, d) the system is initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, e) ℓ≥0 σ V (t k ℓ ) ≤ K 2 V p2 (τ k 0 )+ξ for some p 2 ≥ 0, K 2 ≥ 0 and ξ ≥ 0. If p = min{p 1 , p 2 } < 1 then V (τ k 0 ) ≤ δ(γ, ξ)
, where δ(γ, ξ) is a function that can be made arbitrarily small by increasing the value of γ. The system is practically weakly stable with R = α -1 (δ(γ, ξ)).

Remark 1: Since the Lyapunov function is exponentially decreasing on the Ω k phases, assumption (d) in Proposition 1 means that the system is initialized on Ω 0 sufficiently far from the moment when the trajectory X nc (•) leaves the admissible domain and is therefore not stringent. The practical stability is very useful because attaining asymptotic stability is not an easy task for the unilaterally constrained systems described by (1) especially when n ≥ 2 and M (q) is not a diagonal matrix (i.e. there are inertial couplings, which is the general case). Precisely, the practical weak stability is characterized by an "almost decreasing" Lyapunov function V (x(•).•) as shown in Figure 4.

III. CONTROLLER DESIGN

In order to overcome some difficulties that can appear in the controller definition, the dynamical equations (1) will be expressed in the generalized coordinates introduced by McClamroch & Wang [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. We suppose that the generalized coordinates transformation holds globally in Φ, which may obviously not be the case in general. However, the study of the singularities that might be generated by the coordinates transformation is out of the scope of this paper. Let us

consider D = [I m . . . O] ∈ R m×n , I m ∈ R m×m the identity matrix. The new coordinates will be q = Q(X) ∈ R n , with q = q 1 q 2 , q 1 =    q 1 1 . . . q m 1    such that Φ = {q | Dq ≥ 0}.
The controller used here consists of different low-level control laws for each phase of the system. More precisely, the switching controller can be expressed as

T (q)U =    U nc for t ∈ Ω ∅ 2k U J t for t ∈ I J k U J c for t ∈ Ω J k ( 5
)
where T (q) = T 1 (q) T 2 (q) ∈ R n×n is full-rank under some assumptions (see [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]). The dynamics becomes:

8 > < > :
M11(q)q1 + M12(q)q2 + C1(q, q) q + g1(q) = T1(q)U + λ M21(q)q1 + M22(q)q2 + C2(q, q) q + g2(q) = T2(q)U q i 1 ≥ 0, q i 1 λi = 0, λi ≥ 0, 1 ≤ i ≤ m Collision rule [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF] where the set of complementary relations can be written more compactly as 0 ≤ λ ⊥ Dq ≥ 0. In the sequel U nc coincides with the fixed-parameter controller proposed in [START_REF] Slotine | Adaptive manipulator control: A case study[END_REF]. First, let us introduce some notations: q = qq d , q = qq * d , s = q + γ 2 q, s = q + γ 2 q, qr = qdγ 2 q where γ 2 > 0 is a scalar gain and q d , q * d represent the desired trajectories defined in the previous section. Using the above notations the controller is given by T (q)U 8 > > < > > :

Unc = M (q)qr + C(q, q) qr + G(q) -γ1s U J t = U J c , t ≤ t k 0 U J t = M (q)qr + C(q, q) qr + G(q) -γ1s, , t > t k 0 U J c = Unc -P d + K f (Pq -P d ) (7)
where γ 1 > 0 is a scalar gain, K f > 0, P q = D T λ and P d = D T λ d is the desired contact force during persistently constrained motion. It is clear that during Ω J k not all the constraints are active and, therefore, some components of λ and λ d are zero.

In order to prove the stability of the closed-loop system (6)- [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] we will use the following positive definite function:

V (t, s, q) = 1 2 s T M (q)s + γ 1 γ 2 qT q (8)
IV. TRACKING CONTROL FRAMEWORK In this paper we treat the tracking control problem for the closed-loop dynamical system ( 5)-( 7) with the complete desired path a priori taking into account the complementarity conditions and the impacts. In order to define the desired trajectory let us consider the motion of a virtual and unconstrained particle perfectly following a trajectory (represented by X nc (•) on Figure 1) with an orbit that leaves the admissible domain for a given period. Therefore, the orbit of the virtual particle can be split into two parts, one of them belonging to the admissible domain (inner part) and the other one outside the admissible domain (outer part). In the sequel we deal with the tracking control strategy when the desired trajectory is constructed such that:

(i) when no activated constraints, it coincides with the trajectory of the virtual particle, (ii) when p activated constraints, its orbit coincides with the projection of the outer part of the virtual particle's orbit on the surface of codimension p defined by the activated constraints (X d between A ′′ and C in Figure 1), (iii) the desired detachment moment and the moment when the virtual particle re-enters the admissible domain (w.r.t. p ≤ m constraints) are synchronized. Therefore we have not only to track a desired path but also to impose a desired velocity allowing the motion synchronization on the admissible domain.

The main difficulties here consist of:

• stabilizing the system on ∂Φ during the transition phases I J k k and incorporating the velocity jumps in the overall stability analysis;

• deactivating some constraints at the moment when the unconstrained trajectory re-enters the admissible domain with respect of them; • maintaining a constraint movement between the moment when the system was stabilized on ∂Φ and the detachment moment.

A. Design of the desired trajectories

Throughout the paper we consider

I J k k = [τ k 0 , t k f ],
where τ k 0 is chosen by the designer as the start of the transition phase

I J k k and t k f is the end of I J k k .
We note that all superscripts (•) k will refer to the cycle k of the system motion. We also use the following notations:

• t k 0 is the first impact during the cycle k, • t k ∞ is the accumulation point of the sequence {t k ℓ } ℓ≥0 of the impact instants during the cycle k (t k f ≥ t k ∞ ), • τ k 1 will be explicitly defined later and represents the instant when the desired signal X * d reaches a given value chosen by the designer in order to impose a closed-loop dynamics with impacts during the transition phases, • t k d is the desired detachment instant, therefore the phases

Ω J ′ k 2k+1 can be expressed as [t k f , t k d ].
It is noteworthy that t k 0 , t k ∞ , t k d are state-dependent whereas τ k 1 and τ k 0 are exogenous and imposed by the designer. To better understand the definition of these specific instants we present in Figure 1 a simplified representation of the signals introduced in the previous section.

A A'' B A' C Φ ∂Φ X * d (t) = X d (t) X d (t) X * d (t) X nc (t) = X * d (t) = X d (t) X nc (t)
Fig. 1. The closed-loop desired trajectory and control signals

The points A, A ′ , A ′′ and C in Figure 1 correspond to the moments τ k 0 , t k 0 , t k f and t k d respectively. On the other hand in Figure 1 we see that starting from A the desired trajectory

X d (•) = X * d (•) is deformed compared to X nc (•).
In order to reduce this deformation τ k 0 and implicitly the point A must be close to ∂Φ. Taking into account just the constraints J ′ k \ J k+1 we can identify t k d with the moment when X d (•) and X nc (•) rejoin at C.

B. Design of q *

d (•) and q d (•) on the phases

I J k k
During the transition phases the system must be stabilized on ∂Φ. Obviously, this does not mean that all the constraints have to be activated. Let us consider that only the first p constraints (eventually reordering the coordinates) define the border of Φ where the system must be stabilized. On [τ k 0 , t k 0 ) we define q * d (•) as a twice differentiable signal such that q * d (t) approaches a given point of the normal cone N Φ (q p = 0) on

[τ k 0 , τ k 1 ]
. Precisely, we define q * d (•) such as:

q * d (τ k 0 ) = q nc (τ k 0 ), q * d (τ k 0 ) = q nc (τ k 0 ) q * d (τ k 0 + δ) = q nc (τ k 0 ), q * d (τ k 0 + δ) = 0 (q i d ) * (τ k 1 ) = -ϕV 1/3 (τ k 0 ), ( q * d ) p (τ k 1 ) = 0, i = 1, . . . , p (9)
where ϕ > 0, and δ > 0 is a small constant introduced in order to assure the twice differentiability of q * d before the first impact. The rationale behind the choice of q * d (•) is on one hand to assure a robust stabilization on ∂Φ, mimicking the bouncing-ball dynamics; on the other hand to enable one to compute suitable upper-bounds that will help using Proposition 1 (hence V 1/3 (•) terms in [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF] with V (•) in ( 8)).

Remark 2: Two different situations are possible. The first is given by t k 0 > τ k 1 (see Figure 2) and we shall prove that in this situation all the jumps of the Lyapunov function in [START_REF] Glocker | Concepts for Modeling Impacts without Friction[END_REF] are negative. The second situation was pointed out in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] and is given by t k 0 < τ k 1 . In this situation the first jump at t k 0 in the variation of the Lyapunov function is positive, therefore the system can be only weakly stable.

B C A t t k 0 t k 1 t k d (q * d ) i (t) q i 1 q i 1 (t) τ k 0 t k 0 > τ k 1 Ω 2k I k Ω 2k+1 Ω 2k+2 A ′ t k f τ k 1 (q * d ) i (t) (q * d ) i (t) -ϕV 1/3 (τ k 0 )
Fig. 2. The design of q * 1d on the transition phases

I J k k
During the transition phases

I J k k we define (q d ) n-p (t) = (q * d ) n-p (t).
Assuming a finite accumulation period, the impact process can be considered in some way equivalent to a plastic impact. Therefore, (q d ) p and ( qd ) p are set to zero on the right of t k 0 .

V. DESIGN OF THE DESIRED CONTACT FORCE DURING CONSTRAINT PHASES

For the sake of simplicity we consider the case of the constraint phase Ω J k , J = ∅ with J = {1, . . . , p}. Obviously a sufficiently large desired contact force P d assures a constrained movement on Ω J k . Nevertheless at the end of the Ω J 2k+1 phases a detachment from some surfaces Σ i has to take place. It is clear that a take-off implies not only a well defined desired trajectory but also some small values of the corresponding contact force components. On the other hand, if the components of the desired force decrease too much a detachment can take place before the end of the Ω J 2k+1 phases which can generate other impacts. Therefore we need a lower bound of the desired force which assures the contact during the Ω J k phases. Dropping the time argument, the dynamics of the system on Ω J k can be written as

M (q)q + F (q, q) = U c + D T p λ p 0 ≤ q p ⊥ λ p ≥ 0 (10)
where F (q, q) = C(q, q) q+G(q) and D p = [I p . . . O] ∈ R p×n .

On Ω J k the system is permanently constrained which implies q p (•) = 0 and qp (•) = 0. In order to assure these conditions it is sufficient to have λ p > 0. In the following let us denote

M -1 (q) = [M -1 (q)] p,p [M -1 (q)] p,n-p [M -1 (q)] n-p,p [M -1 (q)] n-p,n-p and
C(q, q) = C(q, q) p,p C(q, q) p,n-p C(q, q) n-p,p C(q, q) n-p,n-p where the meaning of each component is obvious. Proposition 2: On Ω J k the constraint motion of the closedloop system (10),( 5),( 7) is assured if the desired force is defined by

(λ d ) p β - Mp,p (q) 1 + K f [M -1 (q)] p,p C p,n-p (q, q)+ [M -1 (q)] p,n-p C n-p,n-p (q, q) + γ 1 [M -1 (q)] p,n-p s n-p (11)
where Mp,p (q) = [M -1 (q)] p,p -1 and β ∈ R p , β > 0.

Proof: For the sake of brevity we give here only the idea of the proof. The result is based on the solution of the LCP derived combining ( 10) and [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF], which is:

0 ≤ D p M -1 (q) -F (q, q) + U nc -(1 + K f )D T p (λ d ) p +(1 + K f )D p M -1 (q)D T p λ p ⊥ λ p ≥ 0 (12)
It is worth to precise that the LCP ( 12) has a unique solution since

(1 + K f )D p M -1 (q)D T p > 0.

VI. STRATEGY FOR TAKE-OFF AT THE END OF CONSTRAINT PHASES Ω J 2k+1

In this section we are interested in finding the conditions on the control signal U J c that assure the take-off at the end of constraint phases Ω J 2k+1 . As we have already seen before, the phase Ω J 2k+1 can be expressed as the time interval

[t k f , t k d ). The dynamics on [t k f , t k d
) is given by (10) and the system is permanently constrained, which implies q p (•) = 0 and qp (•) = 0. Let us also consider that the first r constraints (r < p) have to be deactivated. Thus, the detachment takes place at t k d if qr (t k+ d ) > 0 which requires λ r (t k- d ) = 0. The last pr constraints remain active which means λ p-r (t k- d ) > 0. To simplify the notation we drop the time argument in many equations of this section. We denote the LCP matrix as:

(1 + K f )D p M -1 (q)D T p = A 1 (q) A 2 (q) A 2 (q) T A 3 (q) with A 1 ∈ R r×r , A 2 ∈ R r×(p-r) and A 3 ∈ R (p-r)×(p-r) .
Proposition 3: For the closed-loop system (10),( 5),( 7) the decrease of active constraints from p to r < p, is possible if at the instant

t k d (λ d ) r = " A1 -A2A -1 3 A T 2 " -1 `br -A2A -1 3 bp-r -A2C2 ´-C1 (λ d ) p-r = C2 + A -1 3 " bp-r -A T 2 (λ d ) r " (13)
where

b p = b(q, q, U nc ) D p M -1 (q)[U nc -F (q, q)] ≥ 0 and C 1 ∈ R r , C 2 ∈ R p-r such that C 1 ≥ 0, C 2 > 0 Proof:
The result follows solving the LCP [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. Proposition 4: The closed-loop system (10),( 5), [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] 13) where

C 1 is replaced by C 1 (t -t k d ). (ii) On [t k d , t k d + ǫ) q * d (t) = q d (t) = q * r (t) q nc n-r (t)
, where q * r (t) is a twice differentiable function such that q

* r (t k d ) = 0, q * r (t k d + ǫ) = q nc r (t k d + ǫ), q * r (t k d ) = 0, q * r (t k d + ǫ) = qnc r (t k d + ǫ) (14)
and q * r (t k d ) = a > max 0, -A 1 (q) ((λ d ) r (t k- d ) . Proof: The uniqueness of solution of the LCP [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] assures us that [START_REF] Morȃrescu | Passivity-based tracking control of multiconstraint complementarity Lagrangian systems[END_REF] and ( 13) agree if C 1 < 0. In other words, replacing C 1 by C 1 (tt k d ) in [START_REF] Slotine | Adaptive manipulator control: A case study[END_REF] we assure a constrained motion on [t k f , t k d ) and the necessary conditions for detachment on [t k d , t k d + ǫ). The second item simply says that the detachment is assured by a sufficiently large positive desired acceleration q * r . Relations ( 14) assure the twice differentiability of q d and q * d .

VII. CLOSED-LOOP STABILITY ANALYSIS

In the case Φ = R n , the function V (t, s, q) in ( 8) can be used to prove the closed-loop stability of the system ( 6), [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] (see for instance [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF]). In the case studied here (Φ ⊂ R n ) the analysis becomes more complicated as shown in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF].

To simplify the notation V (t, s(t), q(t)) is denoted as V (t). In order to introduce the main result of this paper we make the next assumption, which is verified in practice for dissipative systems.

Assumption 1: The controller U t in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] assures that the sequence {t k ℓ } ℓ≥0 of the impact times possesses a finite accumulation point t k ∞ i.e. lim ℓ→∞ t k ℓ = t k ∞ < +∞, ∀k. Theorem 1: Let Assumption 1 hold, e n ∈ [0, 1) and (q * d ) p defined as in [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF]. The closed-loop system (5)-( 7) initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, satisfies the requirements of Proposition 1 and is therefore practically weakly stable with the closed-loop state x(•) = [s(•), q(•)] and R = Ke -γ(t k f -t k ∞ ) where γ = 2γ 1 /λ max (M (q)) and K > 0 is a real constant.

Proof: The proof consist of verifying the conditions b), c) and e) of Proposition 1. The details can be found in [START_REF] Morȃrescu | Passivity-based tracking control of multiconstraint complementarity Lagrangian systems[END_REF].

VIII. ILLUSTRATIVE EXAMPLE

The numerical simulations are done with the Moreau's time-stepping algorithm of the SICONOS software platform [START_REF] Acary | An introduction to Siconos[END_REF]. The choice of a time-stepping algorithm was mainly dictated by the presence of accumulations of impacts which render the use of event-driven methods difficult. The influence of different parameter as γ 1 , γ 2 , e n , τ k 0 or the time-step, is studied in [START_REF] Morȃrescu | Passivity-based tracking control of multiconstraint complementarity Lagrangian systems[END_REF] by simulating the behavior of a planar two-link rigid-joint manipulator in presence of one unilateral constraint.

Let us consider in the sequel a planar two-link rigid-joint manipulator with two constraints. Precisely we impose an admissible domain Φ = {(x, y) | y ≥ 0, 0.7-x ≥ 0}. Let us also consider an unconstrained desired trajectory X nc whose orbit is given by the circle {(x, y) | (x -0.7) 2 + y 2 = 0.5} that violates both constraints. In other words, the two-link planar manipulator must track a quarter-circle; stabilize on and then follow the line Σ 1 = {(x, y) | y = 0}; stabilize on the intersection of Σ 1 and Σ 2 = {(x, y) | x = 0.7}; detach from Σ 1 and follow Σ 2 until the unconstrained circle re-enters Φ and finally take-off from Σ 2 in order to repeat the previous steps. It is noteworthy that the task presented above is not of type (4) since after a constraint phase (when the end-effector is attached to Σ 1 ) follows a transition phase instead of a detachment. However, as we have pointed out in Section II the succession of phases is not determinant and the manipulator can accomplish the task under consideration in a weakly stable way.

Let us consider in this case that a cycle is

Ω 2k ∪ I 1 k ∪ Ω 1 2k+1 ∪ I 2 k ∪ Ω 2 2k+1
where Ω 2k is the free-motion phase and I i k , Ω i 2k+1 are the impacting transients and the constrained phases associated to the surface Σ i . The numerical values used for the dynamical model are again l 1 = l 2 = 0.5m, I 1 = I 2 = 1kg.m 2 , m 1 = m 2 = 1kg and the restitution coefficient e n = 0.7. We impose a period of 10 seconds for each cycle and we simulate the dynamics during 6 cycles. Setting the controller gains γ 1 = 15, γ 2 = 15 we see in Figure 3 that the desired trajectory is accurately followed. The same conclusion can be deduced looking at the variation of the Lyapunov function plotted in Figure 4. The simulation are done imposing a constant contactforce λ 1 during the motion on the surface Σ 1 (see Figure 5 (left)) and a decreasing contact-force, that allows a smooth detachment, during the motion on Σ 2 (see Figure 5

IX. CONCLUSIONS

In this paper we have proposed a methodology to study the tracking control of fully actuated Lagrangian systems subject to multiple frictionless unilateral constraints. The main contribution of the work is twofold: first, it formulates a general framework and second, it provides a complete stability analysis for the class of systems under consideration. It is noteworthy that even in the case of only one frictionless unilateral constraint the paper presents some improvements with respect to the existing works in the literature. Precisely, the stability analysis result is significantly more general than those presented in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] and [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF] and, each element entering the dynamics (desired trajectory, contact force) is explicitly defined. Numerical simulations for an illustrative example, using the SICONOS software platform [START_REF] Acary | An introduction to Siconos[END_REF] complete the study.

O

Fig. 3 .

 3 Fig. 3. The trajectory of the system during 6 cycles.

Fig. 4 .

 4 Fig. 4. Up: Lyapunov function-zoom on the transition phases during the first cycle; Down: The variation of the Lyapunov function on the last 5 cycles.

  (right)).

Fig. 5 .

 5 Fig. 5. Left: Variation of the contact force during the motion on Σ1; Right: Variation of the contact force during the motion on Σ2