Augustin-Alexandru Saucan

Thierry Chonavel
email: thierry.chonavel@telecom-bretagne.eu

Christophe Sintes
email: christophe.sintes@telecom-bretagne.eu

Jean-Marc Le Caillec
email: jm.lecaillec@telecom-bretagne.eu

Interacting Multiple Model Particle Filters for Side Scan Bathymetry

Keywords: side scan sonar, bathymetry, DOA estimation, tracking, bootstrap filter, multiple model, IMM, Markov jump systems

In this paper we propose a multiple sea floor model based approach to improve bathymetry estimation with tracking algorithms. Traditionally interferometry is used to estimate the phase difference of signals received by two sensors, implicitly the direction of arrival (DOA) of the wave impinging both sensors. In our approach, we employ a state space model to describe data collected by a multi-sensor side scan sonar, and the evolution of the underlying DOA angle. The challenge with space state models is choosing the right model, and detecting the switch between models. We propose the use of several models that describe different sea-floor patterns and merge them within the framework of the interacting multiple model (IMM). Since the sonar array processing problem is non-linear and non-Gaussian, we propose an IMM particle filter algorithm to provide robust tracking while not sacrificing performance. Also an interesting new application is the swath segmentation, which appears as a side result implied by calculating the different model probabilities.

I. INTRODUCTION

In a recent paper [START_REF] Saucan | Enhanced sonar bathymetry tracking in multi-path environment[END_REF] we successfully used a DOA tracking technique for sea-floor angle of arrival estimation in the case of side scan sonars. The tracking algorithm exploits noisy array data to recover the sea-floor profile and, more importantly, to resolve multiple echoes interfering with the main sea floor echo. The multiple interfering echoes are caused by the existence of multiple paths involving sea surface reflections. Moreover, due to the non-linear and non-Gaussian state space model, the implementation of the proposed tracker was based on a particle filter, more precisely the bootstrap filter [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF].

