Xuan Wang

Valérie Pommier-Budinger

Yves Gourinat

Aurélien Reysset

A Modified Preisach Model for Hysteresis in Piezoelectric Actuators

Keywords: hysteresis, Preisach model, piezoelectric actuator, damper I

Piezoelectric actuators (PEAs) exhibit hysteresis nonlinearity in open-loop operation, which may lead to undesirable inaccuracy and limit system performance. Classical Preisach model is widely used for portraying hysteresis but it requires a large number of first-order reversal curves to ensure the model accuracy. All the curves may not be obtained due to limitations of experimental conditions, and the detachment between the major and minor loops is not taken into account. This paper aims to propose a modified Preisach model that demands relatively few measurements and that describes the detachment. The modified model is implemented by adding a damper in parallel with the classical Preisach model. The parameter of the damper is adjusted to an appropriate value so that the measured and predicted hysteresis loops are in good agreement. Experimental results prove that the proposed modified Preisach model can characterize hysteresis more accurately than the classical model.

INTRODUCTION

Piezoelectric actuators (PEAs) are commonly used in a variety of applications [START_REF] Ghasemi-Nejhad | Finite element method of active vibration suppression of smart composite structures using piezoelectric materials[END_REF][START_REF] Mokaberi | Compensation of scanner creep and hysteresis for AFM nanomanipulation[END_REF][START_REF] Gaul | Model-based piezoelectric hysteresis and creep compensation for highly-dynamic feedforward rest-to-rest motion control of piezoelectrically actuated flexible structures[END_REF][START_REF] Agrawal | Modeling and Control of Piezoelectric Embedded Elastic Plates[END_REF]. However, PEAs also exhibit adversary hysteresis effect in open-loop operation. Hysteresis, which is observed between the voltages applied to the actuator and the obtained displacement as seen as in Fig. 1, is a nonlinearity that the obtained displacement depends on the present input voltage as well as on how the inputs were applied previously. Hysteresis severely limits system performance, since it leads to a positioning uncertainty of up to 15% of motion range [START_REF] Moriyqmq | Development of a precision diamond turning machine for fabrication of asymmetric aspheric mirrors[END_REF].

To compensate hysteresis by open-loop control, a model precisely modeling hysteresis is needed. Several models have been developed and reported. Amongst those models Preisach [START_REF] Ge | Modeling hysteresis in piezoceramic actuators[END_REF][START_REF] Ge | Generalized preisach model for hysteresis nonlinearity of piezoceramic actuators[END_REF] and Prandtl-Ishlinskii (PI) [START_REF] Janaideh | Generalized Prandtl-Ishlinskii hysteresis modeling and its inverse for compensation in smart actuators[END_REF] models are formed by a weighed superposition of infinite elementary hysteresis operators. Duhem model [START_REF] Visintin | Differential Models of Hysteresis[END_REF] is in the form of nonlinear differential equation. These models do not describe the physical mechanism of hysteresis, thus then are referred to phenomenological models. In comparison, the model proposed in [START_REF] Yeh | Modeling and identification of hysteresis in piezoelectric actuators[END_REF] is constructed by physical consideration in which the basic physical elements including springs, massless carts, and Coulomb frictions are employed to portray hysteresis.

Classical Preisach model is widely used for portraying hysteresis, but the model accuracy is dependent on the number of measured first-order reversal curves of hysteresis loop and, due to the limitations of experimental conditions, it is possible that the number of obtained curves is not enough to accurately model hysteresis. Moreover, the detachment between the major and minor loops is not taken into account in the classical Preisach model. To overcome these limitations, this paper proposes a modified Preisach model to accurately characterize hysteresis in PEAs. This model requires relatively few measured first-order reversal curves and describes the attachment. In the proposed model, a damper is added in parallel with the classical Preisach model. The parameter of the damper is adjusted to an appropriate value so that the measured and predicted hysteresis loops are in good agreement. Experimental results show a significant reduction of the relative error between the measured and predicted hysteresis response curves by using the modified Preisach model.

II. CLASSICAL PREISACH MODEL

A. Model description and implementation

The basic idea of the Preisach model lies in the description of the hysteresis through an infinite number of operators .

Consider a simple two-position relay operator , in which the input-output relationship at each instant of time t as depicted in Fig. 2 can be described by: 0 ( )

( ( )) 1 ( ) ( ) ( ) if u t u t if u t u t if u t ≤ = ≥ < < (1)
where numbers and correspond to "up" and "down" switching values of input, respectively (assumed in the sequel that ≤ ). This research is supported by "the Fundamental Research Funds for the Central Universities" of Civil Aviation University of China (CAUC-ZXH2012D009). Along with the relay operator , consider a weighting function ( , ) . Then the mathematical form of the Preisach model can be written as follows [START_REF] Mayergoyz | Mathematical Models of Hysteresis and Their Application[END_REF]:

( ) ( , ) ( ( )) x t u t d d ≥ = (2)
where ( )

x
Assume that at an initial time 0 t , the input 0 0 ( ) 0 u t u = < . This means that the output of every relay operator is 0, and 0 ( ) P t P -= , 0 ( ) 0 P t + = (Fig. 3(b)). Next as the input is monotonically increasing to some maximum value at is generated as part of boundary (Fig. 3(d)). Then additional horizontal or vertical boundary segments are generated as input reversals.

From above analysis, it can be seen that the output of the Preisach plane is determined by the boundary link between P - and P + . Due to its staircase structure, the boundary link is fully captured by all corner points, which correspond exactly to the past dominant maximum and minimum input values and that are stored and updated in numerical implementation of the Preisach model.

A displacement function ' '
( , ) is introduced in the identification of Preisach model, as follows:

' ' ' ' ' ( , ) ( ) ( , ) x x = - (4)
where -' '

( , ) represents the change in the output ( ) x t as the input ( ) u t changes from ' to ' ; -' ( )

x
is the output value on the limiting ascending branch of the first-order reversal curve and -' '

( , )

x
is the output value on the descending curve, which is attached to the limiting ascending branch at the point ' ( ) x (Fig. 4). ' and ' represent the maxima and the minima of the input ( ) u t .

A geometrical interpretation for ' ( ) x and ' '

( , ) x can be illustrated in the Preisach plane shown in Fig. 5. It can be seen that the result of the input ( ) u t decreasing from ' to ' leads to subtracting a triangle ' '

( , ) T from the region P + . Thus Equation ( 4) can be written as:

' ' ' ' ( , ) ( , ) ( , )
T d d = (5)
If the hysteresis loop contains several extrema, the region P + of the Preisach plane is composed of several trapezoidal parts k P , as shown in Fig. 6. All the extrema ' k and ' k that depend on the past values of the input voltage ( ) u t are stored in the history. For the part 1 P , the following equation is deduced:

1 ' ' ' ' 1 0 1 1 ( , ) ( , ) ( , )
P d d = - (6)
The other parts are calculated in the same way. Because the integration on the region P + is the sum of the integrations on all the parts k P , the total displacement of the PEA ( ) x t for an input voltage ( ) u t is determined thanks to (3), depending on the current slope of ( )

u t : ( ) 0 u t > ' ' ' ' ' 1 1 ( ) [ ( , ) ( , )] ( ( ), ) N k k k k N k x t u t - = = - + (7)
( ) 0 u t < 1 ' ' ' ' ' ' ' 1 1 1 ( ) [ ( , ) ( , )] ( , ) ( , ( )) N k k k k N N N k x t u t - - - = = - + - ( 8
)
where N is the number of maxima ' k and minima ' k that are stored.

In order to compute the values ' '

( , ) , a mesh of and is created within P . The reference value ( , ) , are measured on the PEA for all and of the mesh and stored at each corresponding node as depicted in Fig. 7 ( , )

i j + +
as shown in Fig. 7. The output displacement is determined using either [START_REF] Ge | Generalized preisach model for hysteresis nonlinearity of piezoceramic actuators[END_REF] or [START_REF] Janaideh | Generalized Prandtl-Ishlinskii hysteresis modeling and its inverse for compensation in smart actuators[END_REF].

The wipe-out property is needed in the Preisach model. It allows erasing the pair ' ' 1 ( , )

N N -
from the history once ( ) u t exceeds ' N . Similarly, the pair ' ' ( , )

N N can be erased from the history once ( ) u t becomes smaller than ' N . This avoids the excessive growing of the stored values.

From above procedure of identification, it is obviously seen that the accuracy of the model depends on the mesh density of the Preisach plane. The denser mesh means more data pairs to collect and store, which requires more measured first-order reversal curves.

B. Experiments

The model implementation and validation need many experiments that are realized in open-loop. The main components of experimental setup are shown in Fig. 8. A ringtype piezoelectric actuator (NCE41, Noliac) with 400V maximum operating voltages is used as study object. The displacement sensor is a noncontact capacitive sensor (D510.100, PI) enabling the capture of the PEA performances. Fig. 9 shows a functional diagram of experimental system.

