F Schoefs

Risk analysis of structures in presence of stochastic fi elds of deterioration: Flowchart for coupling inspection results and structural reliability

Inspection by non-destructive testing techniques of existing structures is not perfect and it has become a common practice to model their reliability in terms of receiver operating characteristic (ROC) curves. This paper suggests a method of building ROC curves in case of random fi elds of defects on a structure by using polynomial chaos decomposition. Knowledge of spatial distribution of defects allows a reliability analysis to be performed. When selecting stochastic fi nite element analysis to solve this problem, the format is the same as the one chosen for modelling inspections results. The paper shows how to link these quantities (ie. reliability and inspection results) in a risk analysis by using polynomial chaos decomposition as a common language. 1

INTRODUCTION

Reassessment of existing structures generates a need for up-dated materials properties. Commonly, onsite inspections are needed and in some cases visual inspections are not suffi cient for an accurate sizing or detection. For example, non-destructive testing (NDT) tools are required for the inspection of coastal and marine structures where marine growth acts as a mask. In these fi elds, the cost of inspection can be prohibitive and an accurate description of the onsite performance of NDT tools must be provided. Inspection of existing structures by a NDT tool is not perfect and it has become a common practice to model its reliability in terms of probability of detection (PoD), probability of false alarms (PFA) and receiver operating characteristic (ROC) curves. These quantities are generally the main inputs needed by owners of structures in view to achieve Inspection, maintenance and repair plans (IMR) [START_REF] Sheils | Development of a two stage inspection process for the assessment of deteriorating bridge structures[END_REF]. The assessment of PoD and PFA comes either from inter-calibration of NDT tools or from the modelling of the probability density functions of the noise caused by the NDT imperfection and the signal. In this last case, if the noise and the signal depend on the location on the structure then PoD and PFA are spatially dependent. The decomposition on polynomial chaos (PC) is then used to represent the underlying stochastic fi elds. It is then natural to perform reliability analysis in the presence of a stochastic fi eld of defects by using the same format. It is shown that PC decomposition offers in this case a common language for the two inputs of the RBI: the modelling of inspection results and the computing of the structural reliability. This paper presents fi rstly how to defi ne PoD and PFA when damage and detection are stochastic fi elds, ie. spatially dependent (section 2). Section 3 deals with PC decomposition of the marginal distributions of the stochastic fi elds obtained at given locations along the structure. It is shown how to identify the coeffi cients of the PC decomposition from data.

However, it has been shown [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF] that PoD and PFA cannot be directly linked with consequences when an IMR policy has been defi ned. From the decision theory and a Bayesian description of inspection results, it is shown in section 4 how to introduce new events for risk-based inspection (RBI). Expressions of their likelihood as functions of PoD, PFA and γ, the probability of defect presence are presented. Results are then generalised to use ROC curves for RBI. By knowing the spatial stochastic fi eld for the signal and the noise, ROC curves are deduced in each point of the fi eld. This section is illustrated with the ultrasonic measurement of the spatially distributed corrosion of a sheet pile in a marine environment.

Finally, section 5 gives the last output needed for risk assessment based on reliability analysis: the computing of the probability of failure. From loading modelling and by knowing the stochastic fi eld of defects, a stochastic computation can be performed in view to estimate the spatial distribution of probability of failure. This section ends with a complete application of risk analysis of a corroded sheet-pile in steel.

ASSESSMENT OF STOCHASTIC FIELDS OF DEFECTS THROUGH NDT TOOLS

Hazards and uncertainties when assessing defect sizes on structures

Let us consider the assessment of a defect d through a NDT tool. Inspection is generally not perfect. We consider the case where harsh conditions affect the measurement of the d. Then, depending on the environment and the location of the inspected area, the measurement of the size is affected by hazards.

We model this quantity as a stochastic fi eld d(X, t, θ), where X denotes the spatial coordinates, t the time and θ the hazard. When assessing this defect, the chain of measurement (accessibility of the inspected area, calibration of NDT tool, experience and tiredness of the operator) introduces uncertainties, called noise in the signal processing theory [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. The result of the inspection, ie. the measured defect size, is a stochastic fi eld denoted ˆ( , , ) d X t . In section 2.2, d(X, t, θ) and ˆ( , , ) d X t are written as d and d for simplicity.

Modelling capability of NDT tools in a probabilistic context

The most common concept that characterises inspection tool performance is the PoD. Let a d be the minimal defect size, under which it is assumed that no detection is done. Parameter a d is called the detection threshold in the following. Thus, the PoD is defi ned as:

PoD d P d a (1)
The detection threshold a d is a deterministic parameter or a random variable. This defi nition implies that PoD is a monotonic increasing function.

Detection theory allows for the introduction of another useful quantity, the PFA. From the signal theory, d is called the complete signal or "signal + noise" with a probability density function (pdf) f SN . The real defect size d is unknown and it is deduced from the knowledge of the noise η by the relationship d d . The pdf of noise is denoted f N . Noise depends on environmental conditions, human interference and the nature of what is being measured. Let's assume that noise and signal amplitudes are independent random variables, then [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF].

PoD and PFA have the following expressions (2) and

(3):

ˆPoD d SN a f d d (2) PFA ( ) d N a f (3)
where δ denotes the increment. Figure 1 illustrates the corresponding pdf in the case where variables are normally distributed and the computation of PFA and PoD for a given detection threshold.

For a given detection threshold, the couple (PoD, PFA) allows the performance of a NDT to be defi ned; this is the ROC. This couple can be considered as the coordinates of a point. Let us consider that a d is unknown and takes values in the range [-∞; +∞], this point belongs to a curve called ROC curve. It is a parametric curve defi ned by equations ( 2) and ( 3) with parameter a d .

The ROC curve, which is plotted in figure 2, is computed using the pdf presented on figure 1. From a theoretical point of view, this is a convex curve corresponding to a monotonically increasing function, always lying above the diagonal line, where the fi rst derivative is closely linked to the sensitivity of the receiver [START_REF] Arques | Decisions en traitement du signal[END_REF][START_REF] Fücsök | Human factors: the NDE reliability of routine radiographic film evaluation[END_REF].

The diagonal line running from the lower left to the upper right (curve "PoD=PFA") is the line of no "performance", since in that case the inspection PoD equals PFA whatever the detection threshold [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF].

inspections is reasonable and allows a rich data base of information to be accessed.

We follow the second approach in this paper. Note that the PFA is named PFI (probability of false indication) too. The defi nition of discrete values for PFI is generally expressed as a percentage of false indications on the inspected length [START_REF] Silk | Estimates of the probability of detection of fl aws in TOFD data with varying levels of noise[END_REF][START_REF] Rudlin | Inspection reliability and periodicity for rail axle inspection[END_REF]. Finally, as the number of samples is limited, authors provide confi dence bounds: 90% PoD for example [START_REF] Barnouin | Underwater inspection reliability trials for offshore structures[END_REF][START_REF] Rudlin | Inspection reliability and periodicity for rail axle inspection[END_REF].

In some cases, the performance of NDT tools depends on the location of the point to be inspected on the structure. For example, when locations of defects are on welded joints in the sea, the corresponding PoD and PFA should be changed according to the access, the luminosity and the wave shaking for instance. When defects are continuous fi elds on the structure, the PoD and PFA should be indexed by the coordinates X of the inspected point M. We consider here that the defect is induced by a deterioration mechanism indexed by time t and can be modelled by a stochastic fi eld: d(X, t, θ).

Defi nitions of PoD and PFA for stochastic deterioration model

As described in the previous section, after inspection with a NDT tool, the size of the defect d(X, t, θ) becomes ˆ( , , ) d X t . Then the noise η(X, t, θ) is defi ned from the knowledge of these two stochastic fi elds by equation (4).

, , , , , , X t d X t d X t (4)
For each location X and at each time t, the PoD and PFA are, respectively, computed from equations ( 2) and (3). Thus the PoD and PFA, as well as the ROC curves, are functions indexed by X and t. We denote PoD(X, t) and PFA(X, t) as functions of time and space. The fi eld d at a given time t being assessed from inspection, the computation of PoD(X, t) and PFA(X, t) requires the knowledge of one of the other stochastic fi elds in equation ( 4): d(X, t, θ) or η(X, t, θ). Two situations can be considered: (i) The noise is known because it is constant whatever the location of the NDT tool on the structure or because it is constant on given areas on the structure. It is generally time invariant and zero mean. (ii) The real size is known because it has been measured before on-site inspections as in the ICON project [START_REF] Barnouin | Underwater inspection reliability trials for offshore structures[END_REF] or because an assumption is made.

In both situations, the defi nition of continuous spatial functions needs the complete characterisation of the stochastic fi elds by their joint distribution. Practically, no models of joint distribution are available for onsite inspection and almost all NDT tools give data on Looking for the best detection performance, the PoD should always take larger values than the PFA (low noise sensitivity). We have then PoD ≥ PFA. When analysing ROC curves, one must keep in mind that the PFA depends on the noise and detection threshold only. It does not depend on defect unless the noise depends on defect size. That is the case for instance if the operator adjusts the device to detect smaller defects when the current adjustment does not give any signal. The PoD is a function of the detection threshold, the defect size and the noise. Thus, for a given detection threshold, the PFA is a constant, but the PoD is an increasing function of the defect size. The ROC curve is a fundamental characteristic of the NDT tool performance for a given defect size. A perfect tool is represented by a ROC curve reduced to a single point whose coordinates are (PoD, PFA) = (1, 0). The distance between this "best performance point" and the ROC curve is a measure of the NDT ability [START_REF] Schoefs | Multiple inspection modeling for decision making and management of jacket off-shore platforms: effect of false alarms[END_REF]; we call it NDT-PI (NDT Performance Index). Theoretical ROC curves are presented in [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF], where each one is obtained for several signal/ noise ratios of a NDT tool.

When assessing a ROC curve, two approaches can be considered:

• the statistical approach, which aims to provide discrete values for PFA and PoD and requires the knowledge of the real size; this is performed within specifi c inter-calibration projects such as the ICON project in the fi eld of offshore platform analysis [START_REF] Barnouin | Underwater inspection reliability trials for offshore structures[END_REF], when on-site measurements on real structures are costly or when on-site conditions are affected by a lot of factors as for the detection of corrosion pitting [START_REF] Pakrashi | An Image Analysis Based Damage Classification Methodology[END_REF]) • the probabilistic approach, which requires a modelling of the (signal + noise) and noise pdf; it is generally preferred when the cost of on-site specifi c locations; thus, only marginal distributions can be assessed. Moreover, the owner generally does not base its inspection planning on RBI methods; then the inspection campaign aims to give a global overview of the state of the structure. Thus the distance between measurements is generally larger than the distance of correlation and additional assumptions on the structure of correlation for the stochastic fi elds are needed.