In all state space model descriptions a crucial role is played by the state equation, which represents an added prior upon the temporal evolution of the internal state. For example kynematic state models [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Algorithms and Software for Information Extraction[END_REF] such as the nearly constant velocity (NCV), wiener acceleration process (WAP) or coordinated turn (CT) models are employed when tracking dynamical targets. However, when dealing with maneuvering targets often one state model isn't sufficient and several multiple model (MM) algorithms where developed, the most popular being the IMM, first proposed in [START_REF] Blom | The interacting multiple model algorithm for systems with markovian switching coefficients[END_REF]. In our previous article, since interfering echoes were weakly structured, the internal state was chosen to be the DOA of the sea floor echo. However the problem with such an approach is that one state model cannot describe all possible terrain reliefs and corresponding angle variation. In this sense, this paper aims at presenting the adaptation and application of an IMM Bootstrap filter in the case of bathymetry tracking. The proposed algorithm is based on the IMM particle filter introduced in [START_REF] Boers | Interacting multiple model particle filter[END_REF]; several batymetrical state models are employed on simulated signals and real side-scan sonar data. Also the choice of the state variable between DOA angle and bathymetric height is presented.

II. MARKOVIAN JUMP MODEL FOR BATHYMETRY

TRACKING

We consider a dynamic system with the internal state represented by the stochastic process x t , representing either the sea floor DOA angle θ t or the bathymetric height h t , both being linked by a non-linear transformation. The scope is to estimate the inner state based on the system observations, represented by the stochastic process y t , i.e. the signal received by the sensor array. let us consider equations 1 and 2, that represent the state and measurement equations of the model:

x t = F t (x t-1 , v t , m t ) (1) y t = H t (x t , n t , m t ) (2)
[Π] ij = π ij P (m t = j|m t-1 = i), ∀i, j ∈ {1, . . . , M } (3)
In the context of array processing, the observation is usualy linked to the DOA angle θ t (also indirectly to the bathymetric height h t ) in the following way:

y t = A(θ t )s t + n t (4)
where we considered one narrow-band source located in the far-field of the uniform linear array (ULA) comprised of P receivers. Ideally this source signal will represent the sea floor back scattered signal. y t ∈ C [P ×1] represents the signal received over the array. s t ∈ C represents the source signal, n t ∈ C [P ×1] represents the additive noise component and A ∈ C [P ×1] represents the steering vector, which is a function of the internal state. In the case of a well calibrated array, the steering vector is of the form:

A(θ t ) =      1 exp(i2π d λ sin(θ t )) . . . exp(i2π d λ (P -1) sin(θ t )      (5)
where θ represents the direction of arrival of the source i.e. the back scattered sea-floor echo, d the receiver spacing and λ the wavelength. Note that θ t is the angle between the echo DOA and the plane orthogonal to the array axis. In [START_REF] Saucan | Enhanced sonar bathymetry tracking in multi-path environment[END_REF] it was shown that the observation distribution closely obeys a multivariate Laplace law [START_REF] Eltoft | On the multivariate laplace distribution[END_REF], thus in this article we will utilize this multivariate probability distribution. Since we know the distribution of the whole observation p(y t ) and not just the noise distribution p(v t ) we are able to avoid the nuisance parameter s t estimation. Furthermore the Laplace distribution proves to be robust in the sense that it yields a good overlap between the prediction pdf and the observation pdf thus making possible the application of the bootstrap filter.

The multiple model framework that we propose for bathymetry tracking is described by the jump Markov system state equation 1 and the measurement equation 4. The model state m t is used to switch between different state prediction models either for the DOA angle or directly for the bathymetric height. The measurement equation is unaffected by m t , however in the case of bathymetric height driven state equations an additional transformation is required.

It is well known [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Algorithms and Software for Information Extraction[END_REF], [START_REF] Ristic | Beyond the Kalman Filter: Particle filters for Tracking Applications[END_REF] that the optimal filtering algorithm in the switching multiple model framework, that is, the algorithm that recursively computes the posterior distribution p(x t |y 1:t ), requires a number of mode matched filters that exponentially increases with t. This branching of modes phenomenon renders impractical the optimal filter. Thus, several sub-optimal methods have been proposed. Arguably the most successful sub-optimal filter, the IMM has a distinctive similarity to the optimal filter while imposing a constant number of mode matched filters throughout the filtering process, it has M filters running in parallel [START_REF] Bar-Shalom | IMM estimator versus optimal estimator for hybrid systems[END_REF]. If all the modes in the Markovian switching system are linear and Gaussian, each of these mode matched filters is a Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. However, since the bathymetry tracking system has a non-linear and non-Gaussian measurement equation 4, we propose to use a particle implementation of the mode matched filters. The first practical particle filter, the bootstrap (BF) was proposed in [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF] and approximates the posterior p(x t |y 1:t ) distribution with a weighted set of particles. For a thorough introduction into particle filtering methods see for instance [START_REF] Candy | Bayesian Signal Processing: Classical, Modern and Particle Filtering Methods[END_REF].

III. IMM BF FOR BATHYMETRY TRACKING

In this section we present the proposed IMM based particle filter aimed at estimating the internal state of the bathymetric jump Markov system. For a system comprised of M possible models, the IMM is composed of M parallel mode matched filters, each for one possible model. The probabilities of each model, p(m t |y 1:t ) are recursively computed at each time step and used to merge the individual filter estimates. The mode matched filters are implemented using a bootstrap filter, with each filter having its own set of particles, as opposed to the single shared set implementation proposed in [START_REF] Mcginnity | Multiple model bootstrap filter for maneuvering target tracking[END_REF], also presented into more details in [START_REF] Doucet | Sequantial Monte Carlo Methods in Practice[END_REF]. The IMM algorithm is generically decomposed into 4 stages: mixing, model conditional filtering, model probability update and merging of individual filter estimates. The proposed IMMBF algorithm is illustrated in Algorithm 1.

The key step of the IMM algorithm is the mixing stage, where the input of each mode matched filter j at time instant t is obtained from the mixing of the mode matched filter estimates from time instant t -1, effectively achieving an interaction of the M models. The weightings used to perform the mixing are called mixing probabilities, µ i|j and are defined as the probability that mode i was in effect at time t -1 given that model j is in effect at time t. In the case of linear-Gaussian systems the mixing distribution p(x t |m t = j, y 1:t ) is a weighted sum of Gaussian distributions with the weights µ i|j [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Algorithms and Software for Information Extraction[END_REF]. Note that the original IMM approximates the mixing distribution by a single Gaussian with the same mean and covariance as the Gaussian mixture. When dealing with nonlinear and non-Gaussian systems, Monte Carlo approximations of distributions as weighted particles are employed [START_REF] Doucet | Monte carlo methods for signal processing: A review in the statistical signal processing context[END_REF]. In such conditions, the mixing distribution for m t = j is obtained by sampling from the µ i|j weighted sum of model matched posterior densities obtained at t -1, i µ i|j p(x t-1 |m t-1 = i, y 1:t ).

The mixing distribution computed for each model j is utilized to reinitialize the mode matched particle filters in the model conditioned filtering step. Furthermore filtering is carried out in the usual two step procedure: prediction of the reinitialized particles through the state transition equation 1 and weights updating of the reinitialized particles taking into account the new observation. We only perform resampling when effective number of particles drops below a threshold value [START_REF] Liu | Sequential monte carlo methods for dynamical systems[END_REF]. At the end of this step the mode conditioned posterior distribution p(x t |m t = j, y 1:t ) is obtained for each mode j as a weighted particle approximation.

The model conditioned likelihoods p(y t |x

(n)j t

) computed for each particle (n) and model j are utilized to update the model probabilities. As suggested in [START_REF] Driessen | Efficient particle filter for jump markov nonlinear systems[END_REF] we consider an average model likelihood L j t as the mean of the individual particles likelihoods in each of the models. Based on the newly computed model probabilities the overall posterior density p(x t |y 1:t ) is estimated as the weighted sum of the model conditioned posterior distributions p(x t |m t = j, y 1:t ). This operation is again achieved by sampling.

P (m t = j|y 1:t-1 ) = i π ij µ i t-1
Mixing weight:

µ i|j P (m t-1 = i|m t = j, y 1:t-1 ) = π ij µ i t-1 /µ j t|t-1
Mixing density:

p(x t-1 |m t = j, y 1:t-1 ) = i µ i|j p(x t-1 |m t-1 = i, y 1:t-1 ) Sample {(x (n) t-1|mt=j , w(x (n) t-1|mt=j
)) n=1,...,N } ∼ p(x t-1 |m t = j, y 1:t-1 )

Model conditional filtering

for n = 1 to N do Particle propagation:

x(n)j t ∼ p(x t |x (n) t-1|mt=j , m t = j) Particle likelihood: q (n)j t = p(y t |x (n)j t ) Particle weight update: w(x (n)j t ) = q (n)j t w(x (n) t-1|mt=j ) end Compute mode likelihood: L j t = 1 N n q (n)j t
Normalize particle weights:

w(n)j = w(x (n)j t ) n w(x (n)j t )
Compute effective number of samples:

N j ef f = 1 n w(n)j if N j ef f < N thresh then Resample {(x (n)j t , w(n)j ) n=1,...,N } to obtain {(x (n)j t , w (n)j = 1 N ) n=1,...,N } else {(x (n)j t , w (n)j ) n=1,...,N } = {(x (n)j t
, w(n)j ) n=1,...,N } end end 3. Model Probability update for j = 1 to M do Mode probability:

µ j t|t P (m t = j|Y 1:t ) = µ j t|t-1 L j t M i=1 µ i t|t-1 L i t end 4.

Model merging

Overall MMSE estimate:

xt|t E{x t |y 1:t } ≈ M j=1 µ j t N n=1 x (n)j t w (n)j

IV. RESULTS: SIMULATION AND REAL SONAR DATA

A. Simulated underwater scene

The simulated underwater scenario is used to assess the ability of the IMMBF algorithm to detect jumps, i.e. model changes in the DOA angle evolution and to track the bathymetry. Thus two models for the bathymetric height were envisaged: a random walk model, RW and a nearly constant slope NCS, defined in a similar fashion as the nearly constant velocity model for dynamical targets. The idea behind the two is that the NCS model is intended to describe flat smooth regions with unknown inclination angles and the RW model is intended for rough, i.e. non-smooth areas. The IMMBF consisted of N = 100 particles in each model. The initial model probabilities and the transition matrix are given by:

µ 0 = 0.5 0.5 Π = 0.98 0.02 0.02 0.98 . ( 6
)
The state variable x t is composed of the bathymetric height and it's derivative with respect to the slant range r S . ∆ S = cT s /2 represents the range bin width, T s the sampling period and c the celerity.

x t = h t ∂ht ∂r S T (7)
The model transition functions are given by:

F (x t-1 , v t , 1) = F RW x t-1 + v RW t (8) F (x t-1 , v t , 2) = F N CS x t-1 + Γv N CS t (9) F N CS =   1 ∆ S 0 ∆ 2 S 2   , Γ =   ∆ 2 S 2 ∆ S   , F RW = 1 0 0 0 (10)
v RW t and v N CS t are white zero mean Gaussian sequences with standard deviations σ RW and σ N CS and represent a bathymetric height noise and a height "acceleration" noise sequence. The simulated underwater scene is presented in Fig. 1 and represents the swath bathymetric profile of one ping line. The tilted sonar array is located 15m above the sea-floor, which exhibits a flat and smooth surface except for a rough patch, a fluctuation of the sea-floor height. Since the scenario is intended only for the switching of models, there is no shadow effects simulated on the backscattered signal. However noise is added to the backscattered signal and the SN R, signal to noise power ratio is considered to drop from 30dB to 10dB along the swath width accounting for the wave spherical spreading effect. Also in Fig. 1 we present in detail the rough patch and also the estimates obtained with the IMMBF algorithm.

In Fig. 2 we compare the angular pseudo-spectrum obtained with MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], the sea-floor profile and the estimated probabilities of both models, RW and NCS, obtained with the IMMBF algorithm. The effect of the rough patch on the DOA of the backscattered wave can be seen in a detail of the pseudospectrum, however, affected by the averaging of the MUSIC algorithm. As expected, we observe that the model active in the flat regions of the sea-floor is the NCS model, while the RW is active in the rough patch area. The IMMBF algorithm successfully detects the switch (jump) of the model best suited to describe the sea-floor height evolution. The variances were set to σ RW = 1m and σ N CS = 1m -2 . Although some design considerations should be taken into account when choosing the value of these variances, the detection of the model jumps was found to be robust for values within the chosen order. The ability of the IMMBF algorithm to detect model jumps is not only important for improving the estimates by using the locally best adapted model but considering the model probability vector also provides an information about the seafloor's nature: segmenting the swat into the classes determined by the models.

Similar results were obtained, in terms of model jump detection, with two random walk sea-floor models with standard deviations: σ RW 1 = 0.01m and σ RW 1 = 1m. The first model describes a smoother surface while the second one describes a rougher surface; coupled with the IMMBF they achieve basically the same segmentation, as shown in Fig. 2.

B. Real side-scan sonar data

The IMMBF algorithm was also tested on data issued from a bathymetric side-scan sonar. The sonar receiving arrays are composed of 8 staves allowing for bathymetry reconstruction as well as classical imagery. The data was collected in a survey of a shallow water harbor, with a relatively flat floor at depth varying around 15m. This particular combination of shallow water, calm sea and pole mounted sonar causes the apparition of multi-path echoes in the angular spectrogram, as shown in Fig. 3 and reference [START_REF] Saucan | Enhanced sonar bathymetry tracking in multi-path environment[END_REF]. The spectrogram in Fig. 3 was obtained with the MVDR or Capon beamformer [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF], and superimposed on this image are the IMMBF estimates. Let us notice the similarity of the DOA arch-like curve similar with the one in Fig. 2 that is typical of a relatively flat seabottom. We used a two model Markovian switching system with model transition matrix and initial probabilities given by 6. For the implementation, we used again N = 100 particles in each model. However both models are angular random walk models with different model noises.

x t = θ t (11) F (x t-1 , v t , m t ∈ {1, 2}) = x t-1 + v mt t ( 12
)
where θ t represents the DOA angle, v mt t is white zero mean Gaussian noise sequences with variance σ mt . The first model has a constant standard deviation σ 1 = 0.1 • and the second one has a time varying standard deviation σ 2 = f (t), starting with 0.1 • at the beginning of the swath and increasing linearly to 2 • at the end of the swath width.

In Fig. 4 we present the IMMBF estimated probabilities for both models. At the beginning of the swath width the second model is better adapted to describe the evolution of the bathymetry. Indeed the DOA angle has a faster changing rate at the beginning of the swath than at the end of the swath and the first model with its small variability does not cope well at small range values. However, it offers better estimates at bigger ranges. To better illustrate the model choice, in Fig. 5 we present the particle clouds of both models. In each mode and at each range bin, there are 100 particles, i.e. angular positions, which together with their associated weights form the empirical posterior density. Although the particle positions are not enough to characterize the variance of particles inside each mode, they offer a visual interpretation of the fluctuation in each mode. Also in Fig. 5 around the 35m range mark we can observe a sudden increase in particle fluctuation in both of the models. This is caused by a discontinuity of the DOA curve and its effect is also visible in Fig. 4 with a shorttransient jump in the model probabilities. Such discontinuities, generally caused by shadows, correspond to an absence of signal and need to be detected beforehand since tracking should not be conducted over shadows.

V. CONCLUSION AND PERSPECTIVES

This paper presents the application of the IMMBF algorithm to simulated and real side-scan sonar data. At the price of increased complexity, multiple model algorithms offer the possibility of robust tracking when the models used to track the dynamics of the target are relevant. The simulated and real data experiments show the potential interest of multiple model algorithms for bathymetry enhancement. Also a new and interesting application, which will be further investigated is the swath segmentation.

VI. ACKNOWLEDGMENT

This work was funded by a grant of Direction Générale d'Armement (French MoD) and the Carnot Institute.

Algorithm 1 :

 1 IMM Bootstrap Filter (IMMBF).

for j = 1 to M do 1 .

 1 Mixing or model conditional reinitializationPredicted mode probability: µ j t|t-1

Fig. 1 .

 1 Fig. 1. Simulated underwater scene and detail of rough patch.

Fig. 2 .

 2 Fig. 2. Comparison of angular spectrum, sea-floor profile and IMMBF model probabilities.

Fig. 3 .

 3 Fig. 3. Capon angular spectrum and IMMBF DOA estimates.

Fig. 4 .

 4 Fig. 4. IMMBF model probabilities.

Fig. 5 .

 5 Fig. 5. Particle clouds for both models.

  F t and H t represent the model conditioned transition and measurement functions and are possibly non-linear. x t together with m t define a hybrid space state system composed of euclidean and discreet valued variables. When the modal state is a Hiden Markov Chain the model is called Markovian jump system. The stochastic processes v and n are mutually independent white noises with known probability densities. v represents the model noise and n the measurement noise.The Markov chain m t dictates the model active during the time interval ]t -1, t], with the transition probability matrix Π defined as:

where t denotes the discrete time and m t , the modal state (model) which is a time-homogeneous (hidden) Markov chain with M states (models).


Khan Le Tran

Cyril Douthe
email: cyril.douthe@ifsttar.fr

Karam Sab

Julien Dallot

Laurence Davaine

K L Tran

A preliminary design formula for the strength of stiffened curved panels by design of experiment method

Keywords: Design of computer experiments, Response surface, Cylindrical curved panels, Stiffeners, Stability, GMNIA

In bridge construction, the use of stiffened plates for box-girder or steel beams is common day to day practice. The advantages of the stiffening from the economical and mechanical points of view are unanimously recognized. For curved steel panels, however, applications are more recent and the literature on their mechanical behaviour including the influence of stiffeners is therefore limited. Their design with actual finite element software is possible but significantly time-consuming and this reduces the number of parameters which can be investigated to optimise each panel. The present paper is thus dedicated to the development of a preliminary design formula for the determination of the ultimate strength of stiffened cylindrical steel panels. This approximate formula is developed with help of a design of experiment method which has been adapted from the current statistical knowledge. This method is first presented and its feasibility as well as its efficiency are illustrated through an

Introduction

The interest of stiffening steel plates or panels to increase their strength under compression has been known for almost a century [START_REF] Lundquist | Comparison of three methods for calculating the compressive strength of flat and slightly curved sheet and stiffener combinations[END_REF]. In the field of structural engineering, the use of such panels is a common practice, for example in bottom flanges of box-girder bridges. Recent developments of the curving process allowed for the use of curved panels in civil engineering structures where they offer attractive aesthetic and aerodynamic possibilities.

The verification of these panels is yet difficult due to a lack of specifications, especially in European Standards: EN 1993-1-5 [START_REF]EN1993-1-5 . Eurocode 3 design of steel structures part 1-5: Plated structural elements[END_REF] gives specifications for flat or slightly curved panels with the condition R ≥ R lim = b 2 /t p (where R is the curvature radius of the panel, b its width and t p its thickness) and EN 1993-1-6 [START_REF]Eurocode 3 design of steel structures part 1-6: Strength and stability of shell structures[END_REF] deals only with revolution cylindrical shells. Nevertheless the curved panels in bridges have characteristics exactly between these two conditions, as illustrated in the case of the Confluences bridge in Angers, France 2011 (Fig. 1), whose radius R = 80 m is much smaller than the limit of EN 1993-1-5: R lim = 1440 m (with b = 4.8 m and t = 16 mm) and for which EN 1993-1-6 is not applicable neither because these curved flanges are not full revolution cylinders. From an academic point of view, the articles related to the buckling theory of curved panels are not so numerous due to the complexity of the studied problem and also due to its late application in the bridge construction. First investigations were conducted in the forties by Batdorf & Schildcrout [START_REF] Batdorf | Critical axial-compressive stress of a curved rectangular panel with a central chordwise stiffener[END_REF] and Schildcrout & Stein [START_REF] Schildcrout | Critical axial-compressive stress of a curved rectangular panel with a central longitudinal stiffener[END_REF] who showed that the stiffeners and the curvature increase the critical buckling strength. A state of the art on curved stiffened panels was then proposed by Becker [START_REF] Becker | Handbook of structural stability. Part VI: Strength of stiffened curved plates ans shells[END_REF] in 1958 in its handbook on structural stability. Based on experimental results (provided by Gall [START_REF] Gall | Compressive strength of stiffened sheets of aluminum alloy[END_REF], Lundquist [START_REF] Lundquist | Comparison of three methods for calculating the compressive strength of flat and slightly curved sheet and stiffener combinations[END_REF] and Ramberge et al. [START_REF] Ramberg | Effect of curvature on strength of axially loaded sheet-stringer panels[END_REF]), he confirmed that, when a stiffened flat panel is bent to a circular curve, its buckling stress is slightly increased (around 6% for the tested specimen which is relatively few compared to the effect of stiffening alone or curvature alone). More recent parametric studies based on numerical examples and the finite element modelling (e.g. Cho et al. [START_REF] Cho | Experimental and numerical investigations on the ultimate strength of curved stiffened plates[END_REF],

Khedmati & Edalat [START_REF] Khedmati | A numerical investigation into the effects of parabolic curvature on the buckling strength and behaviour of stiffened plates under in-plane compression[END_REF] or Park et al. [START_REF] Park | Estimation of buckling and collapse behaviours of stiffened curved plates under compressive load[END_REF]) investigated and quantified the influence of the main parameters on the ultimate strength of curved stiffened plates. They however did not lead to a practical criterion for the evaluation of the resistance of such panels which is therefore still an open question.

In a former study, the authors [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF] had investigated the case of unstiffened cylindrical curved panels under axial compression and established a set of formulas for the evaluation of the ultimate strength (which were confirmed by [START_REF] Martins | Eigenvalue analysis of cylindrically curved panels under compressive stresses -extension of rules from {EN} 1993-1-5[END_REF]). These semi-analytic formulas had been fitted on a total of 524 combinations of the main parameters. Each calculus involved Geometrical and Material Non-linearity with Imperfection Analysis (GMNIA) and required between 5 and 10 minutes depending on the refinement of the mesh. Considering the fact that in the case of stiffened panels the number of parameters is considerably larger, re-employing the same methodology seemed unrealistic. It appeared hence that there is a need for a robust strategy for the choice of the set of tested models and for the measure of the approximated model accuracy. Such a strategy exists for the design of physical experiments as well as for that of computer experiments, they are known as "design of experiments methods".

In the following, the authors present thus first the characteristics of computer experiment strategies. Afterwards the feasibility and ease of use of the methodology as well as its efficiency are illustrated through an application to the reference case of unstiffened curved panels. Then, the case of stiffened curved panels is investigated and a preliminary design formula is developed.

The interest of this formula for early stages of design is finally illustrated by a short example of cost optimisation.

Design of computer experiments

Background of the design of experiments method

Design of experiment (DOE) methods exist since the beginning of scientific experiment. The first formal theory for the design of experiments in a "modern sense" was published by Fisher [START_REF] Fisher | The design of experiments[END_REF] in the 1920s and 1930s, while working on improving agricultural yield. Since the 1940s, various researchers have promoted and developed the use of experiments strategies in many other areas [START_REF] Kleijnen | Design and analysis of computational experiments: Overview[END_REF]. In the late 1970s, the theory of Taguchi [START_REF] Taguchi | Taguchi's quality engineering handbook[END_REF] on quality improvement made the design of experiment widely used in the industrial environment. In the past 20 years, advances in computational power have led to the study of physical process through computer simulated experiments, which tends to replace physical experiments in cases where the number of variables is too large to consider performing a physical experiment or where it is simply economically prohibitive to run an experiment on the scale required to gather sufficient information.

Computer experiments differ from traditional physical experiment in their deterministic character, meaning that the computer produces identical answers for the same set of experimental parameters. The error in computer experiments is no longer due to random effects which derive from the variability in experimental units, the order of experiments or the locations of the tests. However, it was shown that in many cases, the systematic error between a deterministic model and its approximation has a normal distribution, so that standard statistical techniques can still be applied [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Several authors [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF]Design of numerical experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF]] also insisted on the fact that the selection of parameter's values for computer runs is still an experimental design problem of primary importance, especially considering the quantification of uncertainty of the model on a statistical point of view. Indeed, as not every combination of parameters can be tested, uncertainty and hazard enter the deterministic process through the choice of tested combinations. The design of a computer experiment is hence at the border of a physical and a statistical problem which specificities are emphasized in the following section.