To implement the classical Preisach model, a set of piecewise sinusoidal voltages of amplitude range from 10 V to 400 V with 10 V increment (Fig. 10) is applied to the PEA. A frequency of 0.25 Hz is chosen in the experiment so that neither creep nor dynamic effect of piezoelectric actuator affects the output response. In total, 861( 41 ( 41 curves computed by using classical Preisach model do not conform very well to the measured curves in the PEA and that, for both cases, the predicted hysteresis loops are thinner than the measured loops. Fig. 11(c) and Fig. 12(c) depict the relative errors between the measured and predicted responses. It is observed that the maximum relative errors for the two cases respectively approach to 4% and 3.5%.

From Fig. 11(a) and Fig. 12(a), it can also be seen that the predicted hysteresis loop with sinusoid input voltages of amplitude 200 V is closer to the measured loop than the case with sinusoid input voltages of amplitude 100 V. This is because more data pairs are employed for the model computation.

In addition to the influence of the number of measured first-order reversal curves on the model accuracy, another factor is analyzed. Fig. 13 shows major and minor hysteresis loops respectively generated by exciting the PEA with positive sinusoid input voltages of 0.25 Hz frequency and of 200 V and 100 V amplitudes. Fig. 13(a) depicts that the measured hysteresis curve of the minor loop detaches inward from that of the major loop. Hence the corresponding data pairs of these inward detachments are collected and stored for the implementation of classical Preisah model. However, the minor loop attaches to the major loop on the basis of classical Preisach model as depicted in Fig. 13(b). That is why the predicted hysteresis loops look like thinner than the measured loops. Concerning the detachment, it is observed that, while the time rate of the input for the minor loop is zero, the rate for the major loop is nearly at its peak. For this reason, it is postulated that the detachment is caused by rate-related dynamics [START_REF] Yeh | An integrated physical model that characterizes creep and hysteresis in piezoelectric actuators[END_REF]. That is why we propose to add a damper in parallel with the classical Preisach model to describe the detachment and then to get a more accurate model of the hysteresis. 16(c), it is observed that the maximum relative errors for both cases reach to 2%. Compared to Fig. 11(c) and Fig. 12(c), 50% and 42% of the errors vanish respectively. The good agreement found between the measured and predicted curves shows that the modified Presach model can describe the hysteresis in PEAs more precisely and do not require more measurements of the first-order reversal curves.

Finally, Fig. 17 shows the comparison between the measured hysteresis loops and the predicted loops computed by using the modified Preisach model for both studied cases. It is obviously seen that the detachment between the major and minor loops can be precisely described with the modified Preisach model. The damper is used for plumping the shape of the predicted hysteresis loop to conform to the measured loop and for modeling the detachment that may be caused by rate-related dynamics. Experimental results show that the maximum relative errors between the measured and predicted hysteresis responses is reduced by 50% and by 42% respectively for two studied cases by using the modified model. It is evident that the proposed modified Preisach model is an effective mean to accurately describe hysteresis.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Hysteresis in piezoelectric actuators Fig. 2. An elementary two-position relay operator

Fig. 3 .and 3 a

 33 Fig. 3. Memory curve in the Preisach plane

  ) measured from the experiment are collected and stored in the computer. Then the classical Preisach model is implemented.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Discretization of the Preisach plane P into a finite number of rectangles

  predicted

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. A sine piecewise voltage signal used to identify classical Preisach model

mD

  is tuned incrementally from zero until the error between the predicted and measured loops is minimal, as shown in Fig.15(a). Here m D is chosen to 0.003. The case with sinusoid input voltages of amplitude 100 V is still studied here and is compared with the reference case of amplitude 200 V. Fig.15(a)-(b) and Fig. 16(a)-(b) illustrate the measured and predicted hysteresis loops and time responses for both cases. Figures show that the predicted curves computed by using modified Preisach model conform exactly to the measured curves of the PEA for both cases. It means that the value of m D obtained from the reference case can be used for other input voltage amplitudes. From Fig. 15(c) and Fig.

FigFig. 15 .Fig. 16 .

 1516 Fig. 14. Modified Preisach model

Fig. 17

 17 Fig. 17. The detachment in measured and predicted (by using modified Preisach model) hysteresis loops

  Fig. 17. The detachment in measured and predicted (by using modified Preisach model) hysteresis loops IV. CONCLUSION This paper proposes a modified Preisach model that can more accurately portray hysteresis in PEA than classical Preisach model. The proposed model is based on relatively few measured first-order reversal curves and can also characterize the detachment phenomenon between the major and minor loops. The presented modified model is implemented by adding a damper in parallel with the classical Preisach model.The damper is used for plumping the shape of the predicted hysteresis loop to conform to the measured loop and for modeling the detachment that may be caused by rate-related dynamics. Experimental results show that the maximum relative errors between the measured and predicted hysteresis responses is reduced by 50% and by 42% respectively for two studied cases by using the modified model. It is evident that the proposed modified Preisach model is an effective mean to accurately describe hysteresis.

  t is the output displacement of a PEA. The mechanism of memory formation in the Preisach model is usually discussed in the Preisach ( , )

	plane,	P	{( , )	∈	P	|	≥	,	≤	0	,	0} ≥ , as shown in Fig.
	3(a). Each ( , ) P ∈ corresponds to an operator	. It is
	assumed that the weighting function ( , ) is equal to zero
	outside of P . At each time instant, P is divided into two parts:
			P -	{( , )	∈	P	| output of	at t is 0}	,
			P +	{( , )	∈	P	| output of	at t is 1}	,
	so that	( ) P P t P t ( ) -+ =	at any time. Equation (2) can be
	replaced by:
						( ) x t			( , )	d d
									( ) P t

+ =

ACKNOWLEDGMENT

This research is supported by "the Fundamental Research Funds for the Central Universities" of Civil Aviation University of China (CAUC-ZXH2012D009).


S Dubreuil

F Petitjean

M Salaün

A coupled parametric and nonparametric approach for modal analysis of a satellite

This study takes place in the context of dynamical prediction for satellite structures. Aims of such studies are to survey the dynamical response of satellite equipment and components, to check that requirements are correct and to give prediction of vibration levels, which are inputs for the experimental test validation according to the launcher specification. This prediction is done by modal analysis performed on a numerical model built by finite element method. Uncertainties on equipment and components properties lead to random frequency response function (FRF). This paper aims at understanding how modal approach can be adapted to probabilistic framework in order to calculate cumulative density function or, at least, some quantiles of a FRF. Starting point of deterministic modal analysis is the study of a single degree of freedom (DOF) system. C. Heinkelé analytically expresses the probability density function (PDF) of the FRF of an oscillator with a random natural pulsation following a uniform law. We have generalized this work to random natural pulsation following a law of finite variance. Then, the expression of a FRF between two DOF of the structure is a linear function of random oscillators FRF and DOF components of random eigenvectors. Assuming that random eigenvectors are close to their means, we have access to the characteristic function of the random FRF between two DOF as a multi-dimensional integral with respect to the joint PDF of the oscillators FRF. This paper mainly focuses on two major points which are calculation of oscillators joint PDF and inversion of characteristic functions. The first one is tackled by copulas theory. Dependence structure of random eigenvalues are identified and modeled by a copula. Then we apply results on copulas transformations to obtain joint PDF of oscillators FRF. Classical results concern monotonic transformations but we extended these ones to non-monotonic cases. Concerning inversion of characteristic function, several methods are studied to numerically compute the Gil-Pelaez formula. This approach allows to access some FRF quantiles by numerical integration which error can be controlled. Moreover, an interesting point is the flexibility in the identification of the random eigenvalues PDF. This is especially interesting in order to couple parametric identification with nonparametric one, when only few dispersion informations are given for equipments, housed in the satellite primary structure.

OSCILLATOR WITH RANDOM NATURAL PULSATION

Starting point of modal analysis is the single DOF system. In deterministic context its transfer function is

H(ω) = -ω 2 + λ (-ω 2 + λ) 2 + 4ξλω 2 + j -2ξ √ λω (-ω 2 + λ) 2 + 4ξλω 2
where j 2 = -1, ω is the pulsation, ξ the critical damping coefficient and λ the square of the natural deterministic pulsation. Uncertainties on oscillator prop-erties (stiffness, mass) lead to random natural pulsation. Let λ be the random variable modeling the square of random natural pulsation, p λ its PDF and H(ω) the associated transfer function. Following results presented in [START_REF] Heinkelé | Synthèse modale probabiliste[END_REF], we calculate PDF of transfer function presented above (difference with the work in [START_REF] Heinkelé | Synthèse modale probabiliste[END_REF]) is critical damping model and generalization to λ a second order random variable). As methodology for real, imaginary part or module and phase is similar, we illustrate our study on real part of transfer function, which expression reads F : R → 1 4ξω 2 (ξ -1) ; 1 4ξω 2 (ξ + 1)

F(λ ; ξ, ω) = Re(H) = -ω 2 + λ (-ω 2 + λ) 2 + 4ξ 2 λω 2 (1)
It should be noted that this transformation is nonmonotonic and have two inverse images called F -1 1 , F -1 2 . Then, one can expresses PDF of Re(H(ω)) using classical results on random variable transformations. This leads to:

p Re(H) (Re(h)) = p λ (F -1 1 (Re(h))) |J Re (F -1 1 (Re(h)))|
(2)

+ p λ (F -1 2 (Re(h))) |J Re (F -1
2 (Re(h)))| where J Re is the jacobian of transformation (1). Integration of this density is possible since λ follows a second order probability law, which is an acceptable assumption in our context. Numerical integration of this PDF allows to access mean, standard deviation, and quantiles of Re(H). Figure 1 shows various shapes of this function for different ω values when p λ follows a normal law.

-0.0004 -0.0003 -0.0002 -0.0001 0.0000 0.0001 0.0002 0.0003 0.0004

Re(h(ω)) Two difficulties arise during numerical integration. First one is the infinite limit at edges of domain. This problem is solved by using numerical integration on [ 1- 4ξω 2 (ξ-1) , 1- 4ξω 2 (ξ+1) ] and by analytical integration of an equivalent of the density on [

1 4ξω 2 (ξ-1) , 1- 4ω 2 (ξ-1) ] and [ 1- 4ξω 2 (ξ+1) , 1 4ξω 2 (ξ+1) ].
Second one concerns the shape of the density when ω is far from E[λ]. As it can be seen on figure 1 for ω = 100 and ω = 340 the density looks like a single peak with a very strong slope. In fact, density tends to a zero Dirac distribution. Thus, to have a correct numerical integration, discretization of the density must be done according to its derivative.