Finally, note that the knowledge of ageing laws for d allows predicting evolution of ROC curves with time. This can be implemented in a RBI method for optimising the period between inspection, as well as the location of inspections. For example, for corrosion processes, several models are available (Melchers, 2003;Paik et al, 2003a;Paik et al, 2003b;[START_REF] Soares | Reliability of maintained corrosion protected plates subjected to Paik[END_REF]. In the fi eld of uniform corrosion, spatially dependent ROC curves are available (Schoefs et al, 2007a).

MODELLING RANDOM FIELDS OF DEFECTS FROM NDT MEASUREMENTS

Only the spatial dependence is addressed here. We suggest a representation with PC decomposition for modelling marginal distributions or random fi elds.

It allows for the systematic identifi cation of random variables or random fi elds. We choose the estimate of maximum likelihood for the identifi cation of PC decomposition [START_REF] Desceliers | Identifi cation of chaos representations of elastic properties of random media using experimental vibration tests[END_REF]. This method has already been applied for the identification of random variables from structural monitoring (Schoefs et al, 2007b). The question is to identify the coeffi cients d i , ˆi d and η i of the one-dimensional PC decomposition for every random variable. For generality, let us denote x(θ) as the random variable and x i as the PC coeffi cients. Equation ( 5) gives the general form of a one-dimensional PC decomposition.

( ) ( )

p i i i x x xh ( 5
)
where p is the order of the PC decomposition; ξ(θ) is the Gaussian germ, ie. a standardised normal variable; and h i is the Hermite polynomial of degree i. By using the maximum likelihood method, coeffi cients x i are solutions of the optimisation problem:

( ) argmax ( ) L L (6)
where is the vector of components x i ( = [x 0 , …, x p ]) with dimension (p+1), and L is the likelihood function:

1 ( ) ( ); N x j j L p x (7)
where x(θ j ) is the j th realisation (here measurement) of x; and p x (.; ) is the probability distribution of x, parameterised by .

The likelihood function ( 7) takes very fair values close to the numerical precision. Then the problem ( 6) is replaced by equation ( 8):

Log ( ) arg min Log ( ) L L (8)
The algorithm for solving the optimisation problem (8) is detailed in Schoefs et al (submitted) and only the main steps are described here. Basically, the question is to build a solution by stating conditions (9):

0 0 2 1 Var( ( ( )) p i x i x x x (9)
where μ x and σ x are, respectively, the statistical average and standard deviation of variable x(θ).

They are computed from the N measurements x(θ j ).

In equation ( 9), the fi rst condition reduces the number of unknown PC coeffi cients to p and the second one facilitates the search for other coeffi cients on an hyper-sphere with radius σ x . Moreover, by denoting x i * the quantity x i /σ x , the second condition of (9) becomes equation ( 10):

2 * 1 1 p i i x ( 10
)
This new condition facilitates the search for p coeffi cients x 1 * , ... x i * , ... x p * on a hyper-sphere with radius 1. This last condition is interesting for the optimisation process. Let us now consider a PC of order 3. To parametrise the hyper-sphere, we can introduce two angular parameters φ 1 and φ 2 : A two-step optimisation fl owchart is used solving equation ( 8) by knowing equation ( 11):

*
• a first localisation of the minimum is found through Monte Carlo simulations (size 100 and 1000, respectively, for PC of order 2 and 3) • starting from this point, the Nelder-Mead Simplex Method is used (Lagarias et al, 1998).

This process allows us to avoid a pseudo-convergence around minimum values. Figure 3(a) presents the fi tting of distributions in the case of uniform corrosion of piles in harbours for the two variables: measured loss of thickness d and noise η with PC for several orders p. Figure 3(b) plots iso-values of -Log(L) for PC chaos of order 3, in the space (φ 1 , φ 2 ). ROC curves are then deduced from integration of these decompositions. Thus both PoD and PFA are functions of polynomial coeffi cients x i . Denoting ˆ and the vectors, respectively, of i th components ˆi d and η i , expressions of PoD and PFA as functions of PC coeffi cients are presented in equations ( 12) and ( 13).

When d and η are independent, a germ of dimension 2 is needed respectively (ξ = (ξ 1 , ξ 2 ). Note that the PC decomposition of signal d is obtained from equation ( 4) by knowing ˆ , , ξ 1 and ξ 2 . where p ξ is the probability density function of ξ (standard normal pdf). Practically, integrals ℑ i (a d ) in equations ( 12) and ( 13) are computed through Monte-Carlo simulations using 10 6 samples. These quantities are independent of the study: they can be preprocessed once for all and used for each application.

The stochastic fi eld of corrosion for sheet piles in harbours is spatially indexed by the vertical abscissa only (Schoefs et al, 2007a). ROC curves at six depths between -1 m and +2 m along a pile are plotted in fi gure 4 for a PC chaos of order 3. Abscissa (+1 m) and (+2 m) represent points located in the tidal area where others refer to locations in the immersion area.

Note that the identifi cation of a variance-covariance matrix is generally more complex due to the distance of inspected points that is much larger than the correlation length as mentioned previously in section 2.3. In this case we suggest using independent germs ξ i (θ) for each inspected point and to interpolate PC coeffi cients between two inspected points (14).

, ( ) .

i k i j j d X X d X (14)
where φ k (X) are linear functions for interpolation from knowledge of coefficients at coordinate Corroded sheet-pile and profi le of the loss of thickness.

X j . Note that this representation needs at least measurements where the average trajectory of the process takes its minimum and maximum value.

In the case of corrosion of coastal structures, shape of profiles and position of extreme values are generally known and inspections are carried out at these positions (see fi gure 5). In view to compare this method with another one based on the fi tting with a classical pdf, we used the NDT-PI: it is the distance between the best performance point in (PFA, PoD) space of coordinates (0, 1) and the ROC curve.

Table 1 gives a comparison between the NDT-PI obtained with predefi ned pdf (table 1(a PC coeffi cients (order 3) identifi ed at several depths.

We notice that the normal pdf and order 1 PC decomposition are the same, leading to the same result. [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF] further developed this methodology by focusing on the probability that a defect exists after an inspection has been carried out.

RISK-BASED DECISION AID-TOOLS BASED ON INSPECTION RESULTS

Rouhan

The PoD is the probability that an existing defect is detected, whereas the decision is based on the event that a defect exists, given the results of an inspection. We consider the events:

• the real state δ: "δ = 1" means that "defect is present" • "D(δ)" the decision on this state after inspection: "D(δ) = 1" means that "defect is detected".

From Bayesian point of view, PoD can then be defi ned as the conditional probability of the event "(D(δ) = 1|δ = 1)". The conditional probability of the event of interest "(δ = 1|D(δ) = 1)" can be evaluated using Bayes Theorem. It is subsequently introduced into cost functions, which are used to investigate the effect of cost overrun due to inaccurate inspection results. Four events E i (i =1:4) are then defi ned. For example, the expression of the probability of event E 3 "(δ = 1|D(δ) = 0)" is given in equation ( 15). This event is of fi rst importance in risk analysis because it can lead to a failure.

3 1 PoD 1 PoD 1 PFA 1 P E (15)
Other events are:

• E 1 , is the event "absence of defect knowing no detection" or "(δ = 0|D(δ) = 0)" • E 2 , is the event "absence of defect knowing a detection" or "(δ = 0|D(δ) = 1)" • E 4 , is the event "presence of defect knowing a detection" or "(δ = 1|D(δ) = 1)".

They can be associated to costs and subsequently introduced in a complete RBI planning [START_REF] Sheils | Development of a two stage inspection process for the assessment of deteriorating bridge structures[END_REF]. The occurrence of these events also depends on another parameter, γ, the probability of presence of a defect. The probability density function of γ is related to the natural size of existing defects and their spatial variation (which comes from expert judgment or an historical data base).

An alternative approach is to introduce a parameter known as the probability of indication (PoI):

PoI(d) = PoD(d) + (1 -PoD(d))PFI ( 16
)
where PFI is the probability of false indication [START_REF] Straub | Modeling dependency in inspection performance[END_REF]. The parameter PFI must be defi ned relative to a specifi c size of the inspected area [START_REF] Straub | Modeling dependency in inspection performance[END_REF].

From equations ( 12) and (13), P(E i ) and PoI are expressed as functions of spatial dependant quantities d and . Figure 7 gives the evolution of P(E 2 ) and P(E 3 ) for three assumptions of expert judgement and the result after selection of a standard expert judgement. This factor strongly affects the shape of these curves.

RELIABILITY ANALYSIS BASED ON STOCHASTIC FINITE ELEMENT IN CASE OF RANDOM FIELDS OF DEFECTS

The stochastic (or probabilistic) finite element method refers to finite elements methods that account for uncertainties in the geometry or material properties of a structure, as well as the applied loads. Since the 1980s, quite a lot of developments have been made. They are mainly concerned with uncertainties of materials properties or hazards in loading. The question of the random geometry due to defects during building or the effect of ageing laws (corrosion) is still a challenge. Methods such as X-SFEM are effi cient in this context (Nouy et al, 2007;2008a;2008b;[START_REF] Clément | Reliability of corroded historical lift-cranes in coastal area[END_REF]. In this paper, the computation is based on the so-called non-intrusive stochastic fi nite element analysis [START_REF] Berveiller | Stochastic finite element methods: intrusive and non-intrusive approaches for reliability analysis[END_REF]. We consider in this section that the loss of thickness d(X, t, θ) can be introduced in reliability through "macro" variables, ie. quantities such as diameter, section or inertia are explicit function of d. Then no specifi c method that deals with random geometries such as X-SFEM is needed.

Firstly, a quantity of interest or a limit state is selected.

Here, we are interested in displacement of each point of the structure displacements u(X). The limit state is then expressed: G(X) = u(X)u c , where u c is the acceptable critical displacement. The structural computation is made by using deterministic software.

The question is then to fi nd an approximation ( , ) u X of the PC decomposition u(X, ξ). By knowing the decomposition of basic input variables such as d(X, t, θ), on a PC decomposition with germ ξ, ( , ) u X is expressed in equation ( 17) and can be assessed by computing u α through the quadrature formula (18).

( ) ( ) u u u H (17) ( ) ( ) k k k u H u dP E H u H u ( 18
)
where α is a multi-index, P is a measure of probability associated with ξ, and H α is an Hermite polynomial associated to this multi-index. θ k are the elementary events where the deterministic computation is needed (generally called point of integration or Gauss points in the case of normal measure of probability) and ω k is the corresponding weights associated with the measure of probability P.

A post-processing allows the computation of the probability of failure P f to be carried out with very low computational costs. This is of primary interest when carrying out risk analysis where a good assessment of P f is needed.

For illustration, we can consider an application where an explicit form of the solution can be obtained by the static balance of the structure assuming an elastic behaviour of the material (see fi gure 5).

Considering here a limit state relying on yield momentum, the probability of failure with depth is plotted on fi gure 8, with the loss of matter on fi gure 6 and the confi guration of loading on fi gure 5. Evolution of P(E 2 ) and P(E 3 ) with depth (top) for several expert judgements γ (in bold line 0.9, in fi ne line 0.5 and in dotted line 0.1); standard expert judgement (middle); and corresponding evolution of P(E 2 ) and P(E 3 ) with depth with this standard judgement (bottom).

CONCLUSION: FLOWCHART FOR COUPLING RESULTS OF INSPECTION AND STRUCTURAL RELIABILITY

Sections 4 and 5 provide the main steps for a fl owchart that couples the results of inspection and structural reliability with a formulation based on PC decomposition. From inspection results and noise modelling, it is based on the PC decomposition of basic variables describing the risk based modelling of inspection results ˆ( , , ) d X t and η(X, t, θ) on the one hand, and a PC decomposition of geometrical basic variables for the reliability study d(X, t, θ) on the other hand. The corresponding terms of risk analysis P(E 3 ) and P f can be assessed. The fl owchart in fi gure 9 resumes the steps of methodology. When considering only event E3 (there is a defect knowing we do not detect it), which leads to failure and the probability of failure, fi gure 8 presents the result of the integrated risk analysis: the evolution of risk with depth. The PC decomposition is shown to be an interesting common language between reliability analysis and RBI when stochastic fi elds of degradation occur. Note that quite a lot of uncertainties and hazards are introduced:

(i) The natural hazard of the corrosion process: d(X, t, θ)

(ii) The uncertainty in the inspection results from on-site measurements: η(X, t, θ)

(iii) The uncertainty of the previous knowledge on the corrosion process: γ (X, t). Evolution of probability of failure and selected P(E 3 ) with depth (top and middle), and result in terms of risk (bottom).

On site inspections results

Noise

Figure 1 :

 1 Figure 1:Illustration of PoD and PFA (signal and noise normally distributed)[START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF].

Figure 2 :

 2 Figure 2:Receiver operating characteristic curve.

Figure 3

 3 Figure 3(a): Fitting of loss of thickness and noise with polynomial chaos for several orders p (corrosion of piles in harbours).

Figure 3

 3 Figure 3(b): Representation of -Log(L) in the space (φ 1 , φ 2 ) for PC chaos of order 3.

Figure 4 :

 4 Figure 4:Comparison between ROC curves coming from polynomial chaos with order 3 identifi cation (left) and experimental data (right) for several inspections in depth.

  Figure 5:

  )) (Normal, Generalised Extreme Value (GEV), Student pdf) and PC identifi cations (table 1(b)) at a given position along the vertical axis (X = +2 m). The value of NDT-PI deduced from the original curve (experimental) being 0.054, only PC chaos of order 3 allows us to reach this value (0.052) when both lower order PC decomposition and predefi ned pdf leads to values of around 0.075. The tables give the detection threshold a d at the best performance point. Values of coeffi cients are given at each inspected level on fi gure 6 for PC decomposition of order 3.

  Figure 7:

  Figure 8:

Table 1 (

 1

	ROC curves	X = +2 m NDT-PI	a d
	Experimental	0.054	0.25
	PC order 1	0.079	0.26
	PC order 2	0.072	0.026
	PC order 3	0.052	0.026

b): Values of NDT-PI and corresponding detection threshold for PC decomposition (order 1 to 3).

Table 1 (a):

 1 Values of NDT-PI and corresponding detection threshold for predefi ned pdf.

	ROC curves	X = +2 m NDT-PI	a d
	Experimental	0.054	0.25
	Normal	0.079	0.26
	GEV	0.081	0.29
	Student	0.072	0.25

  Flowchart for coupling results of inspection and structural reliability. corrosion and compressive loads", Marine Structures, Vol. pp. 425-445. Lagarias, Reeds, J. A., Wright, M. H. & Wright, P. E. 1998, "Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions", SIAM Journal of Optimization, Vol. 9, pp. 112-147. Melchers, E. 2003, new model for marine immersion corrosion in structural reliability assessment", International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP'03). Nouy, A., Schoefs, F. & Moës, N. 2007, "X-SFEM: a computational technique based on X-FEM to deal with random geometry", Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, special issue on "New trends in the development of the X-FEM", Vol. 16, No. 2, pp. 277-293. Nouy, A., Clément, A. & Schoefs, F. 2008a, "An eXtended Stochactic Finite Element Method to deal with geometrical uncertainties in structural analysis", 8 th World Congress on Computational Mechanics (WCCM VIII) and 5 th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2008), Venice, Italy, 30 June to 5 July. Nouy, A., Clément, A., Schoefs, F. & Moës, N. 2008b, "eXtended Stochastic Finite Element methods for solving partial differential equations on random domains", Computer Methods in Applied Mechanics and Engineering (CMAME), to be published.

	modelling
	Uncertainty on measurement
	Random fields
	'signal + noise',
	'noise'
	PC decomposition	Random fields 'signal'	Natural hazard
	Inspection based decision aid-tools	Structural reliability -SFE
	Spatial dependent
	risk analysis
	9:

ACKNOWLEDGEMENTS

The author is grateful for the support received from European Community within the Interreg IIIB project MEDACHS (2005MEDACHS ( -2007) ) and for the contribution of Alexandre Clement (Nantes Atlantic University).


Hussam El Cheikh

Institut de Recherche en Gén ie Civ il et Mécanique (GeM ) UMR CNRS 6183, LUNAM Un iversité de Nantes, Ecole Centrale de Nantes, France *1 Hussam.elcheikh@univ-nantes.fr; 2 Bruno.courant@univ-nantes.fr Abstract--Direct Laser Fabrication (DLF) is a rapid manufacturing process based on laser cladding technique in which parts are fabricated layer by layer. In this work, a coaxial nozzle is used to deposit 316L stainless steel on a low carbon steel substrate. In order to obtain a 3D simulation with reasonable computation time, the thermal and geometrical problems are decoupled. The clad geometry is supposed to be a half cylinder which corresponds to a suitable geometry to build a part layer by layer. The finite element method is used to simulate the heat conduction in the part during the process with a 3D modelling using the Arbitrary Langrangian-Eulerian (ALE) in (Multiphysics Comsol 3.5). This model allows predicting the temperature field in the part during the process.

Keywords--Laser; Laser Cladding; Direct Laser Fabrication; Moving Heat Source; Finite Element Analysis; Thermal Field

I. INT RODUCTION

Direct Laser Fabrication (DLF) is a manufacturing method issue from rapid prototyping and laser technologies. Many works have been proposed to simu late the thermal field induced in a material heated by a moving laser beam during the 80's and 90's years. In the 2000's years many authors used these models mo re or less adapted to simu late the DLF process.

Rosenthal [START_REF] Cheikh | Mathematical theory of moving sources of heat and its application to metal t reat ments[END_REF] proposed a 3D analytical function to describe the shape of the melt pool in the case of a moving heat source above a semi-infin ite solid. El Cheikh et al. [START_REF] Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF] have shown that the clad geometry is not controlled by the powder distribution into the jet. They propose [START_REF] Cheikh | Prediction and analytical description of the single laser track geo metry in direct laser fabricat ion fro m p rocess parameters and energy balance reasoning[END_REF] a phenomenological circular form model taking into account the surface tension forces between the protecting gases, the laser clad material and the substrate material. Oliveira and al. [START_REF] Oliveira | Analysis of coaxial laser cladding processing conditions[END_REF] present a theoretical and experimental study of the coaxial laser cladding process to understand relations between the process parameters and geometrical characteristics of one laser track. Picasso and Hoadley [START_REF] Picasso | Finite element simu lation of laser surface treat ments including convection in the melt pool[END_REF] proposed a two dimensional fin ite element model to calculate the melt pool shape for a known clad height. He and Mazumder [START_REF] He | Transport phenomena during direct metal deposition[END_REF] assumed a Gaussian distribution during the interaction between the powder feed rate and the laser beam, and Kar and Mazu mder [START_REF] Mazu Mder | Solid solubility in laser cladding[END_REF] proposed a one dimensional study for the dissolution of powder and mixing of the clad. Toyserkani and al. [START_REF] Toyserkani | 3-D fin ite element modeling of laser cladding by powder injection: effects of laser pulse shaping on the process[END_REF], [START_REF] Toyserkani | 3-Dfin ite element modeling of laser cladding by powder deposition: Effects of powder feed rate and travel speed on the process[END_REF] use a transient threedimensional finite element modelling to show the influence of many parameters on the laser cladding process. Bamberger et al. [START_REF] Bamberger | Calculat ion of process parameters for laser Alloying and cladding[END_REF] use a simplified theoretical model of Mie theory to show the influence of the injected part icles on the heat flu x and temperature distribution. Cline and Anthony [START_REF] Cline | Heat treating and melting material with a scanning laser or electron beam[END_REF] found the temperature distribution in the steady-state case using the Green's functions. In this paper, a thermal simu lation, using the Arbitrary Langrangian-Eulerian (A LE) mode, is used to access to the temperature distribution in the work piece and to simu late the clad layer geometry. The clad geometry is supposed to be a half cylinder which corresponds to a suitable geometry to build a part layer by layer. The thermal simulat ion is obtained by the calculation of the conduction equation II. EXPERIMENTAL PROCEDURES An easy way to comply with the journal paper formatting requirements is to use this document as a template and simp ly type your text into it.

A 5 axes speed machining is used associated with a fibre laser and a coaxial no zzle specially designed according to the laser beam characteristics. Argon, with a flo w rate of 3l/ min, carries the powder while a secondary gas, Argon with a flow rate of 5l/ min, is shaping the powder stream. This nozzle is also equipped with a cooling channel to dissipate heat.

The laser beam is focused on the substrate surface while the powder stream focused plan is 5 mm above. The diameter of the focused powder jet is 0.6 mm and the focused laser diameter is 0.53mm at the substrate surface.

A 316L powder, with particles size between 45 and 75 µm, is used to realize single clad tracks. These experiments are realized with various parameters values. Three different powers (P =180, 280 and 360 W), three different powder feed rates (Qm= 0.025, 0.05 and 0.075 g/s) and three different velocities (V= 300, 600 and 900 mm/ min) are tested.

Figure 2 shows the conic powder jet stream geometry under the nozzle. Horizontal and vert ical image analyses are studied. On the upper part of the jet, the horizontal image analysis displays two distributions due to the coaxial form of the nozzle. Fro m the powder focus plan one distribution only remains. Image analysis software is used to study the powder distribution in the jet counting on the light reflection. The obtained functions are not directly the powder distribution in the jet but rather a pro jection of this distribution. Nevertheless, it can be deduced that, in the focus plan, the powder distribution is an almost symmetrical Gaussian type distribution. The vertical line displays the distribution along the z d irection. As shown in the image analysis, the higher concentration of the powder stream is in the focus plan. The coaxial no zzle allows a perfect symmetry of the part icles distribution in the stream and very few part icles outside the main t rajectories. The powder focus plan, at a distance l (l=3.5 mm) fro m the nozzle, is a d isc with a diameter dp (dp=0.6 mm). The distance work L is fixed to 5 mm. Under the powder focus plan the powder distribution remains a Gaussian type distribution. The distance between the focus plan and the substrate surface is g iven by L-l (1.5 mm) and in the interaction plan the powder stream is 2 mm in d iameter. The laser beam is focused on the substrate surface wh ile the focused plan for the powder stream is 5 mm above.

III. PREDICTION OF THE CLAD GEOMET RY

El Cheikh et al. [START_REF] Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF] studied the possibility to analytically relate the clad geo metry and the powder distribution in the jet stream. They have established these relationships in case of three different distributions and have shown that the uniform distribution leads to semi cylindrical geometry. They have also shown that the powder accumu lation induced by the powder distribution can't exp lain the clad geometry. The clad geo metry is controlled by the tension forces which lead in all cases to a partial cylinder geometry. However the semi circular geo metry is a particular case suitable to build a part layer by layer. In order to obtain an analytical description of the clad geometry, the analytical formulat ion established in [START_REF] Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF] is used. The clad height h, in the end zone (interaction between the powder jet and the laser), is proportional to the powder feed rate, and in inverse proportion to the laser speed. Where P e is the powder efficiency, ρ is the material density and r p is the radius of the interaction zone between the powder jet and the substrate.

H being the clad height after the laser nozzle past (Zone 2 in the figure 3), it is given by: IV. THERMAL MODEL Arbitrary Langrangian-Eulerian (A LE) method allows the combination between the Eulerian formu lation when the mesh is fixed and Langrangian formulat ion when the mesh moves with the matter to realize the clad laser.

Figure 3 shows the deformed mesh obtained with the ALE method under the nozzle and after it past. h is the clad laser height under the powder jet (Zone 1). The induced thermal field during the DLF process is determined solving the heat conduction equation. The fundamental behaviour of heat conduction is that a flu x Q (W/m3), representing the laser beam, heats the substrate surface.

The follo wing assumptions are taken into account in the fin ite element model:



The init ial temperature of the substrate is 300K



The heat source moves on the substrate surface at a velocity V  Thermophysical proprieties for the 316L powder and the low carbon steel substrate depend on the temperature.



The latent heat of fusion and vaporization are assumed as 218 kJ/kg and 7600 kJ/kg respectively.

The transient temperature distribution T(x,y,z,t) is obtained from the following equation:

Q T k t T C p = ∇ - ∇ + ∂ ∂ ) . ( . . ρ (3)
Where Cp is the specific heat capacity [J/kg.K], k is the thermal conductivity [W/m.K], t is the time [s].

The boundary condition on the surface area irradiated by laser beam (Zone 1) is given by:

) .( . ) ( ' ) . .( 4 0 4 0 T T T T h Q n T K c - - - - = ∇ - σ ε (4)
And on the substrate but outside the laser beam (Zone 2,3) th is condition is given by:

) .( . ) ( ) . .( 4 0 4 0 T T T T h n T K c - - - - = ∇ - σ ε (5)
The heat source is assumed to be Gaussian on the clad front: The density and thermal conductivity for the solid state are those of the substrate while they are those of the 316L powder for the liquid state. For the specific heat Cp both the solid and liqu id state values are those of the 316L stainless steel because it's main ly the powder that melt and solidify (Figure 4). The latent heat is taken into consideration through the specific heat expression by the following function:

) ( exp '
) , ( 2 ) ( dT T T hs flc C C C C f pl ps ps p - - + = 2 2 2 2 2 2 ) ( ) ( ) ( exp( ) ( ) ( ) ( exp( dT dT T T L dT dT T T L v v f f π π - - + - - + (9)
This effective exp ression is given in [START_REF] Jouvard | Thermal modelling of metal surface texturing by pulsed laser[END_REF]. C ps and C pl are the specific heat of solid and liquid state respectively, L f and L v are the latent heat of fusion and vaporization respectively, T v is the temperature of vaporization.

The first part of this expression ensures the latent heat variation between the solid and liquid states. In order to take into account that heat main ly diffuses in the substrate and considering that the main part of the liquid bath corresponds to the melted powder, C p is the solid substrate one and C pl the melted powder one. This way all nodes of the mesh follow the same expression and so there is no need to make time consuming and complex tests with the Arbitrary Langrangian-Eulerian (A LE) mode. The second part of this specific heat exp ression corresponds to the latent heat involves during the melting and vaporisation phenomena. With the same intention to do not slow down the computation speed and considering that the main part of the liquid bath corresponds to the melted powder L f and L v are the latent heat of fusion and vaporizat ion of the powder. The exponential formu lation is a smooth function numerically more suitable. The thermal simu lation of the laser cladding process shows a local melt ing pool. With the fin ite element simu lation on Co msol, an absorptivity of 75% for the laser beam radiation is found to reproduce the form and dimensions of the melt ing pool. A mong the 27 observed cross-sections, nine can be well simu lated by a half-cylinder. All these experimental conditions have been simu lated to validate our simulat ions.

Simu lating the 3D melting pool during the clad process is a real challenge because it is a mult i-physics complex problem. The melted depth in the substrate is a very sensitive element which can be chosen to discuss the obtained results. Indeed, if one considers for example a laser speed decrease two opposite phenomenon occur. In one hand the linear deposit energy increase and so the melted depth in the substrate is expected to increase but in the other hand the linear deposit powder mass increase so the clad height increase and the thermal source moves away fro m the substrate so the melted depth in the substrate is expected to decrease. This competition can lead to apparent contradictory results.

Experimental and simu lated depths of the melt pool are compared in table II.

For the lower powder feed rate (Table II, layer number a, b and c), experimental results show that the speed increase leads to a lower melted depth. The linear deposit energy phenomenon is predominant. In all other cases experimental and simu lated results show that the laser speed increase leads to a melted depth increase. In these cases the second phenomena is predominant, the linear deposit powder mass decreases so the clad height decreases and the thermal source comes closer to the substrate and so the melted depth in the substrate increases.

Cases b, c, g, h and i appear as very acceptable. In the cases a, e and f d ifferences between 70 and 100 micro meters are observed. Reported to the melted depth in the substrate, these differences are relatively important but if this modelling is used to predict the thermal field in the building part or in the substrate away fro m the melting pool the calculated thermal field can still correct. In case d the difference is real because no melting is predicted. The refined mesh dimensions could be insufficient to obtain the good result in this case.

Figure 5 shows a simulation of the thermal field in the clad during the process in case b. The moving heat source is applied at the end of the depositing clad. Choosing the suitable geometry fro m the beginning, the considered complex problem is limited to a thermal one. The thermal field determination can be solved in the 3D space with fast calculus time taking into account the real geo metry of the deposit clad. The delicate melted depth in the substrate is in most cases well simu lated and the melting pool geo metry can be simulated in 3D space. This paper presents a thermal simulat ion of the laser cladding process using the Arbitrary Langrangian-Eu lerian (ALE) mode to represent the laser clad deposition. The conduction equation is solved to calculate temperature distribution. In most cases results show a good correlation between the calculated and observed melt pool geo metries. The calculated thermal field could be used to optimize the process parameters in case of laser cladding process or direct metal laser fabrication. More directly, it helps to better understand the rule and intensity of the physical effects induced during the process as the laser speed effect on the melted depth. The next step is the prevision of the induced microstructures and constraints in the clad and the substrate.

Figure 1

 1 Figure 1 Model of direct laser fabrication

Figure 2

 2 Figure 2 Powder distributions in the jet

Figure 3

 3 Figure 3 Deformed mesh with ALE method under and outside the laser beam

  Where Q' is the flu x of energy [W/m2] Where hc is the heat convection coefficient [W/m2.K], where n is the normal vector o f the surface, ε is the emissivity, σ is the Stefan-Bo ltzman constant [5.67 × 10-8 W/ m2•K4] and T0 is the ambient temperature [K].The init ial condition for the transient analysis is

r

  is the variance of the Gaussian function.The material change of state is taking into account with a Heaviside function: the solid value, xl the liquid value, Tf the melting temperature, dT is the transition temperature between solid and liquid state, flc2hs is the Heaviside function.

Figure 4

 4 Figure 4 Specific heat with distinction between solid and liquid state and with latent heat of fusion Lf, and vaporization Lv V. RESULT AND DISCUSSION

Figure 5

 5 Figure 5 Melt pool view with the ALE method for the parameters (P=360 W, Qm=0.025 g/s and V=0.01 m/s)

Figure 6

 6 Figure 6 shows in two cases a comparison between the experimental melt pool and the simulated one.

Figure 6

 6 Figure 6 Cross sections of experimental results on the left and simulation results on the right showing the melt pool geometry in cases b and c.

Figure 7 Figure 7

 77 Figure7shows the simulated thermal field in the plane (y,z)

  3D Finite Element Simulation to Predict the Induced Thermal Field in Case of Laser Cladding Process and Half Cylinder Laser Clad Hussam El Cheikh *1 , Bruno Courant 2

TABLE I MATERIAL

 I PROPRIETIES FOR SOLID AND LIQUID STATE

		Densi	Conducti	Specific
		ty	vity	heat
		[kg/m	[W/(m.K)	[J/(kg.K
		3]	]	)]
	Solid state	7980	90	433
	Liquid state	7551	12.29	734

TABLE II EXPERIMENTAL

 II AND SIMULATED DEP TH OF THE MELTED POOL

	Laye r Numbe r	Power (W)	Powder feedrate (g/s)	Velocity (m/s)	Expe rimental de pth (mm)	Simulated de pth (mm)
	a	360	0.025	0.005	0.20	0.10
	b	360	0.025	0.010	0.15	0.17
	c	360	0.025	0.015	0.15	0.15
	d	360	0.050	0.005	0.13	0
	e	360	0.050	0.010	0.21	0.14
	f	360	0.050	0.015	0.22	0.14
	g	360	0.075	0.005	0.07	0.08
	h	360	0.075	0.010	0.14	0.10
	i	360	0.075	0.015	0.17	0.21


Pierre-Yves Hicher

Ching S Chang

P.-Y Hicher

An anisotropic non linear elastic model for particulate materials

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Particulate materials, such as soils or powders, can be considered as a collection of particles of different sizes and shapes. These discontinuous media can be represented by fictitious continuous media whose mechanical properties should depend on geometric arrangement and contact interactions between interacting particles. It has been experimentally shown that the elastic properties of such materials depend on many parameters, such as the material type, the stress condition, the void ratio, the coordination number or the packing structure ͑Hardin and Black 1966; [START_REF] Hardin | Shear modulus and damping in soils: Measurements and parameter effects[END_REF][START_REF] Iwasaki | Effects of grain size and grading on dynamic shear moduli of sands[END_REF][START_REF] Chung | Evaluation of dynamic properties of sands by resonant column testing[END_REF][START_REF] Lade | Modelling the elastic behavior of granular materials[END_REF][START_REF] Hicher | Elastic properties of soils[END_REF]. However, it is not easy to predict elastic moduli for particulate material. The main difficulty lies in the complex behavior of stress dependency leading to nonlinear behavior and stress-induced anisotropy. Furthermore, inherent anisotropic behavior occurs due to the mode of deposition. Recent developments in measuring techniques have allowed us to investigate more accurately the behavior of particulate media at very small strains, giving a more comprehensive insight into their elastic properties along different stress paths.

In previous works on granular mechanics, several empirical equations were proposed in the prediction of elastic moduli ͑Hicks and Monismith 1971; [START_REF] Hardin | Shear modulus and damping in soils: Measurements and parameter effects[END_REF][START_REF] Boyce | A nonlinear model for the elastic behaviour of granular materials under repeated loading[END_REF][START_REF] Brown | Analysis of pavement with granular bases, layered pavement system[END_REF][START_REF] Biarez | Elementary mechanics of soil behavior[END_REF][START_REF] Santha | Resilient modulus of subgrade soils: Comparison of two constitutive equations[END_REF]. These equations were derived from curve fitting to measured data. Therefore, the equations usually work well only in a limited range of stress conditions.

Another approach to the prediction of elastic moduli is based on micromechanics in which the granular material is considered as an assembly of particles. The moduli can be predicted directly from the properties of the constituent particles, for example, as in the work by [START_REF] Chang | Initial moduli of particulated mass with frictional contacts[END_REF][START_REF] Ng | Small-strain response of random arrays of spheres using discrete elements method[END_REF][START_REF] Ng | Small-strain response of random arrays of spheres using discrete elements method[END_REF] This approach provides a better understanding of the pressure sensitive behavior of granular media. However, under deviatoric stress conditions, the prediction is computationally intensive. It does not provide a straightforward equation as a prediction tool.

In this study, a constitutive model is developed from microscale behavior of interacting particles using a static hypothesis ͑Cambou et al. 1995; [START_REF] Chang | Kinematic and static hypotheses for constitutive modeling of granulates considering particle rotation[END_REF][START_REF] Liao | Stressstrain relationships for granular materials based on the hypothesis of best fit[END_REF], which relates the average stress of the granular assembly to a mean field of particle contact forces. Comparisons with discrete element models for particle assemblies have shown that the static hypothesis leads to representative results. Using a static hypothesis, the particle displacements are considered to fluctuate around a mean displacement field that represents the least square fit for the actual displacements ͑Liao et al. 1997͒. This leads to a more relaxed kinematic condition in comparison to the kinematic hypothesis.

The determination of the elastic constants of a granular assembly based on the properties of interparticle contacts will be discussed and the predictions of soil modulus under different loading conditions will be compared to experimental results obtained from different granular materials.

Micromechanical Approach

In recent years, several studies have been devoted to the elastic characteristics of granular materials within the micromechanical framework ͑Walton 1987; [START_REF] Chang | Micromechanical modeling of constructive relations for granular material[END_REF][START_REF] Jenkins | Volume change in small strain axisymmetric deformations of a granular material[END_REF][START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF][START_REF] Chang | Estimates of elastic modulus for media of randomly packed granules. II[END_REF][START_REF] Cambou | Stress-induced anisotropy during loading tests with and without stress[END_REF]. All these studies assume fundamentally the following two relationships: 1. Interparticle contact law-This provides the relationship between force and displacement at an interparticle contact; and 2. Link between micro-and macrovariables. There are two micro-macro relationships: a.

Relationship between interparticle contact forces and the assembly stress; and b. Relationship between interparticle contact displacements and the assembly strain.

Interparticle Behavior

The contact stiffness of an orientation includes normal stiffness, k n ␣ , and shear stiffness, k t ␣ , of the contact plane. The elastic stiffness tensor is defined in

f i ␣ = k ij ␣ ␦ j ␣ ͑1͒
where k ij ␣ = interparticle contact stiffness tensor. For two particles in contact, a local coordinate system can be constructed for each contact with three orthogonal base unit vectors: n is normal to the contact plane; s and t are tangential to the contact plane as shown in Fig. 1.

Let k n be the compressive contact stiffness in normal direction and k s be the shear contact stiffness. Assuming the shear contact stiffness is same in s and t directions and that there is no coupling effect between normal and shear directions, the contact stiffness tensor k qk c can then be expressed in terms of the unit vectors n, s, and t, as

k ij ␣ = k n ␣ n i ␣ n j ␣ + k t ␣ ͑s i ␣ s j ␣ + t i ␣ t j ␣ ͒ ͑ 2͒
For each particle contact, the corresponding auxiliary local coordinate system is related to the global coordinate system according to ͑see Fig.

1͒

n = ͑cos ␥, sin ␥ cos ␤, sin ␥ sin ␤͒ s = ͑-sin ␥, cos ␥ cos ␤, cos ␥ sin ␤͒ t = ͑0, -sin ␤, cos ␤͒ ͑ 3͒
The vector s is on the plane defined by x and n. The vector t is perpendicular to this plane and can be obtained by the cross product of n ϫ s. The rolling resistances between two particles are not discussed in this paper.

The value of the stiffness for two elastic spheres can be estimated from Hertz-Mindlin's formulation ͑Mindlin 1969͒. For sand grains, a revised form was adopted ͑Chang et al. 1989͒ given by

k n = k n0 ͩ f n G g l 2 ͪ n ; k t = k t0 ͩ f n G g l 2 ͪ n

͑4͒

where G g = elastic modulus for the grains; f n = contact force in normal direction; l = branch length between the two particles; and k n0 , k t0 , and n = material constants. For quartz mineral, the Poisson ratio is about 0.1 and k n0 can be approximated by

k n0 = ͱ 3/2G g l

͑5͒

For two spherical particles, the branch length is the particle diameter l = d. Using n =1/3, Eqs. ͑5͒ and ͑4͒ are equivalent to Hertz-Mindlin's contact formulation.

Macro-Micro Relationship

The stress-strain relationship for an assembly can be determined from integrating the behavior of interparticle contacts in all orientations. During the integration process, a relationship is required to link the macro-and microvariables. Using the static hypotheses proposed by Liao et al. ͑1997͒, we obtain the relation between the macrostrain and interparticle displacement ͑here, we do not consider the finite strain condition͒

u ˙j,i = A ik -1 ͚ ␣=1 N ␦ ˙j ␣ l k ␣ ͑6͒
where ␦ ˙j = relative displacement between two contact particles and the branch vector l k = vector joining the centers of two contact particles. It is noted that particle contacts include both direct contacts and indirect contacts of neighboring particles associated with a Voronoi polyhedron as discussed by Cambou et al. ͑2000͒.

For convenience, we let N be the total number of contact orientations. The variables ␦ ˙j ␣ and l k ␣ are defined, respectively, as the averaged values of ␦ ˙j and l k for all contacts belonging to the ␣th orientation. The fabric tensor in Eq. ͑6͒ is defined as

A ik = ͚ ␣=1 N l i ␣ l k ␣ ͑7͒
Using the principle of energy balance and using Eq. ͑6͒, the mean force on the contact plane of each orientation is

f ˙j ␣ = ˙ij A ik -1 l k ␣ V ͑8͒
In Eq. ͑8͒, the stress increment ˙ij can be obtained by the contact forces and branch vectors for contacts in all orientations ͑Christofferson et al. 1981; Rothenburg and Selvadurai 1981͒

˙ij = 1 V ͚ ␣=1 N f ˙j ␣ l i ␣ ͑9͒
Applying the defined contact force in Eq. ͑8͒, Eq. ͑9͒ is unconditionally satisfied.

Stress-Strain Relationship

The problem is defined as follows. Initially, we know the global variables ͑ ij and ij ͒ for the assembly and the local variables ͑f j ␣ and ␦ j ␣ ͒ for each contact orientation. For a given loading increment, which can be stress control, strain control or mixed mode, out of the 12 global variables ͑⌬ ij and ⌬ ij ͒, six of them are unknown. The objective is to determine all global vari- ␣ and ␦ j ␣ ͒ at the end of load increment. For a system with N interparticle orientations, the number of unknowns is 3N for f j ␣ and 3N for ␦ j ␣ . The total number of unknowns is 3N +3N +6.

The following constraints must be satisfied: 1. The local constitutive equation, i.e., Eq. ͑1͒. Since there are three equations for each contact plane orientation, the total number of equations is 3N, N being the total number of interparticle orientations. 2. Static hypothesis between global stress and local forces, i.e., Eq. ͑8͒: the number of equations is 3N. 3. Strain definition between global strain and local displacement, i.e., Eq. ͑6͒. The number of equations is 6 ͑strain is symmetric͒. The total number of unknowns is the same as the total number of equations. Therefore, the solution can be determined.

Using Eqs. ͑1͒, ͑6͒, and ͑8͒, the following relationship between stress and strain can be obtained:

u ˙i,j = C ijmp ˙mp
where

C ijmp = A ik -1 A mn -1 V ͚ ␣=1 N ͑k jp ep ͒ -1 l k ␣ l n ␣ ͑10͒
The summation in Eq. ͑10͒ can be replaced by an integral over orientations. The integral can lead to a closed-form solution for the elastic modulus of randomly packed equal-size particles ͑Chang et al. 1995͒. However, due to the nonlinearity nature of the local constitutive equation, a numerical calculation with iterative process is necessary to carry out the summation in Eq. ͑10͒ ͑Chang and Hicher 2005͒. In order to facilitate the numerical calculation, the orientations are selected to coincide with the locations of Gauss integration points in a spherical coordinate. Summation over these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over orientations. The results were more accurate by using a set of fully symmetric integration points. From a study of the performance of using different numbers of orientations, we found N ജ 74 to be adequate ͑Chang and Hicher 2005͒.

Particulate Materials with Isotropic Fabric

Isotropic Loading

Extensive studies of the elastic properties of soils have been done in the last 20 years, thanks to the development of new experimental techniques, which measure with a higher degree of accuracy very small deformations of particulate materials. The results have shown that elastic behavior exists only for very small strains, typically smaller than 10 -5 ; the elastic behavior is nonlinear, with a dependency upon both the actual state of stress and the actual void ratio.

We will at first assume an isotropic fabric for our particulate material. In these conditions, the parameters of the proposed model are as follows:

• Packing density ͑normalized number of contacts per unit vol-ume͒: = Nl 3 / V; • Mean particle size d; and • Interparticle elastic constants: G g , ␣, and n.

Note that the total number of contacts N per volume V is normalized by the cube of mean branch length l. The packing density is unitless. The branch length is defined as the length between centroids of two contact particles. For round particles, the mean branch length is approximately equal to the mean particle size d. The packing density for an assembly of equal-sized spheres can be related to the mean coordination number n ¯and the void ratio e of the assembly ͑Chang et al. 1989͒

= Nl 3 V = 3n ͑1 + e͒ ͑11͒
A relationship between and e obtained for a regular packing of spheres is plotted in Fig. 2.

For a representative volume consisting of a large number of round particles, in which the contact orientations are isotropically distributed, an analytical form of the shear modulus G and the Poisson's ratio of the granular assembly can be obtained ͑Liao

et al. 2000͒ G = 5k n ␣ 6d͑3 + 2␣͒ ͑12͒ = 1 -␣ 2 + 3␣ ͑13͒
where ␣ = k t0 / k n0 and = packing density. Poisson's ratio depends therefore only on the interparticle stiffness ratio ␣. Using Eqs. ͑4͒ and ͑5͒, under an isotropic stress c , the shear modulus G can be expressed as

G = H g 1-n ͩ c p 0 ͪ n

͑14͒

where p 0 = reference pressure ͑p 0 = 1 MPa in the computations in this paper͒ and H g = function depending on the elastic properties of the grains

H g = 5␣ 2 ͱ 6͑3 + 2␣͒ G g ͩ 3p 0 G g ͪ n ͑15͒
where G g is defined in Eq. ͑4͒.

Eq. ͑14͒ was derived from spherical grains. In a subsequent discussion, the relationship is extended to sands by empirically determining the packing density from experimental results. The influences of void ratio and confining stress on Young modulus have been experimentally studied for different granular materials ͑Hicher 2001͒. The main results are summarized in Fig. 3͑a͒. We can draw the following two conclusions: ͑1͒ for a given grain shape there is no influence of the grain size on elastic modulus. This implies that the modulus is dependent on the packing density ͑normalized by particle size͒, not on the particle size. ͑2͒ The modulus is dependent on both void ratio and grain size distribution; for a given void ratio, the Young modulus is higher for soils with uniform gradation ͑i.e., smaller coefficient of uniformity U c = d 60 / d 10 ͒. Thus the packing density is no longer a unique relationship with the void ratio e as it was shown in Fig. 2 for an assembly of equal-sized spheres. The grain size distribution should also be a factor influencing the packing density.

Fig. 3͑b͒ shows results obtained on different quartzic sands, including some studied in this paper. One can see that the general trend is very similar to Fig. 3͑a͒ and confirms the conclusions stated above. If we assume that the grain properties are roughly the same for all these quartzic sands, we can back-calculate, from the experimental results in Fig. 3, the packing densities and establish an empirical relationship between the packing density and the void ratio e for various coefficients of uniformity U c ͑Fig. 4͒. Mean particle size, d, is not an explicit parameter.

Thus, the parameters of the proposed model are as follows: • Packing density ͑obtained from void ratio e and coefficients of uniformity U c ͒; and • Particle shear modulus G g , stiffness ratio ␣, and exponent n.

One can therefore determine theoretically the shear modulus of an assembly of particles, knowing the elastic properties of the particles, the grain size distribution U c , and the void ratio e.

The approach that we adopted in the subsequent study was to determine the parameters of the model through the following procedure. From the test results, we can determine the ratio ␣ = k t0 / k n0 directly from the Poisson's ratio of the assembly ͓Eq. ͑12͔͒. The packing density can be determined from Fig. 4. The value of G g and the coefficient of nonlinearity n can therefore be determined by fitting predicted and experimental values obtained for the material modulus ͑E or G͒.

Fig. 5 shows the comparison of experimental and numerical results of the Young modulus as a function of void ratio and isotropic stress for glass ballotini ͑d 50 =1 mm͒ and Hostun sand ͑U c = 1.8, d 50 = 0.29 mm͒. The measured values of Poisson's ratio = 0.21 for both glass ballotini and Hostun sand ͑Hicher 1998͒, giving a value of ␣ = 0.36. The value of G g = 2,000,000 N / mm 2 and n = 0.5 was used for both glass ballotini and Hostun sand. By using the chart in Fig. 4 for packing density , the computed elastic modulus has a very good agreement with experimental results. Similar agreement was also found for the predicted and measured elastic moduli for Toyoura sand, Reid Bedford sand, Ottawa sand, and Ticino sand ͑the experimental results shown in Fig. 3͒. Results on Toyoura sand and Ticino sand will be discussed in the following sections. This demonstrated that elastic moduli can be predicted accurately for granular materials under an isotropic loading, provided that one can assume the initial isotropy of the packing. This hypothesis will be discussed.

Anisotropic Loading

When the applied state of stress is no longer isotropic, the model will predict a stress-induced anisotropy due to the dependency

␥ 31 · = ΄ 1 E 11 -21 E 22 -31 E 33 0 0 0 -12 E 11 1 E 22 -32 E 33 0 0 0 -13 E 11 -23 E 22 1 E 33 0 0 0 0 0 0 1 G 12 0 0 0 0 0 0 1 G 23 0 0 0 0 0 0 1 G 31 ΅ Ά 11 22 33 12 23
31 · ͑16͒

We will consider here only the case of a simple loading path: starting from an isotropic state of stress 0 , the stress in direction 1 is increased, while the stresses in directions 2 and 3 remain constant. After numerical integration of all interparticle stiffness, the elastic constitutive matrix is found to be cross anisotropic.

The evolution of the elastic constants is plotted in Fig. 6 as a function of the stress ratio 1 / 2 . All the components of the elastic matrix increase with the stress ratio, albeit at different rates. The material anisotropy, represented for example by the ratio E 1 / E 2 , increases with the stress anisotropy.

Particulate Materials with Anisotropic Fabric

One of the main advantages of the model is its capability of taking into account the structural anisotropy. In fact, the decomposition of the constitutive equations along a set of planes allows us to formulate the dependency of the parameters with respect to the orientation. We can thus obtain a response, which depends on the loading orientation in relation to the direction of material anisotropy. In the case of an orthotropic material, the anisotropy can be characterized by a fabric tensor, defined by

͓F ij ͔ = F ave ΄ 1 + a 0 0 0 0 1 - a 0 2 + 3a 22 0 0 0 1 - a 0 2 -3a 22 ΅ ͑17͒
where F ave = ͑F 11 + F 22 + F 33 ͒ /3. a 0 and a 22 = fabric constants describing the packing anisotropy ͑Chang et al. 1995͒. By representing the distribution of contact orientations with a continuous spherical harmonic expansion in three dimensions, it was found ͓see Chang and Misra ͑1990͔͒ that the fabric tensor matches the distribution described by a truncated form of the expansion consisting of second-order terms, given by

E ¯͑␥,␤͒ = F ave 4 ͫ 1 + a 0 4 ͑3 cos 2␥ + 1͒ + 3a 22 sin 2 ␥ cos 2␤ ͬ ͑18͒
The meaning of angles ␤ and ␥ is given in Fig. 1.

For a cross-anisotropic material, a 22 becomes equal to zero and Eq. ͑18͒ is reduced to

E ¯͑␥,␤͒ = F ave 4 ͫ 1 + a 0 4 ͑3 cos 2␥ + 1͒ ͬ ͑19͒
In three dimensions, the inherent anisotropy can be represented by a distribution whose major axis often coincides with the vertical direction. An example of distributions with different values of a 0 is shown in Fig. 7. The anisotropy of the packing structure is responsible to the anisotropy of the mechanical behavior.

When deposited in the gravitational field, granular materials comprise an anisotropic packing structure ͑or anisotropic fabric͒, often termed inherent anisotropy. Recent experimental studies ͑Belloti et al. 1996; [START_REF] Jiang | Inherent and stress-state-induced anisotropy in very small strain stiffness of a sandy gravel[END_REF][START_REF] Hoque | Anisotropy in elastic deformation of granular materials[END_REF] have demonstrated the influence of inherent anisotropy on the mechanical response of granular materials. The elastic matrix due to this inherent anisotropy can be considered as cross anisotropy. Assuming the vertical and horizontal axes are two principal material axes, Eq. ͑16͒ becomes

Ά ⌬ xx ⌬ yy ⌬ zz ⌬␥ yz ⌬␥ zx ⌬␥ xy · = ΄ 1 E h -hh E h -vh E v 0 0 0 -hh E h 1 E h -vh E v 0 0 0 -hv E h -hv E h 1 E v 0 0 0 0 0 0 1 G vh 0 0 0 0 0 0 1 G vh 0 0 0 0 0 0 2͑1 + hh ͒ E h ΅ ϫ Ά ⌬ xx ⌬ yy ⌬ zz ⌬ yz ⌬ zx ⌬ xy · ͑20͒

Isotropic Loading

In order to demonstrate the capability of the model to take into account an inherent anisotropy, three examples were selected

E v = ⌬ v ⌬ v and vh = - ⌬ h ⌬ v ͑21͒
When applying a small horizontal stress increment ⌬ h ͓i.e., ⌬ xx in Eq. ͑20͔͒ while keeping v constant ͓i.e., ⌬ zz =0 in Eq. ͑20͔͒, Eq. ͑20͒ leads to the following relations:

E h = ͑1 -hh ͒ ⌬ h ⌬ h and 2 hv 1 -hh = - ⌬ v ⌬ h ͑22͒
One can observe from Eq. ͑22͒ that the two equations are not sufficient to determine the three unknowns: E h , hv , and hh . The assumption made by Hoque and Tatsuoka ͑1998͒ was hh = vh under isotropic stress states. Furthermore, due to the limited accuracy of lateral strain measurements, Hoque and Tatsuoka believe that the measured values of hv were not reliable, and thus were not reported. This shows the difficulties of measuring with sufficient accuracy all the elastic constants for an anisotropic specimen.

In this study, the same assumption is made concerning hh . Although it leads to some uncertainty concerning the exact value of E h , it will not affect the discussion of model capabilities in reproducing the overall behavior of anisotropic granular materials.

Toyoura sand and Hime gravel are both uniformly graded materials with U c Ͻ 2. The values of vertical moduli E v ͑Fig. 8͒ are in very good accord with the results presented in Fig. 3. The values of horizontal moduli E h are invariably smaller for Toyoura sand, Hime gravel, and Chiba gravel. One can use the chart in Fig. 4 to determine the packing density corresponding to the void ratio of the tested specimens ͑e = 0.65͒. If we do so, a packing density = 2.5 is obtained for Toyoura sand.

A Poisson's ratio vh = 0.17 was measured during isotropic loading and the same value was assumed for hh . Estimated from Eq. ͑12͒, the corresponding value of ␣ = 0.45. The ratio between horizontal and vertical moduli E h / E v was found constant during isotropic loading and equal to 0.9, indicating the packing anisotropy. The anisotropy comes from many sources: the particle shape, the packing arrangement, the curvature of interparticle contact, etc. We have no information on these factors. In order for the model to capture this inherent anisotropy, an orientational dependent was assumed, which has the same material axis as the fabric. For a cross-anisotropic material, Eq. ͑19͒ can be used to define the orientational dependent ͑␥ , ␤͒, in which F ave = and a 0 = a . The anisotropy of the packing density a = 0.24 was determined in order to match the ratio E h / E v . Table 1 summarizes the parameters for Toyoura sand.

According to void ratios for Hime gravel and Chiba gravel ͑U c Ͻ 2͒, the values of are 4 and 4.2, respectively. The measured ratio between the horizontal and vertical moduli were E h / E v = 0.6 for Hime gravel and 0.45 for Chiba gravel. To match the E h / E v ratios, the anisotropy for the packing density was a = 1.03 for Hime gravel and 1.5 for Chiba gravel. The anisotropy for ␣ is also introduced in order to obtain anisotropic Poisson's ratio values in accordance with the measured ones ͑for Hime gravel, hh = vh = 0.15, and for Chiba gravel hh = vh = 0.25͒. For a cross-anisotropic material, Eq. ͑19͒ can be used to define the orientational dependent ␣͑␥ , ␤͒, in which F ave = ␣ and a 0 = a ␣ . Table 2 summarizes the model parameters for these two materials.

Comparisons between numerical and experimental results show that the model can represent accurately the influence of inherent anisotropy under isotropic loading. The results also demonstrate that the inherent anisotropy does not change during isotropic loading. In particular, the ratio between horizontal and vertical moduli remains constant. We will, in a second stage, study the evolution of the elastic constants during anisotropic loading.

Anisotropic Loading

The three previous materials were also subjected to different anisotropic loading. For Toyoura sand and Hime gravel, test programs corresponding to different stress paths lead to various values of the ratio h / v ͑Hoque and Tatsuoka 1998͒. The experimental results showing the evolution of the elastic properties with the state of stress are given in Fig. 9 along with the results of the numerical simulations. The parameters determined from isotropic loading were used for these calculations ͑Tables 1 and 2͒. The experimental results suggest that the Young moduli E v and E h are 9 that the evolution of E v and E h with v and h are also influenced by the value of the stress in the perpendicular direction ͑i.e., h for E v and v for E h ͒.

Numerical simulations show the same pattern and demonstrate the capability of the model to reproduce the observed behavior along anisotropic stress paths. One can notice, however, that the model gives a more pronounced influence of the transverse stress on vertical and horizontal Young moduli. This influence is due to the mode of integration along a set of planes of different orientations. Each plane contributes to the assembly behavior and, therefore, the stresses in all directions affect this behavior. Similar results were obtained on Chiba gravel along different stress paths. As a consequence, the stress induced anisotropy, which can be expressed by the influence of the v / h ratio on the E v / E h ratio, is less marked in the model compared to experimental evidence ͓Fig. 10͑a͔͒. The evolution of Poisson's ratio vh also shows an influence of the stress state which is captured by the model, although with less rate of change than the measured one ͓Fig. 10͑b͔͒. An extensive study of the elastic behavior of Ticino sand was performed by Bellotti et al. ͑1996͒. The tests were performed in a large calibration chamber and elastic properties were measured by means of propagating seismic body waves inside the specimens. With this technique, it is possible to measure wave propagation velocities in more than two directions ͑vertical and horizontal͒, which opens the way to determine all the elastic constants. In comparison with the previous measurement technique, the advantage is that all the elastic constants can be directly measured. The specimens were prepared by air pluviation, leading to an initial fabric which can be considered, as for the examples examined before, as cross anisotropic in vertical and horizontal axes. Several constant stress ratio tests were performed over a range of values of stress ratio K = h / v from 0.5 to 2. Results under isotropic loading conditions confirmed that the inherent anisotropy can be considered as cross anisotropic with isotropy in the horizontal plane since the wave velocities propagating in this plane were independent of the direction of propagation. The wave velocities either propagated along the z axis or polarized in the vertical plane were always lower than those confined in the horizontal plane. These last results lead to lower vertical than horizontal stiffness, contradicting the results discussed previously. Explanation for this phenomenon must await further study.

Considering the physical characteristics of Ticino sand: U c = 1.6 and d 50 = 0.55 mm, values of E v obtained along isotropic stress paths are in agreement with the results presented in Fig. 3. We can therefore apply the chart in Fig. 4 to determine = 1.2 for a given mean unit weight 1.505 MN/ mm 3 ͑relative density D r = 41%, e = 0.79͒. Values of E h / E v = 1.21 and n = 0.45 could be determined from the test results. Measured values of hh and hh were used to determine the values of ␣ in vertical and horizontal directions. Table 3 summarizes the set of model parameters.

Comparison between experimental results and numerical simulations is presented in Fig. 11͑a͒. An overall agreement was obtained for the different elastic constants. Constant values of Poisson's ratios are predicted by the model along the isotropic loading path, contrary to the experimental results which showed a slight increase of hh and a slight decrease of vh and hv . This could be the result of a certain evolution of the fabric during isotropic loading due to strain induced anisotropy, even if the authors estimated that strain induced anisotropy could be neglected in this study. This issue cannot be answered by the model in its present form. When examining experimental data obtained for different stress ratios, one cannot observe any clear trend in the evolution of elastic constants ͑see for example in Fig. 11, the evolution of Poisson's ratios for different stress ratios͒, which could also be due to a change of fabric anisotropy. Therefore, the model parameters governing the initial anisotropy were assigned for each loading condition, given in Table 4. The evolution of the packing density anisotropy a is an indicator of the evolution of the level of anisotropy with the stress ratio.

The results are presented in Fig. 11. Satisfactory agreement is found between experimental results and numerical simulations for moduli. However discrepancies are found in the comparisons for Poisson's ratio. One should be aware that Poisson's ratios for an anisotropic material cannot be experimentally determined with accuracy as discussed previously. Therefore, comparisons of measured Poisson's ratio with numerical simulations are difficult to assess. In any case, the comparisons show that the model captures the general trend for the evolution of elastic constants along stress paths with constant principal stress ratios.

Rotation of Principal Stress Directions

In all examples treated above, the principal stress directions were maintained constant and coincided with the axes of inherent anisotropy. We will now study the evolution of elastic constants along rotational stress paths. A recent study by Geoffroy et al. ͑2003͒ sheds new light on this evolution by means of a new apparatus, a torsional hollow cylinder device. Measurements of local strains by LVDT permit us to determine the following strain tensor expressed in the axes of initial anisotropy ͑vertical z and horizontal r, ͒:

͓͔ = ΄ rr 0 0 0 ␥ z 2 0 ␥ z 2 zz ΅ = ΄ r 0 0 0 ␥ 2 0 ␥ 2 z ΅ ͑23͒
Controlling axial, horizontal, and torsion shear stresses allowed different loading conditions to be applied with and without stress rotation

͓͔ = ΄ rr 0 0 0 z 0 z zz ΅ = ΄ r 0 0 0 0 z ΅ ͑24͒
Therefore, 16 terms of the flexibility matrix can be determined as follows:

Ά r z ␥ ͱ 2 · = ΄ M rr M r M rz M r␥ M r M M z M ␥ M zr M z M zz M z␥ M ␥r M ␥ M ␥z M ␥␥ ΅΄ r z ͱ 2 ΅ ͑25͒
Measurements of wave propagation velocities by bender and extender elements were used for this purpose in addition to LVDT measurements. The tests were performed on Hostun sand prepared by air pluviation. Two types of stress paths were applied on the specimens: a triaxial compression loading ͑Type C͒ and a torsion after axial loading corresponding to a ratio h / v = 0.5 ͑Type K͒. The authors did not give any data concerning the inherent anisotropy due to the mode of deposition and assumed that the specimens were close to initial isotropy. We will make the same assumption for the purpose of this study, even if the results presented in the previous section have demonstrated that the material is mainly cross anisotropic in the case of air pluviation preparation. In this condition, the parameters for Hostun sand were adopted.

In their paper, Geoffroy et al. present experimental values only for the two last columns ͑eight elastic constants͒ of the elastic matrix in Eq. ͑25͒. Comparisons between experimental results and simulations are presented in Fig. 12, where the elastic constants are normalized to the value under initial loading condition. For the triaxial loading case, the elastic constants are plotted against normalized stress z / z0 . For the torsional tests, the elastic constants are plotted against stress ratio / p mean . During triaxial loading ͑Case C͒ the terms in the last column M r␥ , M ␥ , M ␥z , and M z␥ should be null and stay null all along the loading. Despite some discrepancy in the measurements, the experimental data remain close to a mean value equal to zero in accordance with the model prediction ͑Fig. 12͒. During a torsion shear test, these four terms will evolve with the change in the principal stress directions. In the model, these two relations are always verified: M r␥ = M , M ␥z = M z␥ . Comparisons with experimental data show that the model is capable of capturing reasonably well the evolution of these four parameters during torsion, considering the dispersion of experimental data.

The evolutions of the terms M rz and M z during triaxial and torsion tests are also presented in Fig. 12. These two terms are always equal in the model prediction and one can see that, despite the dispersion of the experimental data, the model can predict correctly their evolution during triaxial tests, while their values remain practically constant during torsion shear tests.

Fig. 12 shows the evolution of the terms M zz and M ␥␥ . It can be seen that the model can capture the decrease of these two terms during triaxial loading, which corresponds to an increase in Young and shear moduli with stress ratio. During shearing, these two terms slightly increase with the shear amplitude, corresponding to a decrease in Young and shear moduli. Here this tendency can also be predicted by the model.

In conclusion, the model is able to take into account with reasonable accuracy the evolution of the elastic matrix due to stress induced anisotropy, which may be caused by either stress amplitude or rotation of principal stress axes.

Summary and Conclusion

An elastic model for particulate materials was developed based on micromechanics considerations. A particulate material is consid- ered as an assembly of particles. The stress-strain relationship for the assembly can be determined from integrating the behavior of interparticle contacts in all orientations using a static hypothesis which relates the average stress of the granular assembly to a mean field of particle contact forces. This stress-strain relationship leads to an elastic nonlinear behavior of the particulate material.

The elastic coefficients at the level of the assembly depend on the parameters used to define the contact law at the level of the particles. These parameters comprise the contact law itself, i.e., the normal stiffness, k n0 , and the shear stiffness as a function of the normal stiffness, k r0 = ␣k n0 , as well as the packing density ͑normalized number of contact per unit volume͒ along a given plane orientation. The dependency of the normal stiffness on the normal force gives a pressure-dependent behavior for the assembly, while the interparticle stiffness ratio ␣ determines the value of Poisson's ratio. The influence of the void ratio of the assembly on the elastic properties is taken into account through the packing density . An empirical relationship was proposed in order to link values of to the void ratio. This relationship was based on results obtained for different granular materials with different grain size distributions.

Inherent anisotropy due to the mode of deposition can be assessed by taking into account a dependency of the parameters defining the normalized number of contacts in a given plane on the orientation of this plane. The nonlinear character of the contact stiffness leads to a stress-induced anisotropy with regard to the overall behavior of the grain assembly. This stress induced anisotropy includes the effects of the principal stress amplitudes as well as the rotation of the principal stress directions.

The predictions of the model for different stress conditions were compared to the results obtained during several experimental studies on different granular materials. The comparisons showed that the model is capable of taking into account very precisely the influence of the inherent anisotropy. The stress induced anisotropy along the different stress paths is also taken into account, but the model allows the transverse stresses to have a more pronounced influence on Young and shear moduli. This is the result of the way the local equations are integrated. However, it must be pointed out that this discrepancy becomes noticeable for elevated stress ratios. In this regard, several studies have demonstrated that the change in the fabric due to strain-induced anisotropy could no longer be neglected ͑Kuwano and Jardine 2002; Dano and Hicher 2003͒. Under these conditions, a mechanism of strain-induced anisotropy needs to be introduced in the model, which was not the subject of this study. This could be achieved in future studies, given that an elastoplastic model based on the same numerical concepts as the ones exposed here has recently been developed ͑Chang and Hicher 2005͒, which allows nonreversible strains to be computed that could be used to take into account the fabric evolution and the related evolution of the elastic properties.

Fig. 1 .

 1 Fig. 1. Local coordinate at interparticle contact

Fig. 2 .

 2 Fig. 2. Relationship between void ratio and coordination number

Fig. 5 .Fig. 6 .Fig. 3 .

 563 Fig. 5. Comparison between measured and predicted values of Young's modulus for glass ballotoni and Hostun sand

Fig. 7 .

 7 Fig. 7. Distribution of contact orientations with different values of a 0

Fig. 8 .

 8 Fig. 8. Isotropic loading on anisotropic Toyoura sand, Hime gravel, and Chiba gravel

Fig. 9 .Fig. 11 .

 911 Fig. 9. Stress dependency on the horizontal and vertical Young's modulus for Toyoura sand and Hime gravel

Table 1 .

 1 Elastic Parameters for Toyoura Sand ͑e = 0.65͒

		a	G g	␣	n
	2.5	0.24	2,100,000	0.45	0.5

Table 2 .

 2 Model Parameters for Hime Gravel and Chiba Gravel dependent on the vertical and horizontal stress v and h respectively. Other experimental works ͓see for example Belloti et al. ͑1996͒ and Kuwano and Jardine ͑2002͔͒ reached the same conclusions. However, one can see in Fig.

	Material		a	G g	␣	a ␣	n
	Hime gravel	0.53	1.03	2,100,000	0.56	0.05	0.5
	͑e = 0.5͒
	Chiba gravel	0.36	1.5	2,100,000	0.33	-0.07	0.52
	͑e = 0.49͒

mainly

Table 4 .

 4 Sets of Parameters for Ticino Sand for Different Stress Ratios

	K		a	␣	a ␣	n
	0.5	0.32	0.11	0.5	0.48	0.45
	1.5	0.34	-0.65	0.48	0.67	0.45
	2.0	0.40	-0.78	0.38	0.23	0.45

Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes, Univ. of Nantes, Nantes, France.

Professor, Dept. of Civil and Environmental Engineering, Univ. of Massachusetts, Amherst, MA 01002 ͑corresponding author͒. E-mail: chang@ecs.umass.edu