General progress of the design of computer experiment method

Schematically, a numerical model can be considered as a process: the user specifies the combinations of (input) variables to the computer simulator from which the responses (output) are generated. Fig. 2 illustrates this process in the simple case where there are only two input values (X 1 and X 2 ) and one response Y . Each variable can take a value from "low" to "high". The set of all domains of variation forms the "region of interest". In correspondence with each input variable (X i 1 , X i 2 ), the computer program will provide one result Y i . A set of n responses will then generate by extrapolation a response surface. In practice, the explicit formula for this surface is not known. The aim of DOE method is to provide approximated models (response surfaces) that are sufficiently accurate to replace the true response and can be used to facilitate design space exploration, optimisation or reliability analyses.

The general steps of computer experiments are generally similar to those encountered in classical experiments [START_REF] Montgomery | Design and analysis of experiments[END_REF] and can be summarized as follow:

•

Step 1: Statement of the problem.

• Step 2: Choice of the model for the response surface.

• Step 3: Selection of the input data points.

• Step 4: Evaluation of the approximated model.

•

Step 5: Validation of the accuracy of the response model.

• Step 6: Selection of most significant terms and conclusion. This process is not necessarily linear and could be applied iteratively if the predicted model does not meet the desired accuracy. In the above process, the selection of the input variables (step 3) and the technique for approximating the response (step 4) are the two main issues that differ between physical and numerical experiments due to the deterministic property of computer experiments. These two issues will be developed in the following paragraphs. Within the available methods of sampling [START_REF] Franco | Design of numerical experiments in exploratory stage for complex phenomena[END_REF], the following three are the most common and efficient: the Monte Carlo method (MC), the Latin Hypercube Sampling method (LHS) and the Quasi-Monte Carlo methods (QMC) which can be viewed as deterministic versions of MC methods because they use deterministic points rather than random samples. Blatman et al. [START_REF] Blatman | Quasi random numbers in stochastic finite element analysis[END_REF] showed that QMC overperforms MC and LHS, when used with polynomial response surfaces with a mean computational gain factor of 10 in order to reach a given accuracy. The QMC methods are also termed as low discrepancy procedures: sampling points are selected in such a way that the error bound is as small as possible. There are many ways to construct a QMC sequence but the Sobol' sequences are the most widely used because they are quick to construct and fast to converge [START_REF] Krykova | Evaluating of path-dependent securities with low discrepancy methods[END_REF]. They also have the advantage of preserving the uniformity of the distribution when the dimension increases: a Sobol' sequence can be constructed from a shorter one by adding points to the shorter sequence, on the contrary to LHS, where the entire sampling process must be run again.