Figure 2 shows results of this methodology in determination of 2.5%, 50% and 97.5% quantiles in comparison with Monte-Carlo method (one million runs). Numerical integration on regular domain is performed by trapezoidal rule, number of points depends on the pulsation ω, for reasons explained previously. Results presented by figure 2 allow to be confident in this integration methodology.

In order to adapt deterministic modal analysis to probabilistic framework, we now present a generalization of these results to a multi-degree of freedom system.

ADAPTATION TO MULTI DOF

Modal synthesis

Let us consider a finite element model of a mechanical structure. We denote [M], [K], [D] the mass, stiffness and damping matrices of the mean model (model for which input parameters are set to their means). In frequency domain, vibration equation reads

(-ω 2 [M] +jω [D] + [K]) y(ω) = f (ω)
in which y(ω) is the vector of the DOF and f (ω) is the vector of loads input to the system. Modal synthesis is based on solving following eigenvalues problem

(-ω 2 [M] + [K]) y = 0.
Aims of modal analysis is to select a basis [ϕ] of eigenvectors which dimension N is far less from the number of DOF. We denote by λ k , k = 1 • • • N , corresponding eigenvalues. This approach relies on assumption that low frequency comportment of structure can be expressed on a reduced basis composed by first eigenvectors. In this basis, mass and stiffness matrices are diagonal and a current assumption is to consider that the damping one is diagonal too. This leads to a diagonal transfer functions matrix

[H] kk (ω) = -ω 2 + λ k (-ω 2 + λ k ) 2 + 4ξ k λ k ω 2 + j -2ξ k λ k ω (-ω 2 + λ k ) 2 + 4ξ k λ k ω 2
which is an oscillator equation described in the first part. Then, FRF between two DOF a and b of a structure is:

[F ] ab = N i=1 [H] ii [ϕ] bi [ϕ]
Φ Re(F ab ) (v) = R N exp(jv N i=1 Re(h i )[ϕ] bi [ϕ] ai )) p Re({H}) (Re({h}))dRe({h}) (3)
Once calculation of this characteristic function is done, inversion formulas like Gil-Pelaez one [START_REF] Gil-Pelaez | Note on the inversion theorem[END_REF] or Poisson summation [START_REF] Abate | The Fourier-series method for inverting transforms of probability distributions[END_REF] allows to access CDF or PDF of Re(F ab ). As equation ( 3) is a multi dimensional integral respectively to the joint density of Re({H}), it is of great importance to study dependence structure of this random vector.

Independant case

The simplest assumption is to consider that random variables Re(H k ) for k = 1, • • • , N , of Re({H}) are independent and of marginal density described in the first section. Then, equation (3) becomes a product of N integrals in one dimension

Φ Re(F ab ) (v) = N i=1 exp(jv[ϕ] bi [ϕ] ai Re(h i )) p Re(H i ) (Re(h i ))dRe(h i ), which leads to Φ Re(F ab ) (v) = n i=1 1 |[ϕ] bi [ϕ] ai | Φ Re(H i ) (v)
with an appropriate change of variable. Nevertheless, assumption of independent eigenvalues seems unrealistic. For example, let us introduce the finite element model of satellite TARANIS presented by figure 3.

Input properties of this structure are modeled by independent random variables (input properties are materials ones like Young modulus and Poisson coefficient and geometric ones like thickness of panels etc.). Figure 4 presents dispersion of the two first eigenvalues coming from 1000 runs of this random model. It is clear that these two random eigenvalues are correlated, so modeling and consequence on FRF of this correlation must be investigated.

Dependence model

In this part, we set N = 2 in order to simplify expressions but generalization is straightforward. Moreover, in practice the case N = 2 should be frequent as correlation of eigenvalues needs to be modeled only if oscillators are acting on same pulsations range. We are interested in modeling the joint probability law of random vector Re({H}) whose marginal laws are described in previous section. To build this model we propose to use copulas theory introduced in [START_REF] Nelsen | An introduction to copulas[END_REF]) and (Clauss ). Indeed, a result known as Sklar theorem allows to build joint PDF of random eigenvalues, and to express it as,

p λ 1 ,λ 2 (λ 1 , λ 2 ) = p λ 1 (λ 1 )p λ 2 (λ 2 )c(F λ 1 (λ 1 ), F λ 2 (λ 2 ))
where c is the density of the copula modeling the dependence structure of eigenvalues and F i the marginal CDF of λ i . In order to express joint density of Re({H}), we use results on random variable transformation and copula transformation presented in [START_REF] Nelsen | An introduction to copulas[END_REF]. These last ones are derived for monotonic transformation. Non monotonic application presented here is the principal difficulty.

Let Λ be the random vector [λ 1 , λ 2 ] of joint density p λ 1 ,λ 2 (λ 1 , λ 2 ) described above and f the function of R 2 with values in

1 4ξ 1 ω 2 (ξ 1 -1) ; 1 4ξ 1 ω 2 (ξ 1 +1) × 1 4ξ 2 ω 2 (ξ 2 -1) ; 1 4ξ 2 ω 2 (ξ 2 +1) . f λ 1 λ 2 = Re(H 1 ) = F(λ 1 ; ξ 1 , ω) Re(H 2 ) = F(λ 2 ; ξ 2 , ω)
where F is the function described by equation ( 1). Finally, one accesses density of Re({H}) in function of previous transformation, marginal density of eigenvalues and copula modeling dependence structure of eigenvalues:

p Re({H}) (Re({h})) = 2 i,j=1 p λ 1 (F -1 i (Re(H 1 ))) |J(F -1 i (Re(H 1 )))| × p λ 2 (F -1 j (Re(H 2 ))) |J(F -1 j (Re(H 2 )))| × c F 1 (F -1 i (Re(H 1 ))), F 2 (F -1 j (Re(H 2 ))) (4)
This density is used to numerically calculate the characteristic function (3).

NUMERICAL RESULTS

Dependence structure of Re({H})

In order to illustrate our methodology, we propose to consider two random oscillators with dependent random eigenvalues :

• λ 1 follows a normal law of mean 220 2 and a variation coefficient of 15%.

• λ 2 follows a triangular law defines on [50000, 75000] with a maximum at 59375

• Eigenvalues dependence is modeled by a normal copula of Kendall correlation coefficient equal to 0.95. Figure 5 shows shapes of marginal densities of λ 1 and λ 2 and figure 6 illustrates joint density of eigenvalues.

Note that in an industrial application uncertainties come from input parameters of model, then construction of eigenvalues stochastic model have to be done by identification. This part is not developed here but several tools are available in order to identified marginal densities and copula.

Once probabilistic model of eigenvalues is built, we are interested in representation of the joint density of random vector Re({H}).

In order to compare an empirical density of Re({H}) to the one obtained by our method, we run a Monte-Carlo simulations at a pulsation ω = 235 rad.s -1 from a sample of 5000 eigenvalues λ 1 , λ 2 drawn in p λ 1 ,λ 2 . Results of these simulations is presented by figure 7.

Density obtained from the expression (4) is presented by figure 8.

Influence of correlation of eigenvalues is illustrated by figure 9 which is a simulation of the joint PDF of Re({H}) if we consider independent eigenvalues. It appears clearly that correlation of eigenvalues strongly influence joint PDF of Re({H}). Moreover, empirical joint PDF of Re({H}) (figure 7) seems really close to simulated one (figure 8) which confirms capability of analytical formulation based on copula transformation.

As we can now express analytically joint PDF of Re({H}), knowing eigenvalues stochastic model, we use it to numerically compute the characteristic function (3).

Computation of the real part FRF CDF

For illustrating our methodology we consider [ϕ] bi = [ϕ] ai = 1 and then we compute the characteristic function (3). Numerical integration of dimension two (equation ( 3)) is performed on the domain This domain is meshed by quadrangles which are tensorial product of marginal p λ 1 , p λ 2 discretizations (according to differential of marginal).

[ 1- 4ξ 1 ω 2 (ξ 1 -1) , 1- 4ξ 1 ω 2 (ξ 1 +1) ] × [ 1- 4ξ 2 ω 2 (ξ 2 -1) , 1- 4ξ 2 ω 2 (ξ 2 +1) ].
In the following F ab stands for real part of random FRF between two DOF a and b and f ab for one sample. Once characteristic function is computed, several post treatments can be performed. First of all, Gil-Pelaez formula (Gil-Pelaez 1951) is used to access the CDF of F ab and this formula gives

P [F ab < f ab ] = 1 2 - 1 π ∞ 0 Im[e -jf ab v Φ Re(F ab ) (v)] v dv
In order to compute this integral, we use the method described by Davies [START_REF] Davies | Numerical inversion of a characteristic function[END_REF] based on Fourier cosine series summation formula.

P [F ab < f ab ] = 1 2 - K k=0
Im Φ Re(F ab ) ((k + 0.5)h)e -j(k+0.5)hf ab π(k + 0.5) .

The author recommends to choose a discretization step h, for which max P (F ab < f ab -2π h ) , P (F ab > f ab + 2π h ) is less than half the maximum allowable error. In our application F ab is bounded and

f ab ∈ [B 1 , B 2 ] with B 1 = [ϕ] a1 [ϕ] b1 4ξ 1 ω 2 (ξ 1 -1) + [ϕ] a2 [ϕ] b2 4ξ 2 ω 2 (ξ 2 -1) and B 2 = [ϕ] a1 [ϕ] b1 4ξ 1 ω 2 (ξ 1 +1) + [ϕ] a2 [ϕ] b2 4ξ 2 ω 2 (ξ 2 +1) . Then it follows that P (F ab > B 1 + 2π h ) = 0 ⇔ B 1 + 2π h > B 2 ⇔ h < 2π B 2 -B 1 .
As series converge fast to zero, truncation until the K term is done when Φ Re(F ab ) ((K + 0.5)h) = 0 + .