Response model regression

After selecting the appropriate experimental points and performing the necessary computer runs, the next step is to choose an approximated model and a fitting method. The approximated model must be simple and represent adequately the response of the studied problem. In recent years, a lot of work has been done on approximated models: polynomial response surfaces, neural networks, kriging or multivariate adaptive regression splines [START_REF] Simpson | Metamodels for computerbased engineering design: survey and recommendations[END_REF]. Despite the variety of approximations that is available, comparative studies of these approaches are limited [START_REF] Simpson | Sampling strategies for computer experiments: design and analysis[END_REF]. Depending on the complexity of the problem, one of the aforementioned method might be suitable. However, the polynomial response surface model is by far the simplest; it has been used efficiently in a wide variety of applications and has provided good approximated solutions to even very complex problems [START_REF] Allen | Constructing metamodels for computer experiments[END_REF]. Beside, the use of polynomial response surface for furthers studies such as reliability [START_REF] Zheng | Improved response surface method and its application to stiffened plate reliability analysis[END_REF] and optimisation [START_REF] Roux | Response surface approximations for structural optimization[END_REF] is relatively easy.

Statistic tools for adequacy checking

As mentioned in the section 2, model adequacy checking is an important part of the data analysis procedure. Indeed it is necessary to ensure that the fitted model provides an adequate approximation of the true system and to verify that none of the model assumptions is violated. In most cases, the regression model is a linear function of some unknown coefficients which are identified thanks to the least square method which will be used here for its simplicity and reliability.

Application to cylindrical curved panels under uniform axial compression

In a former study, the authors [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF] investigated the case of unstiffened cylindrical curved panels under axial compression (see Fig. 4) and proposed a set of formulas for the evaluation of the ultimate strength. These formulas were established following the general European Standards procedure for all kind of stability verification and will be used as a reference case to validate the accuracy and relevance of the methodology proposed in previous section.

The strength of the panel χ was hence given as a function of the relative slenderness λ and three parameters λ 0 , β, α Z depending on the relative curvature:

χ = 2β β + λ + β + λ 2 -4β λ -α Z (λ -λ 0 ) (1)
These simulations, as well as those which will be conducted here, involved non-linear material and second-order analyses with imperfection (GMNIA).

They were conducted with the software Ansys version 13 and the standard quadrilateral 4-nodes element [START_REF] Ansys | User's theory manual[END_REF]. Panels were made of elasto-plastic steel with linear hardening as indicated in EN 1993-1.5 C.6.c) (S355, E = 210GP a, ν = 0.3 and a slope of E/100 ). The cylindrical panels were assumed simply supported on all edges and loaded by a uniform longitudinal compression along the curved edges. An initial imperfection with the shape of the first buckling mode and with a maximum amplitude of 1/200 th of the smallest edge was also added. The study was limited to square panels, so that only the thickness, the width and curvature of the panels were varied.

Step 1: Statement of the problem

The aim of this step is to identify in an exhaustive manner the parameters of the problem and to select among them the ones which will have an influence on the response and which are liable of variations in practical applications.

Here the quantity of interest in the panel response (the output) is the ultimate strength of the panel. Basic structural engineering tells us that it is influenced by the geometry of the panels (including their imperfections), their material properties, their boundary conditions and the nature of the loading. All these parameters could be included in the experimental program, but in this first example, the objective is to validate the method and to illustrate its pertinence, so that the same restrictions as in [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF] will be observed:

• the imperfections are chosen following EN 1993-1-5 [START_REF]EN1993-1-5 . Eurocode 3 design of steel structures part 1-5: Plated structural elements[END_REF] (i.e. their shape is that of the first buckling mode and their amplitude is 1/200th of the width of the panel)

• the steel grade is S355 as generally used in modern bridges;

• the panels are simply supported on all edges;

• the longitudinal compression is uniform along the curved edges.

The only varying input factors are thus the dimensions of the panels: their length a, width b, thickness t p and radius of curvature R (see Fig. 4).

Applying the Buckingham-Vaschy's theorem, it can be demonstrated that the ratio of the ultimate strength of the panel and the yield stress (σ ult /f y ) depends on three independent dimensionless parameters:

σ ult f y = f a b ; t p b ; b R (2)
Figure 4: Cylindrical curved panel under uniform axial compression (after [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF]).

Now that the parameters have been defined, it is essential to define the range in which these parameters will vary. Every feasible configuration has to be included but the range of variation has to be kept as small as possible:

it is directly linked with the precision of the approximated expression found at the end of the process. Here the ranges given in table 1 seem reasonable to cover most applications of such panels in bridge engineering.

.6 ≤ a/b ≤ 1.6 X 1 = 2 a/b -2.2 t p /b Slenderness 0.01 ≤ t p /b ≤ 0.04 X 2 = 66.7 t p /b -1.67 b/R Curvature (Angle) 0 ≤ b/R ≤ 1 X 3 = 2 b/R -1
As the order of magnitude of the variations of these three parameters is different, it is preferable to transform the physical parameters into centred variables X i , ranging from -1 (low value) to 1 (high value). Their comparative influence on the response will hence be easier to catch. The three adimensional parameters X 1 , X 2 and X 3 used in the coming paragraphs are thus given in the last column of table 1.

Step 2: Choice of the response surface

The choice of the response surface is based on two issues: the knowledge of the physics of the problem and the desired accuracy of the approximation.

Here the target response is the ultimate strength of the plate, namely the maximum load that the plate can bear when accounting for the elasto-plastic behaviour of the material. From existing standards (EC3), it is known that the strength of a flat plate is related to the slenderness by a second order polynomial which was first proposed by Winter [START_REF] Winter | Strength of thin steel compression flanges[END_REF]:

σ ult f y = 1 λ - 0.22 λ 2 (3)
In ( 3), the slenderness λ is directly related to a/b and t/b (which means to X 1 and X 2 ) by:

λ = f y 12 (1 -ν 2 ) π 2 E • 1 k a/b • b t p (4)
where

k a/b is a function of a/b k a/b =      a b + b a 2 if a b ≤ 1 4 if a b ≥ 1 (5)
It can thus be concluded that a second order polynomial should provide a good approximation of the strength of a curved plate and that it will be meaningful from a physical point of view. The response surface will thus be investigated in the following form:

Y = σ ult f y = β 0 + 3 i=1 β i X i + 3 j=1 j i=1 β ij X i X j ( 6
)
where Y is the approximated response, X i are the three input variables and β (.) are the ten unknown parameters.

Step 3: Selection of the input data or sampling points

The selection of the input data points covers the choice of the number of points and of their distribution in the investigated domain which is here a hypercube in the three dimensional space. The generation of a set of sampling points using a QMC method in this cube can be made easily using a common statistical tool providing a function generating a Sobol' sequence (Matlab here). Noting that with the Sobol' sequence, the extreme values of the parameters can only be reached for an infinite number of variables, some additional points located at the corners of the cube can be added to the sequence to give more weight to the boundary of the domain.

The key issue is thus the definition of the minimal number of experiments to be conducted to get a response with the desired accuracy. The number of simulations n depends on the complexity of the studied phenomenon as well as of the complexity of the approximated model. Yet there is not a unanimously agreed method relating the number of observations versus the number of independent variables in the model. Some authors suggest 3m + 1 points [START_REF] Bowman | Design and analysis of numerical experiments[END_REF] for a second-order polynomial approximation where m = (p+1)(p+2) 2 is the number of unknown coefficients and p is the number of input variables. Following this suggestion for the present example which has 3 input parameters, a second order polynomial approximation will have 10 unknown parameters which could be evaluated with a good accuracy with 31 experiments. Adding the corner points (in total 7 additional points as the point (-1; -1; -1) is by construction always part of the sequence) of the investigated domain, the total number of sampling points is set to 38. An illustration of such a set is shown in Fig. 5.

Step 4: Parameters evaluation

The approximated response surface is here looked for in the form of a second-order polynomial in the three variables X i given by Eq. 6. The unknown parameters β (.) have to be identified from the numerical experiments (here n = 38) which is here done by the least square method. So, from the 38 sampling points shown in Fig. 5, the ultimate strength of a curved plate can be approximated by the following expression:

Y = 0.879 + 0.002X 1 + 0.212X 2 + 0.052X 3 -0.001X 1 X 2 -0.063X 2 X 3 -0.004X 3 X 1 -0.037X 2 1 -0.100X 2 2 -0.003X 2 3 (7)

Step 5: Evaluation of the accuracy of approximated model

To evaluate the accuracy of the approximated model, conducting an analysis of variance (ANOVA) is very useful. The coefficient of determination is first determined R 2 = 0.977 and then the cross-validation coefficient Q 2 = 0.969. The fact that R 2 is very close to 1 indicates that the regression model fits well the data. The small difference between R 2 and Q 2 indicates that most observations have an influence on the regression equation and that the approximation model predicts well the observations. Moreover, the same model has been identified on a sample without additional points in the corners, leading to R 2 = 0.971 and Q 2 = 0.942, the small diminution of these coefficients is a direct consequence of the diminution of the number of sampling points (31 instead of 38), not to the fact that the discarded points were located in the corners of the domain. It is thus concluded that, the addition of corner points is not necessary to get an accurate estimation of the regression coefficients.

3.6.

Step 6: Selection of most significant terms and conclusion.

The second-order formula presented in Eq. 7 for the evaluation of the ultimate strength of curved steel panels under axial compression provides a good and best possible approximation of the real capacity of the panel.

However, it is remarked that not every coefficient in Eq. 7 have the same order of magnitude. So, rather than trying to explain the model with all its terms, it can naturally be asked if some terms could be excluded from the initial model without altering significantly the accuracy of the whole model.

A criterion of exclusion should hence be fixed. If normality assumptions are verified (as in the present case), the t-test provides a fully reliable criterion as it relates the value of each coefficient to its estimated standard error.

More simple criteria, such as arbitrary thresholds of significance are also very effective. Indeed, as the parameters all vary between -1 and 1, the contributions of the various terms can directly be analysed by comparing the coefficients which might then be neglected if their value is bellow a certain absolute value (for example 2% of the sum of the coefficients absolute values or 5% of the maximum value of the coefficients). Fixing here this threshold to 2 % of the sum of the coefficients, the terms (β 1 , β 12 , β 13 and β 33 ) are found not significant (a criterion based on t-test and a 90% two sided interval would give the same results). The new model is thus given by Eq. 8; it preserves good precision with high value of R 2 = 0.976 and Q 2 = 0.9674.

Y = 0.879 + 0.212X 2 + 0.052X 3 -0.063X 2 X 3 -0.037X 2 1 -0.100X 2 2 (8)
Introducing the physical parameters of table 1 into Eq. 8, the ultimate strength of unstiffened cylindrical curved panels under axial compression is given by:

σ ult f y = -0.09 + 0.326 (a/b) -0.148(a/b) 2 + (40.6 + 0.314Z) (t p /b) -(444 + 8.40 Z) (t p /b) 2 (9)
where Z is the curvature parameter defined by Z = b 2 /Rt p . Eq. 9 is very similar to the classical expression of the stability problem, where the ultimate strength is represented as a polynomial function of the slenderness t/b. Fig. 6 shows how well the expressions given by the DOE method (red squares) and by the semi-analytical method [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF] (green triangles) are able to predict the numerical results (given by F.E. model). For most input values, the two models have less than 5 % of discrepancy (in absolute value) to the true numerical value. However a few observations (No. 19, No. 22 and No. 28) predicted by DOE method have higher discrepancy (from 8 % to 10 %) on the contrary to the semi-analytical method whose error remain below 5.5 %. This might be explained by the fact that, although the semianalytical model is not richer (the calibration of λ 0 , β, α Z involves only 7 parameters), its physical bases are finer which proves the crucial importance of the choice of the parameters and response surface. These values of error should be relativised by the fact that in many complex structures such as stiffened plates, the difference between Eurocode predictions and the results of numerical simulations might reach until 20 %, sometimes in favor of safety, sometimes not [START_REF] André | Effect of initial imperfections in numerical simulations of collapse behaviour of stiffened plates under compression[END_REF][START_REF] Ville | Ultimate load of stiffened compressed plates: Effects of some parameters and discussion concerning the ec3 rules[END_REF]. Moreover, as the socalled characteristic value of a member is obtained by dividing its design value (e.g. Eq. 9) by some safety factor (often taken as γ M 1 = 1.1) the present discrepancy is indeed acceptable.

Application to stiffened curved panels under uniform axial compression

The behaviour of stiffened curved panels is a more complex problem, especially due to the interaction of different parameters (curvature, relative rigidity of stiffeners and plate, imperfection, etc.) for which no semi-analytic expression exists. So, as it has just been shown that the design of computer experiment method is well adapted for studying the stability of curved plates, it will be used for the development of a preliminary design formula (i.e suited for hand-calculation) for the ultimate resistance of stiffened curved panels.

Finite element modelling

The stiffened panels are modelled and analysed using the commercial finite element software Ansys [START_REF] Ansys | User's theory manual[END_REF]. The panels are supposed to be simply supported on all edges of the panel (u r = 0 in the cylinder coordinate system of Fig. 7) but not on the stiffeners (unfavourable condition). For loading conditions, the study is here limited to a uniform compression in the longitudinal direction as it is the dominant loading in bottom flange panels. It is applied not only to the main panel, but also to the stiffeners due to their participation in the overall behaviour of the structure (Fig. 8). In fact, in a bridge, the compressive forces acting on the flange come through the diaphragms and webs that connect the upper and lower panels of the box girder. By construction, the stiffeners, in most cases, are continuous and attached by welding to diaphragms: therefore they are also subjected to the compressive load. The curved panels are meshed with eight-nodes shell elements which use an advanced shell formulation that accurately incorporates initial curvature effects (this element is called SHELL-281 in [START_REF] Ansys | User's theory manual[END_REF]). They are well-suited for linear, large rotation and large strain non-linear applications and offers improved accuracy in curved shell structure simulations and faster convergence than plate elements as one can see in figure 9 which represents the convergence study from [START_REF] Braun | Stability of steel plates under combined loading[END_REF]. A fine mesh with more than 30 elements per panel edges is used to reduce the discretisation error. The panels are all made of steel which is assumed to be elasto-plastic with linear strain hardening as indicated in EN 1993-1-5 C.6 for the material non-linear second-order analyses with initial imperfections (GMNIA). The

Young modulus E and Poisson's ratio are taken equal to 210 GP a and 0.3 respectively. The steel grade is S355 with a yield steels equal to 355 MP a.

Evaluation of the ultimate strength

This study is limited to the case of stiffened curved panels under axial compression with open section stiffeners (simple flat plates) because the curvature makes it difficult to realise a close form section of stiffener (boxed rib). Therefore, the number of input parameters is here restricted to seven as presented in table 2. The ranges of variation of these parameters are chosen in order to cover most panels used in bridge construction. As the orders of magnitude of the parameters variations are different, they are transformed the physical parameters into centred variables X i , ranging from -1 to 1.

≤ a ≤ 6 X 1 = a -5 b Width of the panel 4 ≤ b ≤ 6 X 2 = b -5 t p Thickness of panel 0.01 ≤ t p ≤ 0.02 X 3 = 200 • t p -3 1/R Curvature of panel 0 ≤ 1/R ≤ 0.1 X 4 = 20/R -1 d Distance between stiffeners 0.3 ≤ d ≤ 0.8 X 5 = 4 • d -2.2 h s Height of stiffener 0.1 ≤ h s ≤ 0.2 X 6 = 20 • h s -3 t s Thickness of stiffener 0.01 ≤ t s ≤ 0.02 X 7 = 200 • t s -3
The approximated model is searched in the form of a second order polynomial. The total number of experiments, as suggested in section 3.3, is

n = 3 • m + 1 = 109
where m = 36 is the number of unknown coefficients

(1 constant, 7 linear and 28 quadratic terms). Their distribution in the region of interest is generated by a Sobol' sequences and the coefficients are obtained by the least square method, supposing that errors are independent and normally distributed. Like previously, the selection of significant coefficients is made based on the t-test (the limit value being given for 109 tests and the bilateral 5%-95% fractile). Then all remaining coefficients (here 18 coefficients) are re-evaluated using the least square method a second time.

The resulting approximated model (in MN) is thus the following:

Y = +17.09 (10)
-0.47X 1 + 3.58X 2 + 4.24X 3 + 7.32X 4 -3.87X 5 + 4.83X 6 + 2.33X 7 +1.65X 2 X 4 -1.72X 2 X 5 + 1.71X 2 X 6 + 0.89X 2 X 7 + 1.33X 3 X 4 -0.76X 3 X 5 -1.73X 4 X 5 + 0.81X 4 X 6 -1.18X 5 X 6 + 0.94X 6 X 7
The coefficient of determination R 2 = 0.986 (close to 1) indicates that the value predicted by the model fits very well the data (see Fig. 10). Also, the small difference between the values of Q 2 = 0.978 and R 2 shows that there are few undue influence on the regression equation. The accuracy of the response is also checked by the relative mean absolute error (RMAE):

4.2 % (which indicated that a mean error of 4 % is expected), the relative minimum and maximum error: -14.0 % and +13.5 % respectively. The above equation is thus fairly acceptable in bridge constructions and residual errors can be easily covered by using a safety factor. Considering now every coefficients independently, it is found that the 5 % two-sided confidence interval of the regression coefficients is ±0.20 MN for the constant term, ±0.36 MN for the linear terms and ±0.62 MN for the quadratic terms. As generally observed, the uncertainty on linear terms is almost twice smaller than that on quadratic terms: the direct influences of the parameters are better known than those of their interactions.

These statistical remarks being made, it is remarkable that all parameters are found significant in Eq. 10. In decreasing order, the most significant parameters are the curvature (X 4 ), the height of the stiffeners (X 6 ), the thickness of the panel (X 3 ), the distance between stiffeners (X 5 ), the width of the panel (X 2 ), the thickness of the stiffeners (X 7 ) and finally the length of the panel (X 1 ) whose influence is very limited (no more than ±3 % of the total strength). Quite obviously, increasing the curvature, the thickness of the plate or the height and thickness of the stiffeners increases the strength of the panel, while increasing the length of the panel or the distance between stiffeners decreases it. Then the fact that the strength grows with the width is not so immediate but can be easily understood considering that when the width of the panel increases, the distance between the centre of gravity of the panel and the curved plate increases due to curvature and by there the global inertia of the curved panel increases.

There are then multiple interactions which combine effects are more difficult to analyse. Indeed increasing the curvature, the height and thickness of the stiffeners or the thickness of the panel has always a positive effect on the strength because, for these parameters, the linear term dominates clearly the quadratic terms. Then concerning the distance between the stiffeners, in most cases diminishing it leads to an increase of the strength but not mandatory as for slender panels with small curvature and small stiffeners it might lead to a smaller strength (indeed the coefficient of X 5 is -3.87 -1.72X 2 -0.76X 3 -1.73X 4 -1.18X 6 and varies between -9.26 MN and 1.52 MN). In the same way, in most cases increasing the width of the panel leads to an increased strength but not for slender panels when the spac-ing between stiffeners is too large (the coefficient of X 2 being 3.58 + 1.65X 4 -1.72X 5 + 1.71X 6 + 0.89X 7 , it varies between -2.39 MN and 9.55 MN). Very likely in the last two cases, these changes of the coefficient sign correspond to a change in the buckling mode from column to plate or from global to local.