Applying this methodology to the previous example and still considering ω = 235 rad.s -1 , the figure 10 Gil-Pelaez formula considering or not dependence of eigenvalues. Some conclusions can be drawn from this figure. Firstly, eigenvalues dependence structure must be considered to avoid strong differences in CDF of FRF between two DOF. Secondly, comparison with Monte-Carlo allows to be confident in integration methodology described above. Finally, these conclusions are reinforced by providing the corresponding PDF. This one (figure 11) is computed by Poisson summation formula [START_REF] Abate | The Fourier-series method for inverting transforms of probability distributions[END_REF] which is in this case p F ab (f ab ) ≈ h π K k=0 Φ Re(F ab ) (kh) exp(-jkhf ab ).

-0.0008 -0.0006 -0.0004 -0.0002 0.0000 0.0002 0.0004 0.0006 0.0008 In previous example, ω is fixed to 235 rad.s -1 which is a value between the two oscillators deterministic natural pulsation. First section of this paper shows that marginal eigenvalues densities are strongly influence by value pulsation, ω. Then, we illustrate the methodology on a range of pulsation. Considering the same example and looking for 2.5%, 50%, 97.5% quantiles of F ab , previous methodology is applied and compared with Monte-Carlo method (100000 runs) on the pulsation range [180 -290]rad.s -1 . Figure 12 presents the results. As for previous results, comparison with Monte-Carlo method confirms capability of the two dimensional integration method. Moreover, it is interesting to see that this capability is conserved even at pulsation for which low variation of the FRF occurs. Influence of eigenvalues correlation clearly appears for 50% and 97.5% quantiles. For 97.5% quantiles, independent case reveals a single peak around ω = 230 rad.s -1 , which is not present in the dependent case.

All these numerical examples illustrates that correlation between eigenvalues strongly influence probability law of a FRF between two DOF. In addition, they show that copula is a suited tool to describe joint PDF of eigenvalues and then has to be used in determination of FRF CDF. Finally, integration method set-up in this context allows to consider that computation of the FRF CDF is performed with a negligible error.

This methodology is now applied to the satellite TARANIS (figure 3).

INDUSTRIAL APPLICATION : TARANIS

We consider the finite element satellite structure TARANIS previously introduced in part 2.2. More precisely, we investigate a particular FRF which deterministic value is drawn on figure 13.

This FRF is particularly interesting to apply our methodology as the two first modes of the structure are stimulated and have very close natural pulsation.

As it have been told in part 3.1, in industrial applications uncertainties affect input parameters of the model. Without detailing this modelisation, parameters concerned in this example are Young modulus, Poisson ratio, panel thickness etc. Then, eigenvalues densities and correlation must be identified. In this example, we use a nonparametric identification by kernel smoothing over 100 realisation to get p λ 1 and p λ 2 . It is assumed that dependence structure is normal. Figure 14 shows the Kendall plot between the one hundred sample and the normal copula used to model dependence (Kendall plot is equivalent to quantilesquantiles plot in two dimension, informations can be found in (Clauss ) or (OpenTURNS 2011)). This curve shows that normal copula seems to be correct to model dependence between the two first eigenvalues.

Once identification process is done, one can applies the previous methodology. As above, we compute the 2.5%, 50%, 97.5% quantiles of the real part FRF. In order to have reference results, a Monte-Carlo simulation of 1000 runs is also performed and used to compute the 50% quantile. Note that, in this case, Monte-Carlo methodology consists in solving 1000 eigenvalues problems and compute corresponding FRF. quantiles shapes seems relevant but we cannot confirm these results with only 1000 Monte-Carlo runs. Finally, this industrial example shows that identification of eigenvalues joint PDF is realistic and leads to relevant results.

CONCLUSION

Aims of this study was to present a method of probabilistic modal synthesis based on analytical expressions of the random oscillator PDF. All this article focuses on the importance of correlation between eigenvalues. First of all, modeling of eigenvalues dependence structure is performed by copula. Then we show how to use it in the determination of the joint PDF of the random vector of frequency response function Re({H}). Next, this joint PDF is used to determine the characteristic function of a frequency response function between two DOF and several solutions to numerically compute and make use of it are presented. At each step of this methodology, effect of correlation between eigenvalues is illustrated through a two dimensional example in part 3. These illustrations allow to draw two conclusions. First, the effect of correlation is non negligible and strongly affect the PDF of a FRF between two DOF of the system. Second, comparison with Monte-Carlo simulations confirms that integration methodology is rel-evant and leads to a really good approximation of the PDF. Finally, an industrial application on Taranis satellite shows that eigenvalues joint PDF identification is realistic and allows to access some quantiles of the PDF FRF. Moreover, as it has been mentioned in introduction, this identification may be parametric, like in example part 3, or nonparametric, like in Taranis example (kernel smoothing). This is an advantage when only few information are available on eigenvalues variations (for example for some equipments housed in the satellite).

Perspectives for this work are to continue analytical process in the exploitation of copula properties (determination of frequency boundaries in which correlation have to be modeled). Then, we shall apply this methodology in a sub-structured model for which eigenvalues joint PDF of some sub-structures comes from nonparametric random matrices model such as the ones proposed by [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF] or [START_REF] Adhikary | Eigenvalue density of linear stochastic dynamical systems: Arandom matrix approach[END_REF].

FunctionpFigure 1 :

 1 Figure 1: Shape of the density p Re(H) (Re(h)) when p λ follows a normal law.

  ai . Uncertainties on input parameters leads to a random eigenvalues problem. At this step, we assume that the random eigenvector basis is close to the deterministic basis [ϕ]. Random transfer matrix [H] is diagonal and each term is a random oscillator studied in first part. Let us call Re({H}) the random vector composed by real part of diagonal terms of [H]. As transformation giving FRF between two DOF a and b of a structure is a linear combination between terms of the random vector Re({H}), one obtains the characteristic function of Re(F ab ).

Figure 3 :

 3 Figure 3: Finite element model of TARANIS satellite.

2

  Dependence between λ 1 and λ 2

Figure 4 :

 4 Figure 4: Two first eigenvalues of TARANIS model, 1000 runs.

Figure 5 :2

 5 Figure 5: Marginal densities p λ1 , p λ1 .

Figure 6 :

 6 Figure 6: Random simulation of λ 1 , λ 2 following joint density p λ1,λ2 .

Figure 7 :Figure 8 :

 78 Figure 7: Random simulations of Re({H}), at ω 235 rad.s -1 .

Figure 9 :

 9 Figure 9: Simulation of the density of Re({H}), at ω = 235rad.s -1 , with independent eigenvalues.

Figure 10 :

 10 Figure 10: Comparison between CDF of F ab obtained by Monte-Carlo method (red cross) and Gil-Pelaez formula, taking into account correlation of eigenvalues (blue curve) or not (green curve), at ω = 235 rad.s -1 .

Figure 11 :

 11 Figure 11: Comparison between PDF of F ab obtained by Monte-Carlo method (red histogram) and Gil-Pelaez formula, taking into account correlation of eigenvalues (blue curve) or not (green curve), at ω = 235 rad.s -1 .

Figure 13 :

 13 Figure13: Deterministic FRF study in TARANIS example.