Anyway, it must be recalled that even in these extreme cases, the error in the prediction of the strength is not larger than in other cases (cf. Fig. 10) and that according to all statistical criteria mentioned above, the approximated model given by expression [START_REF] Khedmati | A numerical investigation into the effects of parabolic curvature on the buckling strength and behaviour of stiffened plates under in-plane compression[END_REF]) is able to predict correctly the ultimate strength of stiffened curved panels in the interested domain. It can thus be easily inserted in an optimisation pattern as illustrated in the coming section.

Cost optimisation of curved stiffened panels

The cost optimisation scheme proposed here is based on a cost objective function similar to the ones used by [START_REF] Jarmai | Minimum cost design of a welded orthogonally stiffened cylindrical shell[END_REF]. It assumes that the manufacturing cost of a stiffened curved panel K defined by the parameters X i is the sum of the material costs K m (the steel cost) and of the fabrication costs K f which can be defined as follow:

K(X i ) = K m + K f = k m ρV + k f T i ( 11
)
where ρ is the steel density, V is the total volume of the curved panel, k m and k f are characteristic coefficients of material and fabrication costs. T i denotes manufacturing times:

• T 1 : time for preparing, cutting and assembling the pieces:

T 1 = Θ d κρV (12)
with Θ d a factor characterising the impediment for welding and κ the number of elementary pieces to be welded;

• T 2 : time for welding and T 3 : additional time for maintenance of the machine which might be considered as 0.3T 2 , so that:

T 2 + T 3 = 1.3 C i a 2 wi L wi (13)
where L wi is the length of the i th weld, a wi = max(0.4t s , 4mm) its width and C i a coefficient depending on the welding technique which is here taken equal to 0.2349 for Shielded Metal Arc Welding.

The constraint equation is then given by the stability requirement of the panel:

g(X i ) = N app N ult /γ M 1 -1 0 ( 14
)
where N app is the applied load, N ult the capacity of the panel estimated by Eq. ( 10) and γ M 1 is a safety factor.

The panel which is proposed here for optimisation has fixed overall dimensions: its length a is 6 m, its width b is 4 m and its curvature radius R is 20 m). It is subjected to a uniform axial compression N app = 12 MN.

The objective is thus to determine the parameters (thickness of the panel t p , thickness t s and height h s of the stiffeners and distance between stiffeners d) which will minimize the cost of the panel [START_REF] Park | Estimation of buckling and collapse behaviours of stiffened curved plates under compressive load[END_REF] and verify the constraint equation [START_REF] Fisher | The design of experiments[END_REF]. To make the problem more realistic, it is also considered that the variables are not continuous but discrete (which poses no problem to Matlab optimisation algorithm), so that the solution is looked for in the following domain:

• About the method first, it must be noticed here that the set of optimised parameters corresponding to each ratio k f /k m was obtained almost immediately thanks to the preliminary design formula developed in section 4.2 whereas it would have taken hours or even days using directly finite element simulations. Concerning the results then, following remarks can be drawn:

• The number of stiffeners is higher when only material costs are considered. It is however not maximal (n max = 12) which shows that increasing reasonably the thickness of the plate is very efficient from a weight point of view.

• For higher values of k f /k m , stiffeners becomes logically stiffer to reduce their number and the number of welds.

• It is often more economical to increase the panel thickness than to increase the number of stiffeners which confirms the conclusion of [START_REF] Ville | Ultimate load of stiffened compressed plates: Effects of some parameters and discussion concerning the ec3 rules[END_REF].

Conclusion

Stiffened curved panels in civil engineering structures have high sensitivity to instability phenomenon. Analytical or semi-analytical studies are often not feasible as the problem depends on many parameters such as the panel's curvature or the panel configuration with its stiffener and semi-rigid supports.

There is hence a need for a robust strategy when attempting to develop approximated models for such problems. The proposition of such a strategy was the aim of the first part of the present paper and this, through a turnkey methodology based on the theory of the design of experiment method. The efficiency of the method was first reviewed. Some particular points which differentiate the ordinary physical experiments from computer experiments were discussed. Afterwards this methodology was applied to the case of unstiffened curved panels for which solutions were already available in order to evaluate the accuracy of the method and its relevance. A huge gain of time was noticed when using the DOE method: only 38 simulations were needed in the first application against 524 observations in [START_REF] Tran | Stability of curved panels under uniform axial compression[END_REF] for determining the capacity of a curved panel. Also the general accuracy of the model in the form of a second-order polynomial was comparable to that obtained with more standard heuristic methods. Moreover, as the experiment designer had existing knowledge of the problem, the input values were adequately chosen and the physical interpretation of the results was easy and satisfying, despite the simplicity of the model. The strategy proposed here provides thus a reliable alternative method for the prediction of the ultimate strength of curved panels.

Confident in the methodology, the authors then developed a fully reliable formula for preliminary design of stiffened curved panels. The accuracy of the formula was demonstrated and the influence of various design parameters was discussed. A simple cost optimisation problem was finally presented to illustrate the potential of the formula.

Figure 1 :

 1 Figure 1: Stiffened curved panel of the Confluences Bridge in Angers (France, 2011)

Figure 2 :

 2 Figure 2: Principle of the computer experiment process and response surface.

2. 3 .

 3 Selecting sampling points A good experimental design should minimize the number of runs needed to acquire information with a given level of accuracy. The experimental design techniques were initially developed for physical experiments. Due to the discrepancy associated with physical experimentation, classical DOEs will focus on parameter settings near the perimeter of the region of interest and take multiple data points (replicates) as shown in Fig. 3(Left). Computer experiments are determinists and are not subjected to this necessity. The objective of computer experiments is hence mainly to uniformly distribute the sampling points in the region of interest (such a design is called "space filling") as seen in Fig. 3(Right).

Figure 3 :

 3 Figure 3: (Left)"Classical" and (Right)"Space Filling" designs.

Figure 5 :

 5 Figure 5: Sampling points generated by the Sobol' sequence (circles) and additional corner points (triangles).

Figure 6 :

 6 Figure 6: Comparison of FEM results with DoE and semi-analytical method [12].

Figure 7 :

 7 Figure 7: Boundary conditions of simply supported on all edges

Figure 8 :

 8 Figure 8: Loading condition and scheme of the connection stiffener/diagram by welding

Figure 9 :

 9 Figure 9: Convergence study of SHELL-281 element: linear bifurcation analysis (left), non linear buckling analysis (right)

Figure 10 :

 10 Figure 10: Comparison of the ultimate strength of the FEM and DOE

βa

  Parameter characterising the asymptotic behaviour of the panel β 0 Constant term and average value of the approximated response β i Coefficient characterising the effect of the variable X i β ij Coefficient characterising the interaction of the variables X i and X j χ Reduction factor for the panel buckling according to EC3 λ Relative slenderness of the panel according to EC3 λ 0 Slenderness separating plastic buckling from elasto-plastic buckling ρ Steel density σ ult Ultimate strength of the panel Length Fabrication cost per volume unit K m Material costs (steel cost) k m Material cost per volume unit k a/b Parameter characterising the influence of the aspect ratio m Number of unknown coefficients in the approximated model n Number of simulations or numerical experiments N app Normal force applied to the panel N ult Capacity of the panel p Number of input variables R Curvature radius of the panel T i Manufaturing time of the i th operation t p Thickness of the panel t s thickness of stiffeners V Total volume of the curved panel X i Generic name of the i th input variable X j i j th value of the i th input variable Y , Ŷ Response and approximated response Y j j th value of the response Z Curvature parameter defined by Z = b 2 /Rt p

Table 1 :

 1 Design variables for unstiffened curved panels

	Variable	Description	Variation range	Design variable
	a/b	Aspect ratio	0

Table 2 :

 2 Design variable (dimension in meter)

	Variable Description	Variation range Design variable
	a	Length of the panel	4

  Concerning then the definition of cost coefficients, as no precise data were available for k f and k m , it was decided to present the results in an adimensional form considering different values of the ratio k f /k m . For k f /k m = 0, only material cost is taken into account, while for large values of k f /k m , manufacturing cost prevail (reasonable values in northern countries lie between 1 and 2). The results of the optimisation procedure are shown in table 4.3 (where n is the total number of stiffeners).

	• h s ∈ [0.1; 0.2] by steps of 1 cm;
	• t s ∈ [0.01; 0.02] by steps of 1 mm;
	t p ∈ [0.01; 0.02] by steps of 1 mm;
	• d ∈ [0.3; 0.8] by steps of 5 cm;

Table 3 :

 3 Results of the optimisation procedure for N app = 12 M N

	k f /k m	t p	d	n h s	t s	K
	0.0	0.014 0.55 7 0.16 0.015 3450
	0.5	0.015 0.75 5 0.18 0.016 4250
	1.0	0.016 0.90 4 0.19 0.017 5000
	2.0	0.018 0.90 4 0.17 0.015 6200

Acknowledgement

The authors would like to acknowledge here Prof. B. Sudret from ETH Zürich for the kind advices that he provided during the maturation of this paper.


P.-O Vandanjon

A Coiret

B Muresan Paslaru

A Fargier

R Bosquet

M Dauvergne

A Jullien

D François

F Labarthe

Practical guidelines for Life Cycle Assessment applied to railways project

Current transport systems are not likely to be sustainable, since energy resources may be less accessible, facing a rapidly growing demand and because of environmental constrains. In this context a great interest arose for sharing environmental assessment practices between infrastructures owners and researchers. For this present work a partnership between RFF and Ifsttar was decided to study environmental effects and energy consumption for railways life cycle. It paves the way to a full life cycle analysis (LCA) for environmental burdens of railway infrastructures and aims at identifying the best design practices at the project level related to environmental criteria. This LCA takes into account earthwork, railway structure and the operation phase, i.e., the analysis performed is based on the estimation of energy consumptions due to construction, maintenance and on the traffic. Results of this work could enhance practices for environment preservation for the plan of building 4000 kilometres of new high-speed railways, decided in the French National Transport infrastructure scheme. Applying LCA to the project phase is a rather new approach. It requires defining a global functional unit as a basic hypothesis of any study, different from the previous studies. Then, the lifecycle subsystems are successively investigated (construction, maintenance and use) to determine the best way for project Functional Unit definition and further assessment method application.

1.

INTRODUCTION 1.1 Context and objective

The fact that fossil energy resources are limited and that the demand for these resources is increasing are the basis of the Peak Oil Theory, first introduced as a bell-shaped curve analysis in 1956 by Hubbert [START_REF] Hubbert | Energy resources[END_REF]. The peak oil occurrence should lead to higher energy prices and lower availability. Moreover, rising demand for fossil fuels would continue to drive up energy-related CO 2 emissions, making it all but impossible to achieve the 2°C goal [START_REF] Verbruggen | Views on peak oil and its relation to climate change policy[END_REF]. These projections reinforce the need to model the energy consumption of transport operation phase in the perspective of energy optimization. In France, a workgroup on national environmental stakes -"the Grenelle of the Environment"-has set ambitious goals for transport systems in many areas such as biodiversity, natural resources and climate change. The French National Transport infrastructure scheme was set up in 2011 according to this, and, in particular, it was planned to build 4000 kilometres of high speed railway. In order to optimize these new investments, RFF, the manager of the French railway network, and the French institute of science and technology for transport, development and networks (Ifsttar) launched a common study concerning a complete life cycle analysis (LCA) of railways, including the energy consumption during the operating phase as a parameter for environmental optimization.

The expected results of this optimization are: the reduction of environmental global effects of infrastructures considering geometry and construction; the sparing of natural resources (including energy and water); the reduction of emissions (green house gases, carbon monoxide (CO), lead (Pb), particulate matter (PM),...); the reuse of deconstruction materials in comparison to other possibilities; the development of a method prior to designing the railway line.

Goals and scope of the study are to find out the Functional Unit (FU) best fitting the requirements of such a minimization before starting calculations. It then arose that considering each subsystem of the project (earthworks, rail construction, maintenance and use in relation with the global line geometry) should be investigated in terms of infrastructure key parameters.

LCA of railways infrastructures: specificities and literature review

Since the λ0's, at the international level, Life Cycle Analysis has been considered by road civil engineers and public institutions workers as an efficient way to improve road design. It is a standardized environmental assessment method, which has been developed since 1970 for products manufacturing improvement by USA chemical industries and also for packaging industry in Northern Europe and Switzerland. LCA proved to be well adapted for large numbers of manufactured products. It requires the use of a generic functional unit that accurately defines the product properties and functionalities to consider. LCA proved in road framework [START_REF] Jullien | A comparison of environmental impacts of hot and half-warm mix asphalt[END_REF] to be a method for environmental progress according to results obtained on materials processing and mix design (temperature in hot mix plants, recycling, binder content). Hence examining the part of each process in road construction operation helps to highlight the best practice.

Applied to the present study, the life cycle of railway infrastructure consists of three main phases:

Construction phase including three stages of the life cycle: production of materials (raw materials extraction, manufacturing and processing), materials and machines transport, implementation required for the infrastructure construction (in theory, production of machinery should also be included. In practice, it has been removed after applying cut-off rules); Operation and maintenance phases corresponding to the use of the infrastructure with train traffic and maintenance during life use phase; End of life. At the end of life phase, the infrastructure is demolished. Demolition materials are transported, treated and recycled or stored.

Existing environmental studies on railways are relatively recent [START_REF] Kato | A life cycle assessment for evaluating environmental impacts of inter-regional high-speed mass transit projects[END_REF][START_REF] Facanha | Environmental Assessment of Freight Transportation in the U.S. (11 pp)[END_REF][START_REF] Chester | Life-cycle Environmental Inventory of Passenger Transportation in the United States[END_REF][START_REF] Chester | Environmental assessment of passenger transportation should include infrastructure and supply chains[END_REF][START_REF] Chester | Life-cycle assessment of high-speed railμ the case of California[END_REF][START_REF] Uppenberg | Environmental product declarations (epd:s) for the bothnia line[END_REF][START_REF] Stripple | Life cycle assessment of railways and rail transports -Application in environmental product declarations (EPDs) for the Bothnia Line -IVL Svenska Miljöinstitutet[END_REF] and often aim to compare different modes of transport. Thus, Chester and Horvath [START_REF] Chester | Life-cycle assessment of high-speed railμ the case of California[END_REF] found that for highspeed rail system, energy infrastructure consumption is about 5% of the train consumption according to rate of occupancy. Some authors have rather evaluated environmental impacts of alternative structural materials or alternative solutions of construction to lower consumption of energy, natural resources use and emission of pollutants [START_REF] Aihara | Environmental Evaluation of Railway Material Changes Using LCA[END_REF][START_REF] Kirimura | Railway materials tending to be increasingly environment conscious in Japan[END_REF]. One of the constructive solutions especially analyzed [START_REF] Rozycki | Ecology profile of the german high-speed rail passenger transport system, ICE[END_REF][START_REF] Lee | Comparison of environmental loads with rail track systems using simplified life cycle assessment (LCA)[END_REF][START_REF]Comparing environmental impact of conventional and high speed rail[END_REF] is the ballastless slab track vs. ballasted track. Such ballastless technology does not necessarily increase the total energy consumption of the infrastructure [START_REF] Rozycki | Ecology profile of the german high-speed rail passenger transport system, ICE[END_REF] given its 60 years life and its low maintenance. Nevertheless the NetworkRail study [START_REF]Comparing environmental impact of conventional and high speed rail[END_REF] seems to indicate the opposite, ballastless slab track emitting significantly more emissions than ballasted track. In such case the technological issues are rather studied and concentrate more on structure design Of course this study aims at investigating this issue but it is not the only one. Hence to go further on this topic, and as this study focuses on ecodesign, a functional unit (FU) suitable to this study would be "a double track high speed railway of 400 km with travel time of 1 h 50 mn and with 40 trains per day for each track during 50 years". This functional unit is composed of the following elementary functional units (Fu):

-Fu 1 = about 400 earthwork units composed of cut and fill; an earthwork unit is a 1 km area worked to obtain either cut or fill (subgrade soil, capping layer); -Fu 2 = about 488 civil engineering structures including one tunnel; -Fu 3 = about 400 km of sewer drainage system; -Fu 4 = about 300 km of road rehabilitation; -Fu 5 = about 400 km of railway foundation (sublayer); -Fu 6 = about 400 km of railway equipment (ballast, sleepers, signing equipment, catenaries, substations); -Fu 7 = about 400 km of traffic and maintenance railway. In this paper, preliminary results on Fu1, Fu5 and Fu7 are given.

Outline

The original feature of this present study is then to optimize the design of infrastructures while taking into account earthwork, structure and operation phase within LCA frameworks.

After the introduction, the second section of this article is about the construction phase. This section has two parts. The first one deals with earthwork, a step which has not yet been investigated in the literature through LCA studies. The second one is about the structure (layers above the earthwork). The third section of the article presents the methodology used to assess the operating phase; an experiment on the Rhine-Rhone high-speed railway line is presented. The methodological issues for high speed line assessment are based on a selection of impact categories and will be done as soon as energy contributions to life cycle will be analyzed.

CONSTRUCTION

Concerning public works, the better consideration of environmental issues has led to the progressive incorporation of environmental criteria into public calls for tenders. In this context, the strategy of IFSTTAR was to develop the LCA tool "ECORCE" (a software tool with a specifically developed database from open references, since 2008) [START_REF] Ventura | ECORCEμ a life cycle assessment tool for construction and maintenance of road pavements[END_REF], with respect to both the road code of practice and the environment, and finally to adapt the methodology for railways.

The idea is therefore to define indicators enabling the calculation of parameters to assess the environmental effects; prior to describe the principles to minimise environmental loads. This method shall allow the professions concerned to deduce the best compromise between the different technologies.

Earthwork stage

When considering the potential environmental impacts of the different stages inherent to the railroad construction, the earthwork phase is of crucial importance. As a matter of fact, while the earthwork phase stands for a few percents of total road energy consumption (0.4% of life cycle energy consumption for excavation and 1% for entire construction site energies, [START_REF] Rozycki | Ecology profile of the german high-speed rail passenger transport system, ICE[END_REF]), it may account for a significant fraction of other impact indicators like human chronic toxicity potential (I TP, [START_REF] Huijbregts | Priority assessment of toxic substances in life cycle assessment. Part Iμ Calculation of toxicity potentials for 181 substances with the nested multi-media fate, exposure and effects model USES-LCA[END_REF]). The earthwork phase also affects the ecological state of the adjoining and global environment: i) it suddenly sets a barrier that fragments the habitats, ii) it alters the behaviour of the local micro / macro-fauna through the earthwork yard avoidance and disturbance (e.g. noise, soil, water and air pollutions) and iii) it allows alloctonous species / plants to reach previously inaccessible areas and therefore modify the interactions between authigenic organisms. The magnitude of the induced impacts vastly depends on the geomorphology and thus on the consumed amount of fuel relative to cuts and fills. Based on a previous study devoted to road earthworks for a 8.9 km long area, earthworks units of 1 km were defined as a reasonable unitary entity to deal with for assessment [START_REF] Capony | A earthworks and environmental monitoring-A case study for environmental impacts Modelling[END_REF]. Hence, in order to assess the influence of the geomorphology of soils, the earthwork organization has been monitored for a several-years long period. First, the volumes of cuts and fills were analyzed for every single earthwork items. Then, the fuel consumed by machines was determined as a function of their type, daily work and corresponding earthwork item (see items in Table 1). This was a complex task since the earth engineering machines work cannot be easily isolated at one time, as all the materials excavated and compacted are moved continuously. The obtained figures pointed out that the fuel consumed per cubic meter of moved earth increased linearly as a function of the average crossed distance:

Fuel consumed (L m -3 ) = 3.3 10 -4 Distance (m) + 0.06. (r 2 = 0.97) Finally, the amounts of imported materials (e.g. quicklime, aggregates, etc.) were assessed (Table 1). volume (m 3 ) extrapolated from [START_REF] Capony | A earthworks and environmental monitoring-A case study for environmental impacts Modelling[END_REF] Distances crossed by m 3 of moved earth* (average S.D.) (m) from [START_REF] Capony | A earthworks and environmental monitoring-A case study for environmental impacts Modelling[END_REF] Diesel fuel consumed per m 3 of moved earth (L) [START_REF] Capony | A earthworks and environmental monitoring-A case study for environmental impacts Modelling[END_REF] Equivalent energy to diesel fuel consumed extrapolated for FU 1 (Mj) Based on such obtained relationships, researches are currently led in our laboratory in order to assess the potential environmental impacts and fuel consumptions for various earthwork designs. In details, the previous researches aimed at answering to the following questions. What is the relative contribution of the earthwork phase to the environmental impacts (with respect to these of the railroad structure construction, usage and maintenance)? How should the railroad track be designed in order to minimize the environmental impacts and fuel consumption relative to the earthwork phase? How does the earthwork geometry (e.g. the frequency of cuts and fills, distance crossed by moved earth, slope) and employed construction techniques affect the local environmental and the life cycle assessment for railroad construction?

Railway structure construction Construction scenarios -materials production (FU 5)

The objective of the environmental assessment is to compare different structures. The variable elements must only be those corresponding to changes in structures and their possible consequences. The section covered by the scenario will thus be considered with the same longitudinal section; with the same use so that the design is made for a line carrying the same traffic; on the same construction site; with the same use phase (Fig. 3).

Considered scenarios relate to the construction of high-speed railway infrastructure and have been previously defined in a working session with experts of RFF. Two main types of scenarios of ballasted track are considered and are referred as "Scenario 1" and "Scenario 2." The major difference between these two scenarios lies in the material of the sublayer (Fig. 3): o Graded aggregate (0/31.5) for scenario 1 corresponding to the current conventional railway construction; o Bitumen-bound graded aggregate (0/14) for scenario 2, a material which has been spread on some parts of the French high-speed railway, "LGV East", for a test. This material seems interesting for absorbing vibrations transmitted to the ballast and mechanical performance improvement is expected i.e. ballast life time.

Table 2 -Data for FU5 and capping layer part of FU1 assessmentmaterials energy consumption (EC) and Global-warming potential (GWP) according to the two scenarios.

Indicators EC (MJ) GWP (t.eq.CO2)

Scenario 1 5E+8 6.5E+4 Scenario 2 1.25E+9 6.8E+4

For both scenarios, the ballast layer (31.5/50) will have the same thickness, equal to 35 cm under sleepers. The structure will also be identical and consisting of prestressed concrete monoblock sleepers (1666 sleepers per km); Vignoles 60 E1 rails; fastening systems. The width of the railbed is 15m. The results are given in the Table 2. They are computed with Ecorce from Data given in [START_REF] Jullien | Environmental impacts variability of aggregates production Resources Conservation and Recycling[END_REF]. For the scenarios 1 and 2, the vegetation control scenarios correspond to the "classic" scenario on high-speed railway and are as follows:

-Scenario 1: Three annual passes of weed-killing train are considered. Chemicals used are either foliage-applied herbicides (curative treatments), soil-applied herbicides (preventive treatment) or brush killers. -Scenario 2. This construction scenario may change the "classic" scenario of vegetation control. According to experts from the French National Railway Company, SNCF, the changes adopted are: "a 10 to 14-cm thick bitumen-bound graded aggregate is a naturally resistant support to vegetation development and allows, if applied to the whole platform, to avoid any phyto-sanitary treatment ". This could show the interaction between construction technique and maintenance and has to be assessed in future work

Scenarios of track maintenance (part of FU7)

Two different maintenance scenarios were selected from current practices on the existing network. Thus, we consider a first scenario 1-m1 of maintenance corresponding to a speed level reaching 300 km/h and a second, called 1-m2, for a traffic speed of 320 km/h. According to the background of LGV East, for a traffic speed of 320 km/h, the maintenance operations for ballast tamping apply to 85% of the line each year, instead of 58% for the 300 km/h scenario 1. These additional maintenance operations could cause premature wear of ballast, which would imply a 15% ballast consumption increase.

A first estimation can be examined through simple calculations as follows; Scenario 1-m1 = Impacts of construction x 0.85 at 320 km/h Scenario 1-m2 = Impacts of construction x 0.58 at 300 km/h

Transport scenarios for materials (part of FU 7) For some materials, such as ballast with a large tonnage, different transport scenarios will be considered involving the following modes, Boat; Train; Road. At this time a precise system should be defined to obtain the right order of magnitude of this sub system.

EVALUATION OF ENERGY CONSUMPTION DURING THE OPERATING PHASE

Energy evaluation methodology

To estimate the energy consumption, the train is considered as a point with a mass M. Newton's second law is applied to this point (see the equation ( 1)). This gives the power required by the train. Then the electric consumption is deduced by using a constant ratio which illustrates the efficiency of the traction system, including electric motor and mechanical traction.

kM • = F j -R -M • g • sin(α).
(1) M is the train mass, k the inertia coefficient of rotating masses; the longitudinal acceleration; α the slope, F j the force delivered by the electric motor to the drive wheels. R is the resistance force which is composed of: -Rolling resistance, related to the wheel-rail contact; -Mechanical resistance, consisting of varying viscous friction Fv and dry friction Fs, -Aerodynamic resistance, related to drag coefficient Cx, and weather conditions. This resistance force is detailed in [START_REF] Rochard | A review of methods to measure and calculate train resistances[END_REF] and has, generally, the following form:

R = A + B • V + C • V 2
(2) V is the train velocity; A, B, C: coefficients depending on rolling stock and infrastructure.

Full scale experimental tests for model validation

The trial runs on the new Rhin-Rhone high-speed railway line were carried out during June, July and August 2011. During these tests, geometry, energy and dynamic measurements were performed. Direction and velocity of the wind were recorded too.

The previously presented model (eq (2)) has been validated upon these experimental tests [START_REF] Coiret | Energy consumption induced by operation phase of railways and road infrastructures[END_REF] the average estimated consumption is 79.9MJ per kilometre as the average measured consumption is 75.9MJ per kilometre. This result enables us to compute the assessment of FU7 for the energy consumption (see Table 3). This is consistent, in terms of magnitude with Janic [START_REF] Janic | High-speed rail and air passenger transport: a comparison of the operational environmental performance[END_REF] by example who gives a consumption of 68.4 Mj per kilometre for a TGV (French train) and 79.6 Mj per kilometre for an ICE (German train), given the fact that the actual speed is higher than the speeds given in [START_REF] Janic | High-speed rail and air passenger transport: a comparison of the operational environmental performance[END_REF]. About 80% of consumed energy is used to compensate the resistance force.

Thanks to this model, it is possible:

• To investigate the influence of a chosen geometry path to the energy consumption;

• To calculate the interest to get a rolling stock having energy regeneration brakes;

• To calculate the influence of the speed references to the energy consumption.

Table 3 -Data for FU 7 assessment about 400 km of traffic (40 trains per day for each track during 50 years) considering energy consumption (EC) and Global-warming potential (GWP)

Indicators EC (MJ) GWP (t.eq.CO2) Traffic 1.6E+12 5E+7

CONCLUSION

A possible functional unit was investigated in this paper for high-speed railway ecodesign. Based on former road experience and new train data collected for this study, some orders of magnitude of energy assessment and GWP were determined to validate the framework and some guidelines.

Due to environmental constrains on the existing transport system, several railway projects were launched. For consistency sake, these projects have to be assessed from an environmental point of view. In this paper, we propose a methodology to apply life cycle assessment, to railway projects with a special focus on the construction phase and the operation phase. A method is proposed to quantify the energy spent on the earthwork stage by modeling earth moving. Concerning the structure, two scenarios of construction are presented. The first scenario is classical aggregate subbase, when the second one includes a layer of bitumen bound aggregates. Concerning the operation phase, the electrical consumption due to the traffic on the infrastructure is evaluated. The next step of this work will consist in taking into account the maintenance. The final goal is to propose a tool allowing a global comparison of different projects including different techniques of construction and different geometries.

Figure 3

 3 Figure 3 Synthetic scheme for scenarios selection

Table 1 -

 1 Data for FU 1 assessment Engineered earth volumes (i.e. total cuts volume) and fuel consumptions for moved earth relative to the different earthwork items constitutive of a typical earthwork.

	Earthwork	Engineered earth
	item
	For 400 km

  The moved earth accounts for the part of engineered earth that was not locally processed.

	Revegetation	7.6E+05	50 30	0.08	2E+06
	Treated sub-base layer	9.9E+05	900 260	0.31	1E+07
	Unusable cuts †	4.9E+06	770 190	0.34	7E+07
	Treated fills	3.6E+06	1360 450	0.53	7E+07

*† The unusable cuts are locally stored as lateral earth walls and/or noise barriers

ACKNOWLEDGEMENTS

This work has been achieved with the collaboration of experts from RFF, particularly O. Cazier and A. Grux.