Figure 14 :

 14 Figure 14: Kendall plot of the two first eigenvalues fitted by a normal copula.

  Determination of the 2.5%, 50% and 97.5% quantiles by integration and comparison with Monte-Carlo simulations, p λ follows a normal law.

		0.0004		2, 5%, 50% et 97, 5%, Quantiles Comparison with Monte-Carlo simulation, 1000000 runs
											Integration
		0.0003									Monte-Carlo
		0.0002
	Re(H(ω))	0.0000 0.0001
		-0.0001
		-0.0002
		-0.0003	0	50	100	150	200	ω	250	300	350	400	450
	Figure 2:

  illustrates CDF of F ab . Comparison is done between Monte-Carlo method (one million runs), and

		1.2	CDF, ω = 325rad.s -1 .
			Monte-Carlo
		1.0	Dependent case Independent case
		0.8
	ab < f ab )	0.6
	(F	0.4
	P
		0.2
		0.0
		-0.0008 -0.0006 -0.0004 -0.0002 -0.2	f ab 0.0000	0.0002	0.0004	0.0006	0.0008

ACKNOWLEDGMENTS

This work is supported by Centre National d Études Spatiales (CNES) and Thales Alenia Space (TAS). The authors thank Fabrice Buffe (CNES), Jérôme Buffe and Anne Carpine (TAS).


Thibault Gouache

Joseph Morlier

Guilhem Michon

Baptiste Coulange

OPERATIONAL MODAL ANALYSIS WITH NON STATIONNARY INPUTS

Keywords: Operational, Modal, Analysis, Transient, Phase, Pressure Oscillations, Solid rocket motors

Operational modal analysis (OMA) techniques enable the use of in-situ and uncontrolled vibrations to be used to lead modal analysis of structures. In reality operational vibrations are a combination of numerous excitations sources that are much more complex than a random white noise or a harmonic. Numerous OMA techniques exist like SSI, NExT, FDD and BSS. All these methods are based on the fundamental hypothesis that the input or force applied to the structure to be analyzed is a stationary white noise. For some applications this hypothesis is reasonable. However in numerous situations, the analyzed structure is subject to harmonic and transient forces. Numerous methods and research has enabled to develop methods that are robust to such harmonic contributions. To enable OMA during pressure oscillations in solid rocket boosters, the authors propose to consider transient and harmonic inputs no longer as parasites but as the main force applied to the structure that must be analyzed. This is the case during pressure oscillations in rocket boosters. We propose the use of phase analysis adapted to a transient context to conduct operational modal analysis under a harmonic transient input. This time-based novel OMA method will be exposed. The theoretical developments and algorithmic implementations are exposed. First tests have been conducted on laboratory single degree of freedom setup to validate this new OMA technique and are reported here.

INTRODUCTION 1.Classical operational modal analysis

Most Operational Modal Analysis (OMA) techniques are based on a fundamental hypothesis: the unknown and unmeasured input of the studied system or structure is stationary white Gaussian noise.

The Natural Excitation Techniques method, which uses auto and cross correlation is based on this hypothesis [START_REF] James | The Natural Excitation Technique (NExT) for Modal Parameter Extraction from Operating Structures[END_REF]. The Frequency Domain Decomposition method, based on singular value decomposition is also built on the hypothesis of spatially and harmonically white input [START_REF] Brincker | Modal Identification from Ambient Responses using Frequency Domain Decomposition[END_REF]. Stochastic subspace identification methods (using covariance or subspace technique) [START_REF] Peeters | Stochastic Identification for Operational Modal Analysis : A Review[END_REF] are another example of such techniques. These techniques have been applied to a various number of structures: bridges [START_REF] Luscher | Modal Parameter Extraction of Z24 Bridge Data[END_REF], aircrafts in flight [START_REF] Hermans | Modal Testing and Analysis of Structures Under Operational Conditions : Industrial Applications[END_REF], ship structures [START_REF] Rosenow | Investigation of the Dynamic Behaviour of Ship Structures Using Classical and Operational Modal Testing[END_REF], wind turbines [START_REF] Moller | Stochastic Subspace Identification Technique in Operational Modal Analysis[END_REF], commercial satellites and many other structures [START_REF] Peeters | Industrial Relevance of Operational Modal Analysis -Civil Aerospace and Automotive Case Histories[END_REF].

But, in addition to white noise, the input can contain some dominant frequency components (due to a running engine for example) [START_REF] Peeters | Stochastic Identification for Operational Modal Analysis : A Review[END_REF]. The input of the system can also evolve in time and induce transient effects [START_REF] Hermans | Modal Testing and Analysis of Structures Under Operational Conditions : Industrial Applications[END_REF]. To take into account the limit imposed to OMA techniques by the fundamental stationary white noise input assumption, two general approaches are used. The first one consists in developing methods that are robust and insensible to the violations of this hypothesis. In [START_REF] Cauberghe | Identification of modal parameters including unmeasured forces and transient effects[END_REF], B Cauberghe et al. studied the robustness of their frequency-domain OMA technique against transient effects. The second approach used is to enrich or to relax the fundamental hypothesis of white noise. The OMA methods based on the cesptrum domain relax the classical OMA hypothesis: they only suppose that input is cepstrum short [START_REF] Randall | Extraction Of Modal Parameters From The Response Power Cepstrum[END_REF]. It is also possible to suppose that the structure input is a superposition of Gaussian white noise and a harmonic excitation. In [START_REF] Pintelon | Continuous-time operational modal analysis in the presence of harmonic disturbances[END_REF], R. Pintelon et al. proposed three methods to suppress unknown harmonic time varying pulsations. Generally, enriching the fundamental hypothesis of OMA gives birth to new OMA techniques like those proposed by P. Mohanty and D.J. Rixen [START_REF] Mohanty | Operational modal analysis in the presence of harmonic excitation[END_REF][START_REF] Mohanty | Modified SSTD method to account for harmonic excitations during operational modal analysis[END_REF][START_REF] Mohanty | A modified Ibrahim time domain algorithm for operational modal analysis including harmonic excitation[END_REF][START_REF] Mohanty | Modified ERA method for operational modal analysis in the presence of harmonic excitation[END_REF].

Transient vibrations (solid rocket motors and many other cases)

A large number of transient phenomenon exist and this should be taken into account in OMA techniques [5; 11]. Transient pressure oscillations are often observed in Solid Rocket Motors (SRM) [START_REF] Chedevergne | Detailed Analysis of the Thrust Oscillations in Reduced Scale Solid Rocket Motors[END_REF]. Everyday structures are also subject to transient inputs: plane fluttering [START_REF] Zouari | Détection Précoce de l'Instabilitée Aéroélastique des Structures Aéronautiques[END_REF], wind gusts on bridges, waves on dams and, earthquakes on buildings are some examples. Instead of developing robust OMA algorithms against transient effects like in [START_REF] Cauberghe | Identification of modal parameters including unmeasured forces and transient effects[END_REF], it is proposed here to take into account the presence of transient harmonic excitations in a structure's input to develop a new OMA technique based on transient inputs.

Paper layout

After describing the properties of resonance, the first section of this paper proposes the novel definition of "Transient Phase" (TP). To better understand the new notion of transient phase, a Single Degree Of Freedom (SDOF) transient oscillator model is introduced. The second section is dedicated to the study of transient phase in the SDOF case: a mathematical parametric reduction proves that transient phase only depends on the ratios of four characteristic time quantities and some typical transient phase profiles are exposed. This allows the authors to propose a SDOF OMA technique based on the existence of transient regimes. Section three presents the SDOF experimental setup, the experimental validation of mathematical parametric reduction and the validation of the proposed SDOF OMA technique. Finally, to accommodate this new SDOF OMA technique to multiple degree of freedom structures, an algorithm is proposed based on mixed time-frequency decomposition. This algorithm and the notion of transient phase are then tested on Solid Rocket Motor (SRM) firing tests.

RESONANCE AND PHASE IN A TRANSIENT CONTEXTE

Definition of transient phase

Resonance is identified thanks to the maximum dynamic amplification or thanks to a phase difference between input and output equal to Π/2. In the case of SRM pressure oscillations, it is very difficult to measure the forces applied to the SRM's structure. It is thus impossible to measure phase difference or dynamic amplification between output and input. With SRM pressure oscillations, a second difficulty is the non-stationary nature of the forces applied to the SRM's structure. However the transient nature of the input can be used to detect resonance. Indeed, before a system is subject to resonance it will vibrate at a given phase. But as the forced harmonic regime is progressively imposed to the mechanical system and has the established regime is progressively reached, the phase of the system will change. Though no absolute value of phase (compared to the input phase) can be defined without access to the input, a measure of the evolution of the phase output can be observed. The evolution of output phase will be directly linked to the mechanical properties of the system and its input.

The transient response of a single degree of freedom system will have a variable amplitude. To adapt the notion of phase it is necessary to take this amplitude variation into account. The local phase of the system will thus be defined based on the local amplitude of the response. This is illustrated in Figure1. The transient phase phase, is thus defined by the equation here under where s(t) is the value of the response and local_max is the local maximum of the half-period as illustrated in Figure 1.

!ℎ!"#(!) = asin !(!) !"#$!_!"# (!) (1)
Finally, since phase is always defined based on a reference, and that the input is not accessible it is necessary to define a reference. In this paper, the reference will be a pure harmonic signal whose frequency is equal to the frequency that maximizes s(t) Fourier transform.

Single degree of freedom transient model

To illustrate and test the notion of transient phase, the following transient single degree of freedom model is adopted (presented in Equation 2). x is the position; the dot is the time derivative; t is time; m, c and k are respectively system mass, damping and stiffness; F is the harmonic amplitude; ω is the harmonic input's frequency (in rad/s) and α(t) is a function that enables input amplitude variation.

!! + !! + !" = ! • !(!) • cos !" (2)
In this paper the function α(t) is defined by Equation 3, where α is a coefficient enabling amplitude variation speed modulation. If α is close to 0, amplitude build up is very slow, and if α is very high amplitude build up is very fast.

! ! = 0 ! ≤ 0 1 -! !!" ! > 0 (3)
By dividing Equation 2 by m and by using linear variable changes on t and x, it is possible to show that Equation 2 depends only on the following three parameters shown in Equation 4.

! !! ! ! ! ! ! ! ! ( 4
)
Since phase is defined only by amplitude ratios, it is unchanged by linear variable changes on x and since phase is defined based on a reference it is also unchanged by linear variable changes on t. The transient phase profiles will be function only of the three parameters described in Equation 4.

Transient phase profiles

To validate the conclusions of previous section, numerical simulations of a transient single degree of freedom oscillator were conducted and the transient phase profiles were calculated based on the output of these simulations. the simulations were conducted thanks to SCILAB freeware developed by INRIA. These simulations confirmed the fact that transient phase depends only of three variables described in Equation 4.

Transient phase operational modal analysis

Based on the observations in sections here above, the authors developed an Operational Modal Analysis technique based on transient phase profiles. Indeed, to identify the parameters of Equation 2 and 3, it is possible to extract from operational data the transient phase of the transient response. By using curve fit algorithms it is then possible to identify the three parameters of Equation 4 directly on the transient phase profile of the systems response. Then by using the reference frequency and the real amplitude of the response it is possible to identify ω and F/m.

The limitation of such a modal identification technique is that a correct parametric model of the input must be used. Indeed to identify the three parameters on the transient phase profile, transient profiles are simulated and feed into a curve fit algorithm to find the best fit to the experimental transient phase.

However the model does not need to be highly precise. Indeed one can substitute the amplitude function α(t) defined in Equation 3 by a linear ramp or by other monotonous, continuous increasing functions without modifying significantly the transient phase profiles. The strength of this approach is very generic and that any model can be used to simulate the behavior of the system. If a mechanical system that must be identified is subject to specific inputs (time varying frequency) or has specific properties (non-linearity) it is possible to take these into account in the model.

EXPERIMENTAL VALIDATION

Experimental Setup

To validate the theoretical developments in sections here above and to test the single degree of freedom transient phase operational modal analysis technique proposed, a simple damped-mass-spring experimental setup was used. The experimental setup is composed of a spring attached to a mass and is described in Figure 3. The mass is guided in translation. Four different springs and three masses are used. The mass is connected to a viscous damper. The viscous damper is made of two parallel plates that bath in a cylinder full of oil.

Experimental transient phase profiles dependency on three parameters

First experiments were conducted to verify that experimental transient phase profiles are function exclusively of the three parameters in Equation 4. Figure 4 presents two transient phase profiles (one in red the other in blue), that were obtained experimentally. They were obtained with different masses, spring constants, damping, input frequency, amplitude and transient amplitude variation speed. However the six parameters were chosen in order for the ratios in Equation 4to be all equal. The experimental transient phase profiles are very close to each other (correlation of 0.98). Numerous other tests were conducted and gave similar results. This confirmed experimentally the dependency of transient phase profiles exclusively on the three parameters present in Equation 4.

Correlation between experimental and simulated transient phase profiles

The experimental transient phase profiles were then compared to the numerical simulated ones. A genetic curve fit algorithm was used to find the best fitting simulated phase profile. The algorithm maximizes correlation between the experimental and simulated transient phase profile. Figure 5 shows an example of an experimental transient phase profile and the best simulated fit. The agreement between numerical and experimental transient profile is very good (correlation of 0.91).

Operational modal analysis identification results

Numerous identifications were conducted on the experimental transient data. Table 1 presents the maximum error on identification. The identification of the input and modal frequency are very precise and present highly satisfactory results. The identification of the damping can present significant errors but when compared t typical errors of other operational modal analysis methods, they are reasonable. Errors on α are much more important than errors on the other parameters. However this value is generally not identified in other operational modal analysis algorithms.

These first results are encouraging. It is possible to increase the quality of these results by proposing more efficient curve fit algorithms. These results were generated by a genetic algorithm that maximizes the correlation between the experimental transient phase and the simulated one. It is possible to optimize other measures between the two transient phases or to implement other optimization algorithms.

CONCLUSIONS

Some structures like solid rocket motors used on rockets to launch satellites are subject to violent transient regimes characterized by harmonic inputs with rapidly varying amplitudes. In order to better understand these phenomenon and to better identify the modal properties of these structures, it is necessary to develop operational modal analysis techniques that can cope with transient inputs. This paper has introduced the basic building blocks and experimental and numerical explorations conducted to develop such a technique.

The notion of transient phase was introduced. Using a single degree of freedom model with a transient input, it was shown that transient phase depends only on three parameters rather than six (as does the full time response of such a system). This enables authors to propose an operational modal technique enabling the identification of a structures properties during transient inputs. The first step is the identification on the experimental transient phase of three basic parameters that define the transient phase. The other parameters are then identified by using the original time response of the system. This technique requires modeling the system's input. However the model does not need to be too precise.

Experiments were conducted on mass-spring-damper setup to confirm the theoretical developments and to test the identification method. These experiments confirmed the dependency of transient phase on only three parameters. The experimental transient phase showed close agreement to the simulated ones. despite a very crude curve fit algorithm, the first modal identifications conducted with this novel method showed promising performances.

Based on this first building blocks, the modal identification of data coming from solid rocket boosters has started. The multi-degree of freedom system (that a rocket booster is) is considered as multiple independent single degree of freedom systems. Thus by filtering in the frequency domain, it is possible to use the single degree of freedom model proposed here. First results are promising and will be presented in future publications.

Figure 1

 1 Figure 1 Illustration of local maximum definition with variable amplitude response of mechanical system.

Figure 2

 2 Figure 2 present examples of transient phase profiles obtained through numerical simulations. The transient phase resemble free responses. When the ratio ω/ω0 is above one the transient phase tends to decrease, when it is smaller than one the transient phase decreases. The number of oscillation depends on this ratio and the damping. Here the same damping was used all 4 simulations. Had a smaller damping been simulated, more oscillations would have been observed.

Figure 2

 2 Figure 2Examples of transient phase profiles obtained thanks to the numerical simulations of a single degree of freedom oscillator. The ratio ω/ω 0 for each simulation is equal to 1.01, 1.05, 0.99, 0.95. For these figures damping is equal to 5%.

Figure 3

 3 Figure 3 Schematic of experimental setup. A spring-mass-oil damper system is activated by an electro-dynamic shaker. Acceleration and displacement are measured.

Figure 4

 4 Figure 4 Transient phase profiles (one in red and one in blue) versus time of two different experiments done with m, c, k and input parameters all different but the three ratios of Equation 4 all equal.

Figure 5

 5 Figure 5 Transient phase profiles (experimental one in red and simulated one in blue) versus time. Correlation between experimental and numerical transient phase is equal to 0.91.

Table 1

 1 Maximum error on identification

	Value	Mean dispersion
	ω	0.5 %
	ω 0	1.2 %
	damping	75 %
	α	180 %

Dr, Cornis SAS, thibault.gouache@cornis.fr

Dr, Université de Toulouse/ICA/ISAE, joseph.morlier@isae.fr

Dr, Université de Toulouse/ICA/ISAE, guilhem.michon@isae.fr

Dr, Cornis SAS, baptiste.coulange@cornis.fr

ACKNOWLEDGEMENTS

The authors would like to thank Prof Yves Gourniat from ISAE for allowing us to explore transient harmonic inputs. We would also like to thank Elsa Piollet and Claire Mauris-Demourioux for their help doing the experiments and their useful feed-back and ideas. finally we would like to thank Benoit Geffroy from the European Space Agency and Luc Gonidou from Centre National d'Etudes Spatiales (French space agency) for their helpful ideas and knowledge of pressure oscillations in solid motor rocket boosters.


L Bennani

P Villedieu

M Salaun

Two Dimensional Model of an Electro-Thermal Ice Protection System

In this communication we shall focus on the main governing equations and building blocks of the M.A.D (Anti-icing Deicing Modelling) numerical tool, which is now renamed as INUIT (Integrated NUmerical model of Ice protection sysTems) and part of the new generation of ONERA icing codes. The code simulates the functioning of an electro-thermal de-icing system. We shall also discuss the various improvements and new features we have added, especially a mechanical model of the ice block in order to improve the ice-shedding criterion.

I. Introduction

In the aeronautical world, icing is one of the most serious hazards that can be encountered. Icing is caused by freezing upon impact of supercooled water droplets. Not only does it increase mass but it may also lead to a degradation of aerodynamic performances, blocked air intakes and in the worst case, loss of control of the aircraft. Aircraft manufacturers must therefore comply with certifications and regulations regarding flight safety in icing conditions. In order to achieve that goal several ice protection technologies may be adopted. A popular system is the bleed air anti-icing system, where hot air is taken from the engine to heat the desired protected surfaces so as to prevent the formation of ice. However nowadays aircraft manufacturers are moving towards more economic, green and electrical aircraft. In this context a possible candidate is the electro-thermal ice protection system, which relies on cyclic heater mat activation to detach the ice.

This calls for numerical tools in order to assist the manufacturer during the design phase of these ice protection systems so as to assess their performances and suitability. Moreover it is most important to be able to predict the ice shape and the performance degradation it causes. In this context, a new generation of icing codes are being developed at ONERA. The ONERA 2D icing codes are currently being updated to more modern numerical methods and programming architectures. The work presented here concerns the updating and development of new models in the electro-thermal ice protection system simulation module. Figure 1 gives a schematic view of the new code's intended architecture: In the following sections, we will first describe the way of approaching the different heat transfer problems. Then we will describe how we intend to tackle the ice shedding problem.

II. The electro-thermal system

The system we are trying to model is composed of heater mats installed within a multi-layered material and can be used in anti-icing or deicing configurations. [START_REF] Henry | Development of an electrothermal de-icing/anti-icing model[END_REF][START_REF] Keith | Numerical simulation of icing, de-icing and shedding[END_REF] These mats are activated according to a predefined time dependent cycle in order to melt some regions of the ice at the ice/structure interface. This will reduce the ability of the ice to adhere to the surface and will lead to shedding. At the time being the system is modeled in two dimensions. Figure 2 gives an illustration of the system. As the structure containing the system remains fixed, only the ice block's frontier will be subject to change in time. Thus considering each domain has its own geometrical behaviour in time, it seems reasonable to treat them separately. What we mean by "separately" is that the meshes will be generated differently and separate modular blocks will be assigned in the code to treat each domain. However the physics will be treated simultaneously by a coupling procedure at the ice/structure interface.

III. Heat transfer in the (multi-layered) structure

III.A. Physical process and modelling

Here we must solve the problem of heat conduction, governed by the heat equation (eq 1), in a structure composed of different materials and with localized heat sources provided by the heating mats.

ρc p ∂T ∂t + ∇. -Λ --→ ∇T = S (1)
Thus on the one hand changes in material characteristic constants such as density and specific heat must be handled. On the other hand the heating mats are not explicitly represented as a part of the calculation domain and their presence must be modelled. Thermal conductivity is allowed to be anisotropic, having different conductivities along ξ and η (where (ξ,η) is a curvilinear coordinate system as illustrated figure 2). Therefore the thermal conductivity matrix Λ has the form:

λ ξ 0 0 λ η
We then proceed to a discretization of equation 1 that we solve with a finite volume method. We choose an implicit time marching algorithm due to the fact that some materials in the structure may be very good conductors. Indeed, this would put a harsh limit to the authorised value of the time step had we chosen an explicit formulation. The heating mats are considered to be located at the interface between two layers and are not explicitly meshed. In order to take into account their presence the product ρc p of each adjacent cell is modified according to the following formula, defining an equivalent ρc p : 9 (ρc p ) eq = (ρc p ) j e j + (ρc p ) r e r /2

e j + e r /2 (2)
Where e j and e r are the thickness of the adjacent cell and the heating element respectively. If P is the power of the heating element, the heat provided is represented as a heat source S j = P 2 in each adjacent cell. We assume perfect thermal contact between layers, but this assumption could easily be extended to non perfect contact in a future version of the code. Figure 3 shows an instantaneous shot of the temperature field in the structure taken during a de-icing cycle. We can clearly see the zones where the heater pads are active at that instant.

III.B. Academic validation case

As a first step of validation we consider a simple geometry, a hollow cylinder, made of a single material, which has a fixed uniformly distributed temperature T o on its inner and outer surfaces and a uniform initial temperature field T i (see figure 4). The analytic solution, given in, [START_REF] Vansant | Conduction heat transfer solutions[END_REF] is a temperature field that depends only on the radial coordinate r and time.

T -T o T i -T o = π +∞ n=1 f (R, λ n )e -λn αt ro (3a) f (R, λ n ) = J 0 (λ n R i ) [J 0 (λ n R)Y 0 (λ n ) -J 0 (λ n )Y 0 (λ n R)] J 0 (λ n R) + J 0 (λ n ) (3b) J 0 (λ n R i )Y 0 (λ n ) -J 0 (λ n )Y 0 (λ n R i ) = 0 (3c) R = r r o (3d)
We compare this solution to the solution obtained by the heat transfer module at different times. As shown on figure 5 we have a very good match between the analytic and computed profiles. IV. Heat transfer in the ice block

IV.A. Physical process and modelling

In the ice block the material remains the same over the whole domain, that is to say, water. The specific physical process that arises when dealing with this domain is phase change. This mechanism involves a discontinuous enthalpy-temperature relation. In order to capture this discontinuity many methods exist (eg : level set, phase field). In our case the so called enthalpy method is chosen. This method enables an efficient capturing of the melting front. [START_REF] Alexiades | Mathematical modelling of melting and freezing processes[END_REF] The main idea is to return to the general energy conservation equation :

t+∆t t ∂ ∂t V ρhdV dt = t+∆t t ∂V -- → q • -→ n dSdt (4)
The heat flux -→ q is computed using the temperature field. Knowing the flux then enables us to update the value of the enthalpy. The following relationship (5) then gives us the new values of temperature :

T =      T m + h cs h ≤ 0 T m 0 ≤ h ≤ L T m + h-L c l h ≥ L (5)
and we can update the liquid fraction, which is a crucial element for computing the heat fluxes :

φ L =      0 h ≤ 0 h L 0 ≤ h ≤ L 1 h ≥ L (6)
Here the following convention has been chosen for h :

h =      c s (T -T m ) T < T m solid φ L L T = T m mixed state c l (T -T m ) + L T < T m liquid (7)
In this approach the density ρ is assumed constant (and the same for water in both solid and liquid states). The equations are solved using an explicit time stepping finite volume method. In this case the explicit formulation is justified by the fact that ice and water are not very good heat conductors. Therefore the time step is not too limited by the stability condition. The open source meshing software Triangle is used to generate the unstructured mesh. [START_REF] Shewchuk | Triangle: Engineering a 2d quality mesh generator and delaunay triangulator[END_REF] Figure 6 shows preliminary results of this module. An instantaneous shot of a phase indicator (liquid fraction) field is represented. We clearly see the formation of a liquid water film at the heated surface.

IV.B. Academic validation case

As in the case of the heat transfer in the structure (section III.B) we would first like to know if our numerical procedure correctly solves the mathematical equations used to model our problem. In this case we are going to consider the case of a rectangular slab of ice uniformly heated at one of its surfaces. The other boundaries are assumed adiabatic. The case is illustrated fig 7. The problem reduces itself to a one dimensional one, for which analytic solutions are available. [START_REF] Alexiades | Mathematical modelling of melting and freezing processes[END_REF] In this case the position of the melting front is given by:

x Γ (t) = 2χ √ α l t ( 8
)
Where χ is the solution of the transcendental equation :

St l e χ 2 erf (χ) - St s νe ν 2 χ 2 erf c(νχ) = χ √ π (9)
With

St l = c l (T L -T m ) L St s = c s (T m -T (x, 0)) L ν = α l α s (10)
In order to compare the analytic and numerical results we perform an integration of the numerical liquid fraction over the whole domain. The result is then divided by the thickness of the slab. This gives us the position x Γ of the melting front. The time dependent positions of these fronts are shown figure 8. We observe a good match between analytic and numerical results, although the two curves tend to very slightly drift apart with time.

V. Mechanical models and ice shedding

Once we have determined the regions of the ice block which have melted and the global temperature field, we must predict if ice shedding occurs or not. To do so an empirical model was originally used in the code, as we will see in section V.A. Nevertheless such a criterion is unsatisfactory since it does not take into account aerodynamic and inertial forces (ex: for rotating blades).

Modelling techniques based on continuum mechanics have already been applied to ice. Scavuzzo et al. performed a finite element analysis of the stress distribution due to aerodynamic forces in an accreted ice block. [START_REF] Scavuzzo | Influence of aerodynamic forces in ice shedding[END_REF] More recently Zhang et al. have used a crack propagation and re-meshing technique to study ice break up. [START_REF] Khurram | Fem analysis of in-flight ice break-up[END_REF] As for adhesion, in a general way, there are several identified mechanisms (mechanical interlocking, electrostatic attraction, chemical, adhesion by diffusion [START_REF] Shanahan | Adhesion and wetting: similarities and differences[END_REF] ). The problem is then to identify by which of these mechanisms ice adheres to a surface. Some theoretical models have been proposed to explain and predict ice adhesion. One assumes the existence of a liquid like layer on ice. [START_REF] Jellinek | Ice adhesion[END_REF][START_REF] Jellinek | Liquid-like (transition) layer on ice[END_REF] Another is a model based on the electrostatic adhesion mechanism. [START_REF] Petrenko | Physical mechanisms responsible for ice adhesion[END_REF] In addition, Fortin and Perron have recently developed an ice adhesion model based on electrostatic and mechanical interlocking mechanisms. [START_REF] Perron | Ice adhesion models to predict shear stress at shedding[END_REF] Although the adhesion models cited previously are based on more fundamental physical considerations, in our work we adopt a more global approach based on continuum damage mechanics.

One of the main problems that arises now is to determine what mechanical properties we are going to use in order to characterise atmospheric ice, for unfortunately very few studies on the subject exist. Most studies are interested in the tensile or compressive strength but do not provide many information on mechanical characteristics. [START_REF] Farzaneh | An experimental study on the tensile properties of atmospheric ice[END_REF][START_REF] Laforte | The adhesion of glaze and rime on aluminium electrical conductors[END_REF][START_REF] Farzaneh | Bending strength and effective modulus of atmospheric ice[END_REF] Moreover these experiments are very difficult to conduct due to the vast number of parameters on which those properties depend. Eskandarian 7 reports experimental measurements of Young's modulus and Poisson's ratio for porous ice. These results could be used for rime ice, in which small air bubbles are trapped, but this would require us to know the actual value of the porosity. Thus we see that it is very complicated to correctly describe the mechanical properties of atmospheric ice.

Therefore as a first approach we will use the values given by experiments for natural ice, which are more precise, and more widely studied. According to Schulson and Duval 6 the most precise results were obtained by Gammon et al. [START_REF] Clouter | Elastic constants of artificial and natural ice samples by brillouin spectroscopy[END_REF] In the case of a homogeneous polycrystalline ice aggregate, it may be assumed that the grains (crystallites) are oriented randomly making it elastically isotropic. The parameters characterising homogeneous isotropic elastic behaviour of polycrystalline ice Ih are shown in the following table ( The temperature dependence of these constants is obtained using the formulas given in, [START_REF] Duval | Creep and fracture of ice[END_REF] however in this first approach we will neglect the effects of temperature on mechanical properties.

As for the fracture toughness, according to the data presented in, [START_REF] Duval | Creep and fracture of ice[END_REF] we choose a value of 100kP a/m -1/2 .

V.A. Empirical model

The ice adhesion model that was originally implemented in M.A.D is empirical. It states that if the ratio Lwf /Lt between water film length and total ice block length becomes greater than a certain user defined limit value (80% recommended) then the ice block detaches from the structure. A variant of this criterion based on the water film height is also implemented. These criteria illustrate the fact that ice adhesion is a very delicate mechanisms to model.

V.B. Linear elastic damage model

V.B.1. Presentation of the model

Brittle fracture is one of the most striking phenomena involved in ice shedding. Several approaches to fracture mechanics exist such as XFEM or remeshing. Here we chose to use a method that resembles those of continuum damage mechanics.

In the following paragraphs we are going to re-derive the equations obtained by Miehe et al 5 following an auxiliary (although closely related) method. The starting point is a principle of conservation of energy whose physical motivation is the following: when a solid is deformed by action of external forces it internally stores elastic deformation energy. If, locally, this energy exceeds a certain critical energy then it will go into creating an increase in crack surface. Therefore, if a crack nucleates and/or propagates, a possible mechanism (from a macroscopic point of view) could be a process of energy transfer from the applied forces to elastic deformation energy to crack surface energy.

The idea is to introduce a parameter d that characterises the local state of damage/fracture 12 in the ice block. Let's start by first stating the general equilibrium equations of continuum mechanics. Consider a domain Ω on which is applied a body force f vol and a surface load f surf on Γ 1 . On Γ 2 the displacement field is imposed to u d . Then we have :

-div(σ(u, d)) = f vol in Ω σ.n = f surf on Γ 1 u = u d on Γ 2 (11)
Let ψ(ǫ, d) be the elastic deformation energy per unit volume, and φ(d, ∇d) be the crack energy per unit volume. Therefore the total elastic energy will be E el = Ω ψ(ǫ, d)dV and the total crack energy will be E crack = Ω φ(d, ∇d)dV . In order to evaluate the energy dissipated by a crack Bourdin et al 3 have proposed a regularization of the crack surface density, given by the functional 12. According to Miehe et al [START_REF] Hofacker | A phase field model for rate-independant crack propagation robust algorithmic implementation based on operator splits[END_REF] this is analogous to a phase-field model where d would be the order parameter.

γ(d, ∇d) = 1 2l d 2 + l 2 ∇d.∇d (12)
The length scale l is a spreading parameter for the crack. The energy of fracture is obtained by multiplying this surface by g c , the energy release rate, and is given by ( 13)

E crack = Ω g c [ 1 2l d 2 + l 2 ∇d.∇d]dV (13)
When external forces are applied, the change in energy is equal to the work produced by those forces. Thus we have ( 14)

δE crack + δE el = Ω f vol δudV + Γ1 f surf δudΓ (14)
Using relation (13) and integrating by parts we find :

δE crack = Ω g c l d -g c l△d δddV + ∂Ω g c l∇d.nδddΓ (15)
Using 15 and δE el = Ω ∂ψ ∂ǫ δǫ + ∂ψ ∂d δd dV we find:

Ω ∂ψ ∂ǫ δǫ + ∂ψ ∂d + g c l d -g c l△d δddV + ∂Ω g c l∇d.nδddΓ = Ω f vol δudV + Γ1 f surf δudΓ
Moreover, the mechanical boundary condition gives us :

∂Ω f surf δudΓ = ∂Ω σ.nδudΓ = Ω div(σ)δudV + Ω σ : ∇δudV (16)
As the Cauchy stress tensor is symmetric, we also have σ : ∇δu = σ : δǫ. Therefore we get:

Ω f vol δudV + ∂Ω f surf δudΓ = Ω [div(σ) + f vol ] δudV + Ω σ : δǫdV
Given the mechanical equilibrium relation div(σ) + f vol = 0, relation (14) reduces to :

Ω ∂ψ ∂ǫ -σ δǫdV + Ω ∂ψ ∂d + g c l d -g c l△d δddV + ∂Ω g c l∇d.nδddΓ = 0
Thus we finally obtain (17) :

σ = ∂ψ ∂ǫ in Ω -∂ψ ∂d = gc l d -g c l△d in Ω ∇d.n = 0 on ∂Ω (17)
We see that the stress tensor σ is deduced from the choice of ψ. In the case of linear isotropic homogeneous elasticity it is well established that ψ takes the form:

ψ(ǫ) = λ 2 tr(ǫ) 2 + µtr(ǫ 2 ) (18)
However, considering -∂ψ ∂d is the source term, the driving force of the fracture process, we must also chose ψ according to what part of the elastic deformation energy creates or propagates a crack. It is considered that a crack may nucleate or propagate only under local tension. Therefore ψ is decomposed into the sum of a tensile and a compressive part, which are defined by analogy with 18 and using the eigen values of ǫ as follows:

ψ +/- 0 (ǫ) = λ 2 < ǫ 1 + ǫ 2 > 2 +/-+µ(< ǫ 1 > 2 +/-+ < ǫ 2 > 2 +/-)
As fracture may only be influenced by the tensile energy, ψ is defined as:

ψ(ǫ, d) = g(d)ψ + 0 (ǫ) + ψ - 0 (ǫ) (19)
Where g(d) represents the degradation of tensile energy due to crack formation. Therefore g is a non increasing function and we have g(0) = 1 and g(1) = 0. Moreover, when d = 1 it is asked that the driving force -∂ψ ∂d = -g ′ (d)ψ + 0 (ǫ) be equal to 0. To do so the additional condition g ′ (1) = 0 is imposed. Under these conditions Miehe et al [START_REF] Hofacker | A phase field model for rate-independant crack propagation robust algorithmic implementation based on operator splits[END_REF] proposed the function g(d) = (1 -d) [START_REF] Farzaneh | An experimental study on the tensile properties of atmospheric ice[END_REF] . Using these new relations the complete set of equations becomes :

-div(σ(u, d)) = f vol in Ω σ.n = f surf on Γ 1 u = u d on Γ 2 gc l d -g c l△d = 2(1 -d)ψ + 0 (ǫ) in Ω ∇d.n = 0 on ∂Ω (20)
Equations 20 are nonlinear and describe a stationary damaged equilibrium state compatible with the external constraints. In order to solve these equations Miehe et al [START_REF] Hofacker | A phase field model for rate-independant crack propagation robust algorithmic implementation based on operator splits[END_REF] have proposed an iterative procedure based on the introduction of a history function H. ψ + 0 (ǫ) is then replaced by H in equation 20. The iterative process is solved by a finite element method a and defined by :

• Compute the history field:

H i = max(H i-1 , ψ + 0 (ǫ i-1 )) (21)
• Compute the damage field:

gc l d i -g c l△d i = 2(1 -d i )H i in Ω ∇d i .n = 0 on ∂Ω (22)
• Compute the displacement field:

-div(σ i (ǫ i , d i ) = f vol in Ω σ i .n = f surf on Γ 1 u i = u d on Γ 2 (23)
a As a first approach we have simplified the constitutive law by considering that the damage coefficient affects both the tensile and compressive parts.

If we accept to consider the iterative process as pseudo-unsteady, it can be interpreted as follows: equations 20 only translate conservation of energy. At a given iteration i the tensile energy could locally be inferior to its value at iteration i -1. Thus if, at a given iteration, d = 1, nothing is preventing it from decreasing at any following iteration. But clearly in our case, the problem of crack propagation is irreversible: d should not be allowed to decrease. Thus the history function H serves to take into account the irreversibility of the crack propagation process. It records, locally, the maximum of the tensile energy over all iterations. Intuitively, if at a given iteration there was enough tensile energy to increase the damage variable d, then this information will be contained in H for the following iteration. Nevertheless this is only true in a pseudo-iterative context, further investigation remains to be done in order to determine the exactitude of this interpretation.

V.B.2. Crack propagation test cases

In order to assess the quality of the results obtained via the model presented in section V.B, we performed test cases for which experimental results are available. [START_REF] Bolander | Fracture analyses using spring networks with random geometry[END_REF][START_REF] Aliha | An improved semi-circular bend specimen for investigating mixed mode brittle fracture[END_REF] First we consider a slab of brittle material, which has a vertical notch, on which we apply purely shear loads (mode II). This is one of the test cases presented by Miehe et al . [START_REF] Hofacker | A phase field model for rate-independant crack propagation robust algorithmic implementation based on operator splits[END_REF] We also choose to use the same parameters for this test case. They are summarized in table 2. Figure 9(a) illustrates such a set up. On figure 9(b) we observe that the crack takes a curved path which is the awaited result. Indeed, this result is also observed experimentally. [START_REF] Aliha | An improved semi-circular bend specimen for investigating mixed mode brittle fracture[END_REF] Next, let's consider a slab of brittle material, which has two horizontal notches, on which we apply loads both in tension and shear (mixed mode), as depicted figure 10(a). This time the mechanical properties are those of ice. The imposed displacement has value u d = 5.0 10 -6 m with θ = π/3 and the length scale is set to l = 7.5 10 -5 m The result of the computation is shown figure 10(b). The crack paths are once again those qualitatively observed in experiments using the same setup. [START_REF] Bolander | Fracture analyses using spring networks with random geometry[END_REF]

VI. Simulation

VI.A. Thermal modules

The two thermal modules are explicitly coupled in order to compute the whole temperature field. On figure 11 we can clearly see the formation of the melted region (in white) corresponding to the activation of a heater mat. The snapshot is taken from a fictitious test case and serves only to prove feasibility.

VI.B. Ice shedding feasibility study

Here the goal is to study the feasibility of shedding prediction using the crack propagation model. To do so we consider a simplified situation. We obtain an ice shape and aerodynamic field with the accretion solver of the 2D suite. We then suppress a part of the ice at the leading edge so as to imitate a parting strip. Next we launch a thermal computation to predict the melted regions at the ice/surface interface. We then apply the aerodynamic loads and predict crack propagation. Figure 12 shows a crack that cuts through the ice shape, therefore demonstrating the feasibility of ice shedding prediction using this method.

VII. Conclusion and perspectives

Firstly, we have presented the thermal models and the results they give (coupled or alone). A first step of validation shows that they correctly predict the thermal fields when used individually. A next step will be to perform different academic test cases to span all possible prototype problems. But more importantly we need to compare the coupled computation with experimental results in order to completely assess the validity of the approach.

Secondly, we have demonstrated the feasibility of ice shedding prediction using a new methodology of crack propagation. However this approach remains to be experimentally validated on simple cases as much as on complex de-icing cases.

The next step will now be to integrate these modules into IGLOO2D, the 2nd generation ONERA icing toolbox. Once that will be done, it will be possible to perform a realistic simulation of a complete de-icing cycle in flight conditions. The results of such a computation will be compared with experimental data.

Figure 1 .

 1 Figure 1. The code's architecture

Figure 2 .

 2 Figure 2. Illustration of a leading edge electro-thermal ice protection system

Figure 3 .

 3 Figure 3. Example of temperature field and mesh in the structure

Figure 4 .

 4 Figure 4. Illustration of the hollow cylinder academic case

  (a) At t = 0.2s (b) At t = 2.0s

Figure 5 .

 5 Figure 5. Comparison of non-dimensional analytic and computed temperature profiles

  (a) Phase indicator (b) Mesh

Figure 6 .

 6 Figure 6. Phase indicator and mesh of the ice block

Figure 7 .

 7 Figure 7. Phase change academic case

Figure 8 .

 8 Figure 8. Comparison of melting front positions

Figure 9 .

 9 Figure 9. Computational setup and results

Figure 10 .

 10 Figure 10. Computational setup and results

Figure 11 .

 11 Figure 11. Thermal computation

Figure 12 .

 12 Figure 12. Global view of the fractured ice block

Table 1 .

 1 1) 6 : 17 Elastic constants of interest for isotropic polycrystalline Ice Ih at T = -16 • C

	Property	Units	Value
	Young's modulus, E N.m -2 9.33 10 9
	Poisson's ratio, ν	∅	0.325

Table 2 .

 2 Computation parameters

Acknowledgements

The work presented in this communication was realised in the context of an Airbus/ANRT CIFRE grant. We gratefully acknowledge their support. We also wish to thank F. Dezitter, J. Cliquet and M. Logeais for the helpful discussions we had.