Evgeny Glushkov

Natalia Glushkova

Artem Eremin

Victor Giurgiutiu

ANALYTICALLY-BASED SIMULATION FOR CORROSION DETECTION BY GUIDED WAVES

Keywords: guided waves, diffraction, corrosion detection, layer and beam models

To create a theoretical basis for guided wave detection and identification of corrosion damages, a set of analytically based computer models of various complexity has been developed. The present paper is focused on the simplest and fastest beam model for stepped and notched waveguides, which has exhibited a wide frequency range of reasonable coincidence with the results obtained within more complex integral equation based model for a 2D notched elastic strip.

INTRODUCTION

Guided wave (GW) detection and identification of corrosion damages is still a challenging task for SHM [START_REF] Swenson | Comparisons of analytical and experimental measurements of lamb wave interaction with corrosion damage in aluminum plates[END_REF]. Corrosion areas are less contrast than "conventional" cracks, delaminations and disbondings which give strongly localized response. Therefore, the extraction of corrosion indications from received signals (both scattered by and getting through the corrosion area) requires more comprehensive processing based on fast computer simulation of GW diffraction by surface irregularities. Commercial FEM packages are rather time-consuming, especially with 3D scattering, so it is worthy to select as simple models as possible, which, nevertheless, capture the characteristic diffraction features of areas affected by corrosion.

Motivated by that idea, we have developed a hierarchy of 2D and 3D semi-analytical models of varying complexity, comparing them with each other, with benchmark results of other authors and with experimental measurements, estimating in this way the range of applicability of each model. In descending order of complexity there are 1) laminate element method (LEM ) for 3D scattering by depressions and cavities [START_REF] Ye | The layered element method in the dynamic theory of elasticity[END_REF] ; 2) LEM based models for GW diffraction by arbitrarily shaped notches in 2D elastic strip waveguides; 3) expansion in series in propagating and evanescent normal modes for 2D stepped waveguides [START_REF] Glushkov | Diffraction of normal modes in composite and stepped elastic waveguides[END_REF]; 4) eccentrically butted beams of different thickness as 1D stepped and notched waveguides.

The closest analogue of the first LEM based model is the approach developed by Moreau et al [START_REF] Moreau | Scattering of guided waves by flat-bottomed cavities with irregular shapes[END_REF]. Numerical comparisons with the GW scattering diagrams obtained by that method exhibit a full coincidence. The second and third 2D models have been tested against the numerical and experimental results by Lowe et al [START_REF] Lowe | The low frequency reflection characteristics of the fundamental antisymmetric lamb wave a 0 from a rectangular notch in a plate[END_REF] as well as by checking the boundary conditions at the joint line and the energy balance among the incident, reflected, transmitted and converted modes. In turn, these models were used to estimate the range of applicability of the simplest and fastest beam models. The latter has demonstrated an unexpectedly wide frequency range in which the reflection and transmission coefficients reasonably coincide with those for the fundamental S 0 and A 0 modes in 2D stepped and notched waveguides.

The present paper is focused on this beam model: first, a mathematical model for GW propagation in pristine and damaged Euler-Bernoulli beams is given and quantitative energy characteristics of corresponding waves are introduced. After that a corresponding 2D LEM based approach for a notched plate is briefly described to verify the applicability limits of the simplified model. Finally the results obtained by the both approaches are compared and discussed.

BEAM MODEL

An elastic beam of width b and thickness h occupies the volume |x| < ∞, |y| ≤ b/2 and |z| ≤ h/2 in the Cartesian coordinate system. Its deformation is independent of y, so it is specified by the independent of y displacement vector u = (u x , 0, u z ) lying in the central plane (x, z) (2D plane-strain deformation, Fig. 1a). Within the beam assumptions

x z h/2 -h/2 0 u,w a) b) h 1 B 1 h 1 h 2 B 3 B 2 x x x
u x (x, z) = u(x) -zw ′ (x) (1) u z (x, z) = w(x), (2)
where u(x) and w(x) are 1D functions set on the beam axis -∞ < x < ∞, y = 0, z = 0. With a steadystate time-harmonic oscillation ue -iωt , the vector u is a complex vector of displacement amplitude, and the functions u(x) and w(x) obey the beam equations

u ′′ + ζ 2 1 u = 0 (3) w IV -ζ 4 2 w = 0. ( 4
)
Here

ζ 2 1 = ρ Y ω 2 and ζ 4 2 = ρA Y I ω 2 = 12ρ Y h 2 ω 2 ,
Y is the modulus of elasticity (Young modulus), ρ is density, A = bh is the area of beam's cross-section, I =

∫∫ A z 2 dydz = bh 3 /12 is the moment of inertia, ω = 2π f is angular frequency, f is frequency.
The beam supports longitudinal and flexural guided waves S 0 and A 0 , which are specified by the eigensolutions of Eqs. (3) and (4), respectively:

S 0 : u 0 (x) = S 0 e iζ 1 x
(5)

A 0 : w 0 (x) = A 0 e iζ 2 x (6)
The amplitudes S 0 and A 0 are arbitrary complex constants, ζ 1 and ζ 2 are the wavenumbers of these waves. Besides, the evanescent terms Be ±ζ 2 x also satisfy Eq. ( 4).

Force resultants, e.g., longitudinal force N, moment M and shear force V are expressed in terms of the displacement components:

N(x) = ∫∫ A σ x dydz = YAu ′ (x), M(x) = ∫∫ A zσ x dydz = -Y Iw ′′ (x) (7) V (x) = M ′ (x) = -Y Iw ′′′ (x). (8)
Here σ x (x) = Y ε x = Y (u ′ -zw ′′ ) is cross-sectional traction. The corresponding shear stress τ xz is estimated via the assumption of even distribution of the force V over the cross-section A:

τ xz (x) = V /A = -Y I A w ′′′ (x). (9)

WAVE ENERGY

In a time-harmonic wave field, the energy flux is estimated in terms of quantities averaged over the period of oscillation T = 2π/ω = 1/ f . The averaged density and the direction of energy flux passing through a spatial point x per unit time is specified by the energy density vector e(x) = (e x , e y , e z ). The total amount of energy E, carried by harmonic waves through a certain surface S per unit time (in fact, the power of the energy flux), is determined via the integration of the energy density over S:

E = ∫∫ S e n (x)dS. (10)
Here e n = (e, n) = ω 2 Im(τ n , u) is the normal to S component of vector e; n(x) is the unit surface normal at the current point x ∈ S; τ n is the stress vector at a surface element specified by n; the scalar product of complex vectors assumes the complex conjugation of the second factor hereinafter denoted with asterisk:

(a, b) = ∑ i a i b * i .
For the calculation of wave energy carried by S 0 and A 0 guided waves along the beam, one has to take its cross-section A as the surface S in Eq. (10). At that, the normal n = (1, 0, 0), τ x = (σ x , 0, τ xz ), and E can be obtained as follows

E(x) = b ∫ |z|<h/2 e x (x, z)dz = ω 2 bY Im ∫ |z|<h/2 [u ′ u * + z 2 w ′′ (w ′ ) * -I A w ′′′ w * ]dz = = ω 2 Y Im [Au ′ u * + Iw ′′ (w ′ ) * -Iw ′′′ w * ]. (11)
In view of the energy conservation law, E(x) must be constant in the segments of ideally elastic waveguides free from wave sources and energy drains irrespective of their thickness variation. For travelling waves ( 5) -( 6)

u ′ 0 u * 0 = iζ 1 |S 0 | 2 and w ′′ 0 (w ′ 0 ) * = -w ′′′ 0 w * 0 = iζ 3 2 |A 0 | 2 , thus, E = E S + E A = ω 2 YAζ 1 |S 0 | 2 + ωY Iζ 3 2 |A 0 | 2 . ( 12
)
The parts E S and E A are energy of S 0 and A 0 modes, respectively. They independently contribute into the the total amount of wave energy E transferred through the cross-section x = const per time unit. Remark 1. In the case of waves propagating in opposite directions, they contribute into E with opposite signs:

E S = ω 2 YAζ 1 (|S + | 2 -|S -| 2 ) and E A = ωY Iζ 2 2 (|A + | 2 -|A -| 2 )
for u = S + e iζ 1 x + S -e -iζ 1 x and w = A + e iζ 2 x + A -e -iζ 2 x . Note that mixed terms such as S + (S -) * e 2iζ 1 x appear as complex conjugate pairs. Their sums, being complex conjugate values, do not affect the amount of energy E calculated via Eq. ( 11). Remark 2. The evanescent terms also do not affect the amount of energy carried by GWs. For example, with w = Ae iζ 2 x + Be -ζ 2 , the mixed terms in w ′′ (w ′ ) *w ′′′ w * in Eq. ( 11) are real values of

form 2ζ 3 2 [Im (iAB * e iζ 2 x ) -Re (A * Be -iζ 2 x )]e -ζ 2 x .

BOUNDARY CONDITIONS IN THE DOCKING AREA

Let us consider a beam with a rectangular notch of length ∆x and depth ∆z. It may be treated as a joint of three beams B 1 , B 2 and B 3 of, generally speaking, three different thicknesses h 1 , h 2 and h 3 .

In the case under consideration

h 3 = h 1 , ∆x = x 2 -x 1 , ∆z = h 1 -h 2 ,
and the value e = ∆z/2 is called "the eccentricity" of the non-coaxial beam B 2 (Fig. 1b). The central axis of the latter goes through the point O 2 (x 1 , -e) in the global coordinate system (x, z). Therefore, the displacements in B 2 should be written in the local coordinates (x, y, z 2 ), where z 2 = z + e:

u x,2 = u 2 -z 2 w ′ 2 u z,2 = w 2 (13)
In these coordinates the moment and shear force are

M2 (x) = ∫∫ A 2 z 2 σ x,2 dydz 2 = -Y I 2 w ′′ 2 and V 2 (x) = M′ 2 (x) = -Y I 2 w ′′′ 2 . ( 14
)
It is necessary to consider that M2 is calculated with respect to the local center O 2 , while the calculation of the moment M 2 with respect to the same as in B 1 center O 1 leads to the additional, proportional to the eccentricity e, term:

M 2 (x) = ∫∫ A 2 (z 2 -e)σ x,2 dydz 2 = M2 -eN 2 , N 2 = YA 2 u ′ 2 (x). (15)
Note that the force V 2 is independent of coordinate system, it remains of form (14).

To formulate proper boundary conditions at the joining points x 1 and x 2 , one should equate the displacements, forces and moments at these points. The displacement equality at x = x 1 entails three equalities connecting u 1 and w 1 with u 2 and w 2 :

u 1 = u 2 -ew ′ 2 w 1 = w 2 w ′ 1 = w ′ 2 , x = x 1 , ( 16
)
while the force and moment equalities

N 1 = N 2 , V 1 = V 2 ,
and M 1 = M 2 at x = x 1 yield more three links:

h 1 u ′ 1 = h 2 u ′ 2 h 3 1 w ′′ 1 = h 3 2 w ′′ 2 + 12eh 2 u ′ 2 h 3 1 w ′′′ 1 = h 3 2 w ′′′ 2 , at x = x 1 ( 17
)
The same way, the conditions at x = x 2 are

u 2 -ew ′ 2 = u 3 w 2 = w 3 w ′ 2 = w ′ 3 h 2 u ′ 2 = h 1 u ′ 3 h 3 2 w ′′ 2 + 12eh 2 u ′ 2 = h 3 1 w ′′ 3 h 3 2 w ′′′ 2 = h 3 1 w ′′′ 3 ( 18
)
EWSHM 2014 -Nantes, France Remark 3. Boundary conditions ( 16) -( 17) and (18) assure the energy conservation in the course of its transfer through the butt joints. This can make explicit, e.g., taking into account that in line with the conditions at x = x 1

A 1 u ′ 1 u * 1 = A 2 u ′ 2 u * 2 -A 2 eu ′ 2 (w ′ 2 ) * and I 1 w ′′ 1 (w ′ 1 ) * = I 2 w ′′ 2 u * 2 + A 2 eu ′ 2 (w ′ 2 ) * .
The second terms with the eccentricity e are reduced in the sum of Eq. ( 11), hence, the expression for the amount of energy E 1 at the left of x 1 becomes the same as E 2 at the right side of the joint. Similarly for x = x 2 .

REFLECTION AND TRANSMISSION COEFFICIENTS

The diffraction of an incident A 0 or S 0 wave by the notch gives rise to the reflected and transmitted fields u - 1 and u + 3 . Thus the general solution of Eqs. ( 3), ( 4) in the whole domain B = B 1 ∪ B 2 ∪ B 3 may be represented in the following form

u 1 = u 0 + u - 1 = u 0 + c 1 e -iζ 1 (x-x 1 ) , x ∈ B 1 u 2 = c 2 e iζ 1 (x-x 1 ) + c 3 e -iζ 1 (x-x 2 ) , x ∈ B 2 u 3 = u + 3 = c 4 e iζ 1 (x-x 2 ) , x ∈ B 3 ( 19
)
w 1 = w 0 + w - 1 = w 0 + c 5 e -iζ 2,1 (x-x 1 ) + c 6 e ζ 2,1 (x-x 1 ) , x ∈ B 1 w 2 = c 7 e iζ 2,2 (x-x 1 ) + c 8 e -ζ 2,2 (x-x 1 ) + c 9 e -iζ 2,2 (x-x 2 ) + c 10 e ζ 2,2 (x-x 2 ) , x ∈ B 2 w 3 = w + 3 = c 11 e iζ 2,3 (x-x 2 ) + c 12 e -ζ 2,3 (x-x 2 ) , x ∈ B 3 (20)
The 12 unknown constants c j , j = 1, 2, ..., 12, are obtained from the 12 × 12 system of linear equations resulting from the 6 + 6 joining conditions ( 16) -( 18).

The energy E -and E + of the reflected and transmitted waves u -and u + may be obtained using the same formulas as Eqs. (12) but with the coefficients c 1 and c 5 substituted for A 0 and S 0 in the case of reflected wave energy E -= E - S + E - A :

E - S = ω 2 YA 1 ζ 1 |c 1 | 2 , E - A = ωY I 1 ζ 3 2,1 |c 5 | 2 , ( 21
)
and with the amplitude constants c 4 and c 11 for the transmitted wave energy

E + = E + S + E + A : E + S = ω 2 YA 2 ζ 1 |c 4 | 2 , E + A = ωY I 2 ζ 3 2,2 |c 11 | 2 . ( 22
)
The values E ± S and E ± A are energy transferred by transmitted and reflected S 0 and A 0 GWs. With an ideally elastic waveguide structure, the energy balance is

κ -+ κ + = 1, ( 23
)
where the transmission and reflection coefficients κ ± = E ± /E 0 consist of specific coefficients for the scattered S 0 and A 0 waves:

κ ± = κ ± S + κ ± A = E ± S /E 0 + E ± A /E 0 ( 24
)
If only S 0 mode is taken as the incident field (A 0 = 0), the values κ ± S may be treated as the coefficients of its transmission and reflection (κ ± S ≡ κ ± SS ), while κ ± A are the coefficients of its forward and backward conversion into A 0 modes (κ ± A ≡ κ ± SA ). Similarly, with an A 0 incidence (S 0 = 0), κ ± A ≡ κ ± AA and κ ± S ≡ κ ± AS .

2D NOTCHED WAVEGUIDE

As a more complex model for GW diffraction by corrosion areas, a 2D notched elastic strip (Fig. 1b) governed by the full system of elastodynamic equations

(λ + 2µ)div u + µ∆u + ρω 2 u = 0 (25)
has been also considered; λ and µ are Lamé constants of elasticity. In this statement, the coupling boundary conditions are different from those for the beams. They are the conditions of displacement and stress field continuity at the joint interface lines

[u] m = 0, [τ x ] m = 0, x = x m , -h 1 /2 ≤ z ≤ h 1 /2 -∆z, m = 1, 2 (26)
and the stress-free conditions

τ x = 0, x = x m , h 1 /2 -∆z ≤ z ≤ h 1 /2, m = 1, 2 (27)
at the rest of edges of the thicker strips B 1 and B 3 , not contacting with the thinner intermediate domain B 2 . Square brackets denote here the jump of related vector functions at the cross-sections x = x m . The horizontal sides of the notched domain are also stress-free. The diffraction of an incident guided wave u 0 (x) = a 0 (z)e iζ k x gives rise to the scattered field u sc (x), so that the total wave field in the notched waveguide is u = u 0 + u sc . Here ζ k , k = 1 or 2, is the wavenumber of the fundamental S 0 or A 0 Lamb wave; these values tend to the beam's wavenumbers in Eqs. ( 3) -( 4) as ωh √ ρ/Y → 0; a 0 (z) is an eigenform of the corresponding Lamb wave. The scattered field is derived using the LEM technique [START_REF] Ye | The layered element method in the dynamic theory of elasticity[END_REF] in the form

u sc (x) = ∫ S l(x, ξ )c(ξ )dξ , ( 28
)
where S is the boundary of the notch, c(x) is an unknown potential's density and l(x, ξ ) = [l 1 . . . l 2 . . .

l 3 ]
is the matrix of fundamental solutions for the intact elastic layered structure under consideration. Its columns l j , j = 1, 2, 3, are displacements generated by the point sources δ (x -ξ )i j , where i j are basic coordinate vectors. They satisfy the governing equations and the homogeneous boundary conditions at all plane-parallel surfaces (at the surfaces z = ±h 1 /2 in the case). A substitution of relation (28) into the boundary conditions on S yields a boundary integral equation with respect to the vector function c. Its solution is obtained using boundary element method technique.

NUMERICAL EXAMPLES AND DISCUSSION

To estimate the range of practical applicability of the results obtained within the beam model, a systematic comparison of the transmission, reflection and conversion coefficients κ ± obtained versus frequency within beam and LEM models have been carried out. Preliminary the LEM model has been validated against the known FEM and experimental results [START_REF] Lowe | The low frequency reflection characteristics of the fundamental antisymmetric lamb wave a 0 from a rectangular notch in a plate[END_REF]. Figure 2 gives examples of such comparisons for the amplitude reflection coefficient µ -= |w -/w 0 | of A 0 mode propagating in a steel plate of thickness h 1 = 3 mm. The left subplots show the frequency dependence of µ -for two notches of width ∆x = h 1 and depths h 2 = 0.83h 1 and 0.5h 1 . The right subplot illustrates the influence of the notch width ∆x variation. The reflection coefficient µ -is shown here versus the ratio ∆x/λ , where λ = 5.5 mm is the A 0 wavelength at f = 450 kHz. Figures 3 and4 give examples of beam-to-LEM comparisons for two aluminium samples of thickness h 1 = 1 mm with the notches of length ∆x = 5 mm and depths ∆z = 0.5 and 0.75 mm (h 2 = 0.5 and 0.25 mm); the material properties Y = 71 GPa, ρ = 2700 kg/m 3 and Poisson's ratio ν = 1/3. The figures are for S 0 and A 0 incidence. To avoid cluttering the figures, only transmission (blue) and reflection (red) coefficients κ + and κ -are shown without the coefficients of mode conversion, exhibiting similar rate of coincidence. At relatively low frequencies h f /v << 1 (v is a characteristic wave velocity), the displacement field u in an elastic strip exhibits linear dependence on the cross coordinate z, as it is formulated in beam assumptions (1) -( 2). The fundamental A 0 and S 0 Lamb waves are also well approximated by beam's guided waves (1), ( 2), ( 5) and (6).

On the other hand, linear behaviour of displacement and stress fields relative to z coordinate is violated near the lines of butt junction even in the limit f → 0. Therefore, it is not clear how well the conditions of beam coupling may substitute for conditions (26) -( 27). One more factor differing beam and strip models is that though the S 0 and A 0 modes of the latter become of form ( 1) -( 2), ( 5) -( 6) at low frequencies, their wavenumbers are slightly different from ζ 1 and ζ 2 diverging from them as frequency increases. Nevertheless, the comparison of transmission and reflection coefficients κ ± S and κ ± A obtained for notched beams and strips has shown that at low frequencies the beam plots follow reasonably close to the strip counterparts (Figs. 3,4). The coincidence of results for the S 0 incidence is visibly better then in the A 0 case. In the frequency range up to 300 kHz, it is very good. In the A 0 case the coincidence is not so good. Although, even with the most complicate curve behaviour at h 2 = 0.25 (Fig. 4, right), the beam curves catch the peaks and minima in the range f ≤ 100 kHz, following them with a certain shift in a wider range as well.

One of the reasons for better beam model performance with S 0 incidence is the much better wavelength-to-thickness (λ /h) ratio than with A 0 waves. For example, at the frequencies f = 100 kHz and 500 kHz, these ratios for S 0 are λ S /h = 51.3 and 10.3, while for A 0 they are λ A /h = 9.6 and 2.9, respectively.

CONCLUSION

In spite of apparent simplicity of the Bernoulli-Euler beam equations, it is quite acceptable to use them for the simulation of GW diffraction by step and notch obstacles at low frequencies, at least in the first quarter of the two-mode Lamb wave range. The crucial point here is the formulation of coupling boundary conditions at eccentric butt joints. Only a proper accounting for the eccentricity provides the wave energy conservation across the junction as well as correct energy partition among the reflected and transmitted travelling waves.

Figure 1 :

 1 Figure 1 : Geometry of problem: pristine beam (top) and notched waveguide (bottom).

Figure 2 :

 2 Figure 2 : Examples of LEM validation against FEM and experimental results [5].

Figure 3 :

 3 Figure 3 : Transmission (blue) and reflection (red) coefficients κ + S and κ - S obtained within LEM and beam models (solid and dashed lines, respectively) for the S 0 mode in notched waveguides of web thicknesses h 2 = 0.5 and 0.25 mm; h 1 = 1 mm, ∆x = 5 mm.

Figure 4 :

 4 Figure4: Same as in Fig.3for A 0 incidence.

EWSHM 2014 -Nantes, France

The authors are grateful to Prof. W. Seemann, KIT, Karlsruhe, for the useful discussion of the beam model. The work is partly supported by the Russian Foundation for Basic Research (RFBR) (projects No. 12-01-00320 and 14-08-00370).


Jean-Pierre Monchau

Laurent Ibos

Vincent Feuillet
email: vincent.feuillet@u-pec.fr

DIAGNOSIS OF INSULATED BUILDING WALLS USING PASSIVE INFRARED THERMOGRAPHY AND NUMERICAL SIMULATIONS

Keywords: building monitoring, passive infrared thermography, thermal modeling, thermal resistance

This work presents the thermal monitoring of a multi-layered wall of a restored building (PANISSE platform). Surface temperatures measured by infrared thermography are corrected by taking account of the influence parameters. Then they are compared with thermocouple measurements and numerical simulations.

INTRODUCTION

Infrared thermography is frequently used for building monitoring because of several advantages. It is a non-intrusive measurement method that leads to high spatial and thermal resolutions. Despite that, infrared thermography is not commonly used for quantitative diagnosis because of the difficulty to understand in details all the influence parameters of the thermal scene to be viewed [START_REF] Bareira | Evaluation of building materials using infrared thermography[END_REF].

Currently there is only one international norm that defines an in-situ measurement method of the thermal resistance of insulated building walls based on the use of contact sensors (thermocouples and fluxmeters) [2]. However a new norm based on the use of infrared thermography is in preparation [START_REF] Kato | Method of in-situ measurement of thermal insulation performance of building elements using infrared camera[END_REF]. In literature, most studies involving passive infrared thermography in building monitoring leads to qualitative results. For example insulation defects are detected on building façades [START_REF] Balaras | Infrared thermography for building diagnostics[END_REF] or the influence of the measurement hour in the day is analyzed [START_REF] Ocaña | Thermographic survey of two rural buildings in Spain[END_REF]. Only a few quantitative diagnoses by means of passive infrared thermography exist. A complete methodology based on the infrared thermovision technique has been developed to reach quantitative data of thermal transmittances of building envelopes [START_REF] Albatici | Infrared thermovision technique for the assessment of thermal transmittance value of opaque building elements on site[END_REF]. However authors insist on the high influence of the measurement conditions (weather, managing of building by users, maintenance works, ageing of the building materials). Another interesting approach consists in investigating both visual and IR images in the same study [START_REF] Ribarić | A knowledge-based system for the non-destructive diagnostics of façade isolation using the information fusion of visual and IR images[END_REF] to yield a diagnosis of various façades with a significant reduction of the time analysis per evaluation.

In this study passive infrared thermography measurements on a multi-layered wall (PANISSE platform) are analysed. This building has been restored with an outside thermal insulation in 2011. A set of parameters is measured continuously: outside and inside air temperatures and humidities, surface wall temperatures. A fluxmetric method is used to measure the thermal resistance of the wall. Finally infrared temperature measurements are compared after correction with numerical simulation results.

INSTRUMENTATION AND MEASUREMENTS

Description of the PANISSE platform

The PANISSE platform is a residential building with two floors located in the town of Villemomble, at about ten kilometres in the east of Paris. It was built in the sixties and walls are made of cellular concrete. During a renovation in 2011, a 6cm-thick external insulating layer was fixed onto eastern, northern and western façades. Southern wall is a shared with a symmetrical house. Insulating material used is graphite expanded polystyrene. An air gap of about 1cm thick was kept between concrete wall and polystyrene layer. A schematic description of the wall structure is provided in Figure 1. The house is a lived-in home and its terrace roof is not thermally insulated. In the northern façade, a garage is used as a thermal guard (see Figure 2). A weather station was installed on the building roof. This device provides several parameters: wind direction and speed, air temperature, relative humidity, atmospheric pressure, rainfall amounts. During renovation works, K-type thermocouples were fixed onto several positions of the three free façades: 3 onto the eastern façade, 1 onto the northern façade and 6 onto the western façade (see Figure 2 for details). For each position, one thermocouple was fixed between the polystyrene and façade coating layers and another one was fixed onto the external face of the cellular concrete wall, as presented in Figure 1. Solar cells were fixed onto each of the free façades of the house and on the roof in order to measure the incident solar flux density on each façade. These solar cells were preliminarily calibrated by using a reference Kipp&Zonnen™ pyranometer. Inside wall surface temperature is not measured. However, Warito™ data-loggers are used to record relative humidity and air temperature in several rooms of the house with a sampling period of 10 minutes. These data-loggers are provided with a calibration certificate: uncertainties on temperature and relative humidity are of 0.12°C and 3% respectively. Configuration and recorded data downloading is performed periodically through a USB interface. Some preliminary measurements were performed to characterize optical properties of the coating used on each façade. The solar flux absorption coefficient was in-situ determined using a home-made albedometer. This device is composed of two opposite solar cells calibrated with a pyranometer. The ratio of both reflected and incident fluxes gives the albedo of the façade surface. The obtained value of absorption coefficient is α = 0.37. In the same way, the emissivity of the façades coating was in-situ determined using a portable emissometer recently developed in the CERTES laboratory [START_REF] Monchau | Infrared emissivity measurements of building and civil engineering materials[END_REF]. The emissivity value in the 8-14 µm band is ε = 0.93.

Infrared thermography measurement campaign

A measurement campaign by infrared thermography was performed in 2013 during 12 consecutive days between January 25 and February 6. The camera used is a FLIR™ SC 7300 LWIR camera (7.7-9.2 µm) with a Stirling cooled detector of 320 x 256 pixels. A standard 25 mm lens was used. Thus, the IFOV value is 1.3 mrad. The camera was connected through an Ethernet interface to a portable computer. Images were recorded each 15 seconds during 12 days using the Altair™ software. The infrared camera and computer were placed inside a watertight enclosure (see Figure 3). This enclosure was installed in the garden of the house in order to observe the western façade. An example of thermal image recorded during a sunny day is presented in Figure 3. Air temperature and relative humidity inside the enclosure were continuously monitored using a Warito™ datalogger. A diffusive mirror of rough aluminium of 0.063 emissivity [START_REF] Datcu | Contribution to improvement of building wall surface temperature measurements by infrared thermography[END_REF] and a black surface coated with a Nextel Velvet 811-21 coating of 0.97 emissivity were fixed on the upper left part of the investigated façade. Temperatures of both surfaces were recorded all along the measurement campaign using a numeric thermometer and two K-type thermocouples. The diffusive mirror allows computing the mean radiant temperature of the environment considered as a blackbody [START_REF] Datcu | Contribution to improvement of building wall surface temperature measurements by infrared thermography[END_REF]. The black reference surface allows estimating the atmospheric attenuation.

Thermal image sequences are then corrected in order to take into account all parameters influencing computed values of the external wall surface temperature. These parameters are wall surface emissivity, mean radiant temperature, atmospheric transmission factor and atmospheric temperature [START_REF] Monchau | Mesure d'émissivité pour la thermographie infrarouge appliquée au diagnostic quantitatif des structures[END_REF]. An example of temperature variations recorded during the whole campaign by thermocouple and infrared camera is presented in Figure 4. These measurements concern the measurement zone in the left upper part of observed façade. The differences between infrared camera and thermocouple temperature measurements for three regions of the façade are plotted in Figure 5. Mean difference values obtained in that cases are -0.07°C, 0.15°C and 0.14°C respectively and corresponding standard deviations are 0.44°C, 0.43°C and 0.39°C. Thus, a satisfactory agreement between thermocouple and infrared camera measurements is obtained excepted for some cloudy and sunny alternating short periods.

MEASUREMENT OF THE THERMAL RESISTANCE OF THE WALL

Principle

Two methods are used to measure the thermal resistance of the wall: the first one is based on the ISO-9869 norm [2], the second one uses the knowledge of the wall composition. The insulating layer is made of graphite expanded polystyrene that is 62 mm in thickness and 2 m 2 .K.W -1 in thermal resistance.

Figure 6 shows measurement points used. T 1 is the interior wall surface temperature, T 2 is the surface temperature between the cellular concrete and the air gap and T 3 is the exterior wall surface temperature. These measurements are done with K-type thermocouples. Interior and exterior air temperatures are measured with Warito TM sensors. HFM is a heat flux sensor from Captec TM (sensitivity 18 µV.W -1 .m 2 ).

Measurement with ISO-9869 norm

The measurement consists in computing time averages of the wall surface temperatures and of the heat flux. Time averages are ideally calculated on several days. One must check that accumulated energy in the wall is negligible compared to the energy transferred through the wall. A simple experimental approach is used: thermal resistance R th is computed every 24h and the obtained value is confirmed if differences between two consecutive days are less than 5%. ISO-9869 describes a dynamic approach allowing taking account of the wall thermal inertia. The presented measurement lasted 48h. This could have been insufficient if meteorological conditions have not been suitable: cloudy sky, no sunshine, quite stable temperatures between day and night and equivalent conditions during the days before the measurement.

The total thermal resistance R th is given as follows:

R th = T 1 ! T 3 " S (1)
with Φ s is the heat flux per unit surface area.

Method using the wall composition

The thermal resistance of the insulating layer is used to estimate the heat flux:

! S = T 2 " T 3 R ins + R gap (2)
with R ins is the thermal resistance of the insulating layer and R gap the thermal resistance of the air gap. Then Φ s is used to compute R th with eq. ( 1).

Results

The average values of the two-day measurements are given in Table 1: The computed heat flux with the second method gives Φ s =3.46 W.m -2 . Thermal resistances obtained with the two methods are respectively: R th1 =3.83 m 2 .K.W -1 et R th2 =3.44 m 2 .K.W -1 . The difference between the two methods is about 10%, which is reasonable for a fluxmetric method.

COMPARISON BETWEEN MEASUREMENTS AND NUMERICAL SIMULATIONS

Principle

In this part measurements are compared with simulation results obtained with a modeling of the heat transfer in the wall. Analysis concerns the measurements on the first floor of the western façade. The numerical simulation is achieved with commercial software Comsol TM that solves the studied thermal system thanks to the finite element method.

First the heat transfer modeling is described in terms of domains and boundary conditions. The different parameters of the modeling are also defined. The simulated wall is presented in Table 2 along with its geometrical and thermophysical properties. An air gap is put between the cellular concrete and the expanded polystyrene. Figure 7 presents a schematic description of the wall considered in the thermal modeling. The wall is submitted to heat exchanges both on interior and exterior surfaces, respectively with constant global heat transfer coefficients h int and h ext and air temperatures T int (t) and T ext (t). We consider mono dimensional heat transfer by conduction in the zdirection through the layers of the wall (thickness e i , thermal conductivity λ i , density ρ i , heat capacity c i ). Solar heat flux Φ 0 (t) is absorbed by the exterior wall surface by considering the absorption factor α of the surface alpha. The air gap is modeled by taking account of a thermal resistance R gap between the concrete and the insulating layer. Indeed the investigation area is taken in the center of the wall, distant from windows and thermal bridges. Heat transfer parameters from literature [START_REF]Réglementation Thermique[END_REF] or measurement (coefficient α) are reported in Table 3.

Results

The exterior wall surface temperature measured by IR thermography is compared with the temperature calculated with modeling presented in 3.1. The measured meteorological parameters are used in the modeling. The Figure 8(a) shows measured and calculated temperatures T mea and T cal for 7 days of measurement. The Figure 8(b) shows the corresponding differences T mea -T cal .

We note a correct similarity between measurement and modeling temperature profiles for the considered period. However a bias is observed along the investigated period, the measured temperature is globally inferior to the calculated temperature. The temperature difference average is -0.71°C while the standard deviation is 0.63°C. The most important discrepancies appear when rapid variations of the solar flux are observed. The next step of the study will consist in using a parameter estimation method to decrease these differences by adjusting some modeling parameters.

CONCLUSION

This work presents the instrumentation of a lived-in home allowing the thermal monitoring of the walls of the building. The first step of the study consists in correcting the wall surface temperatures measured by infrared thermography by taking account all the influence parameters (wall surface emissivity, mean radiant temperature, atmospheric transmission factor and atmospheric temperature). The proposed correction method is evaluated by comparing with K-type thermocouple measurements. The following step shows a correct similarity between measured exterior wall surface temperatures and modeling results calculated with the finite elements method.

The next step will consist in estimating unknown parameters in the modeling in order to determine the thermal resistance of the wall. This value will be compared with the one obtained with the fluxmetric method. This identification is achieved by solving the corresponding inverse problem formulated in the least-squares sense. It consists in minimizing the functional that links the measured and calculated temperatures by using Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for the least squares estimation of non linear parameters[END_REF].

Figure 1 :

 1 Figure 1: Schematic description of the PANISSE walls structure

  All sensors are connected to a National Instruments™ data acquisition device controlled by a LabView™ application. Data are recorded with a sampling period of 15 seconds since the beginning of 2013. Parameters values are stored in daily time stamped ASCII files.

Figure 2 :

 2 Figure 2: Views of the PANISSE platform; eastern and northern façades (left), northern and western façades (right); red points indicate the positions of thermocouples

Figure 3 :

 3 Figure 3: View of the thermographic instrumentation (left) and example of thermal image recorded during a sunny day (right)

Figure 4 :

 4 Figure 4: Comparison between thermocouple and corrected infrared measurements for the left upper zone of the western façade

Figure 5 :

 5 Figure 5: Temperature differences between thermocouple and corrected infrared measurements for three zones of the western façade

Figure 6 :

 6 Figure 6: Schematic representation of the instrumentation used to measure the thermal resistance wall

Figure 7 :

 7 Figure 7: Schematic representation of the wall considered in the thermal modeling

Figure 8 :

 8 Figure 8: Exterior wall surface temperatures -first floor of the western façade (a) Measured and calculated temperatures T mea and T cal , (b) Differences T mea -T cal

Table 1 :

 1 Average value parameters to measure the thermal resistance of the wall

		Average values of the two-day measurements
	T int (°C)	T 1 (°C)	T 2 (°C)	T 3 (°C)	T ext (°C) Φ S (W.m -2 )
	19.66	18.34	16.00	6.42	6.28	3.11

Table 2 :

 2 Geometrical and thermophysical properties of the wall studied in the modeling

	Layer i	Thickness e i (cm)	Thermal conductivity λ i (W.m -1 .K -1 )	Density ρ i (kg.m -3 )	Heat capacity c i (J.K -1 .kg -1 )
	Cellular concrete (i=1)	20	0.14	400	1000
	Expanded polystyrene (i=2)	6.2	0.031	20	1450
	Coating (i=3)	0.3	0.3	833	1000

Table 3 :

 3 Heat transfer parameters used in the modeling

	h int (W.m -2 .K -1 ) h ext (W.m -2 .K -1 )	α	R gap (m 2 .K.W -1 )
	25	8	0.37	0.15

July 8-11, 2014. La Cité, Nantes, France Copyright © Inria (2014)

EWSHM 2014 -Nantes, France


P Pons
email: ppons@laas.fr

H Aubert

E Debourg

A Rifai

M Olszacki

M Matusiak

I Augustyniak

P Knapkiewicz

J Dziuban

D Lavielle

C Chatry

CHIPLESS PASSIVE SENSOR FOR WIRELESS MONITORING OF HIGH RADIATION DOSES IN NUCLEAR INFRASTRUCTURES

Keywords: Radiation dosimeter, polymer radiolysis, MEMS

The dosimetry is one of the crucial techniques that are needed to assure personal safety and facilities security in the areas of high radioactivity as nuclear power plants or powerful experimental research infrastructure (LHC, XFEL, ITER). Until now, the common electronic dosimeters are based on the silicon diodes or field effect transistors what limits the measured dose value up to few Mrad, as the higher levels lead to the device saturation. The other techniques that allows for higher dose measurement (up to 100 Mrad and more) are based mostly on the variation of physic-chemical parameters of different materials but need the post measurement treatment of the sensor that has to be taken out of the monitored zone and use of the special equipment what makes the continuous measurement impossible. Moreover, in case of nuclear reactors the access to the radioactive zones is limited what leads to the measurement frequency in range of a few times per year. The solution that may overcome that problem is the use of the recently developed innovative concept of the passive sensors that are interrogated by the radar technique. Such a solution uses the small (millimeter sized) micro-resonators designed for the ultra high frequency (about 30GHz) that are connected to the antenna and serve as its load. Any change of the resonant frequency of the resonator may be detected as the shift in the radar echo of the antenna. The sensor is interrogated by the directional radar and the reflected signal is analyzed. The associate transduction principle is that the resonator changes its resonant frequency due to the variation of the specific measured physical value like temperature, pressure or humidity. In our case, we proposed to use the known principle of Hydrogen-pressure dosimeters (HPD) where the polymer material dehydrogenates under irradiation. However, using the MEMS technology we want to miniaturize it and combine with the RF pressure transducer, in order to obtain the miniature (about 4 mm square) passive dosimeter. The main advantages of such a system are the passive (battery less) operation mode and radar interrogation what assures the distant (few ten of meters) signal reading and real time monitoring.

INTRODUCTION

The dosimetry is one of the crucial techniques that are needed to assure personal safety and facilities security in the areas of high radioactivity (10kGy -10 MGy). Until now the existing solutions do not allow continuous monitoring as they cannot make in-situ readout. The objective of this study is to develop an innovative concept based on miniaturized HPD (Hydrogen Pressure Dosimetry) measurement technique coupled with passive RF pressure transducers that can be interrogated wirelessly with a radar (Figure 1) [START_REF] Jatlaoui | Working Principle Description of the Wireless Passive EM Transduction Pressure Sensor[END_REF]. The transducers principles are described in Figure 2. The irradiation will generate the outgassing of a polymer inside a micro-chamber. This overpressure leads to a deflection of a silicon membrane that modifies the resonant frequency of the RF resonator that can be wireless interrogated thanks to a radar. . Actually, the radiation levels achieved by now are small in comparison to the levels that will be achieved in the future when new generation of nuclear reactors (including nuclear fusion) or SLHC experiments will be launched. Such radiation doses cause the serious damage in the constructing materials and thus, the monitoring of the absorbed dose is crucial in order to assure the safe operation of these facilities.

According to the recent data, the typical 1GW NPP is monitored by ~ 10000 sensors, which are mostly: 2500 pressure transmitters, 150 gamma radiation monitors, 100 thermocouples and some 20 neutron detectors. The measurement equipment is basically divided into four classes: nuclear, process, radiation monitoring and special. The measured radiation doses are totally different depending on the place and purpose of measurement. Basically, five classes of radiation doses measurement are made : Personal monitoring (low doses), Coolant loop monitoring (normally low doses), Out-core infrastructure (normally low doses), In-core doses and reactor vessel structural materials (high doses), Accident measurement (extremely high doses).

The high energy physics facilities have totally different operating mode than the nuclear plants, but also the extreme high doses occurs during the experiments. We have to take into account that such facilities does not work in continuous mode what means that most of the time they are not operational. Nevertheless, the doses of radiation emitted during the colliding cycle are extreme and may last a short time. Again, monitoring of such doses is very often needed in the places where structural materials should be protected for the safety operation of the facility.

In Table 1 is listed the typical parameters that have to be addressed depending on the high doses application of the sensor.

Table1 : Dependency of main sensor parameters in a function of the application. The solution for high dose monitoring exists widely on the market since a long time. There are a lot of methods, where most of them rely on the variation of physic-chemical parameters of the different materials : Polymer Alanine Dosimetry (PAD) , Radiophotoluminescent Dosimetry (RPLD), Hydrogen Pressure Dosimetry (HPD), Thermo Luminescent Dosimetry (TLD) [START_REF] Brabcová | Dosimetry in Heavy Ion Beams Using Various Detectors[END_REF][START_REF] Ravotti | Radiation Monitoring in Mixed Environments at CERN: From the IRRAD6 Facility to the LHC Experiments[END_REF]. But the main problem with all these techniques is their readout. Even if now, we have all miniaturized equipment like spectrometers or other analyzers, we cannot make the in-situ readout. All of the sensors based on different crystals need to be taken from the measurement area and the read out must be done outof-site. The online monitoring is then impossible

Application

MATERIAL STUDY

The transducer is done using 4 different materials (Figure 2) : silicon as the structural material for movable membrane, metal for passive microwave elements (resonator, antenna), glass as the structural material for microwave component holder and finally the polymer that will outgas due to the radiolysis phenomena.

Polymer

The most important material that will be used for the sensor operation is a polymer that outgases due to the radiolysis phenomenon. Due to the ionizing radiation there are some effects on polymers that include: main chain rupture, side chain rupture, crosslinking, degradation, radiolysis, oxydation. The first four effects will affect mainly the mechanical properties of the polymer and the last one can be avoid using inert atmosphere. The radiolysis effect is generally characterized by the radiation yield G equals to number of gas molecules that are released by 1 g of polymer for the 100eV of energy absorbed due to radiation. This G factor is strongly dependant on polymer composition (0.1< G < 6) and decrease generally with dose [START_REF] Kinderman | Nuclear Radiation Effects on Structural Plastics and Adhesives, Astia Document[END_REF].

After studying publication data, we chosen High Density PolyEthylene (HDPE) as it is done in classical HPD dosimetry.

-PE has a quite high radiation yield (G between 4 to 6 depending on the composition) and G factor is energy independent what means that the amount of the released gas is dose dependent only. -PE releases mainly hydrogen what is caused by its simple structure. Then the outgassing effect should be not affected by the interaction between other gases and no hydrocarbons will be present inside the chamber. Generally hydrocarbons have non linear response. -Water absorption by PE is quite low (< 0.01% for HDPE) -PE has a high decomposition temperature (about 235 °C) that is interesting for the sensor fabrication process. -The neutron interaction should not affect a lot the radiolysis effect even if PE is a good neutron stopper (as it contains a lot of hydrogen).

-PE dissolves easily in many solvents. So it should be possible to dose it precisely by using the pipette technique, that could be useful in our case. -PE is the most common polymer with a simple composition and is then extremely cheap.

Glass

The main effect ofgamma irradiation on glass is the densification that modify the mechanical properties and can then induce additional stresses and may affect the silicon membrane deflection, the quality of silicon/glass bonding or cause cracks inside the glass. For 10 MGy gamma dose (maximum for targeted applications), Pyrex densification is around 0.04%. Most of the glasses used in microfabrication contains Boron. This element( 10 B) reacts very easily with the thermal neutrons and then induce also material densification (0.08% for 10 16 n/cm 2 thermal neutron). Concerning RF glass properties, we analyzed the shift of resonant frequency of a mmwave planar resonator used for the transducer and fabricated on glass substrate. The resonator (Figure 3) is obtained through 1µm thick aluminum layer deposited by PVD (Physical Vapor Deposition) on a 1mm thick Corning-Pyrex 7740 substrate. Resonators have been irradiated by high energy electron (6 MeV) with doses ranged from 11 kGy to 166 kGy and only slight variation of spectral response has been noticed. In a first approximation, the effect of irradiation on glass can be then neglected.

Silicon

Silicon, from the mechanical point of view, is highly resistant to gamma and neutron irradiation. The only effect of neutron irradiation on silicon is a transmutation to phosphorus that induces a slight variation of Young modulus.

TECHNOLOGICAL DEVELOPMENTS

The main critical step for the transducer fabrication concerns the silicon / glass anodic bonding with the HDPE. Indeed temperature around 450°C is generally used that is not compatible with HDPE properties. A test structure has been designed and fabricated (Figure 4) to validate the technology. High Density PolyEthylene (HDPE) is included inside a micro-chamber fabricated with silicon / glass technology. This cavity is connected to a silicon membrane. The H 2 gas generated by irradiation will then create an over pressure inside the cavity that deflects the silicon membrane and lead to the membrane breaking.

The fabrication process is the following. A mono-crystalline n-type double-side polished (100) silicon is first etched in KOH (80°C, 10M KOH) to fabricate HDPE cavity, flow-channel and membrane. Following this, first anodic bonding is performed on back side with glass substrate (Borofloat 33, Schott) in air at 450°C, 1000 V. After that, small portion (few mg) of HDPE is precisely placed inside the cavity and cured at 180°C in N 2 . Finally anodic bonding is done with top side glass (Borofloat 33, Schott) cover in pure nitrogen (1 bar) at 300°C, 1200 V. Test structures have been fabricated with different silicon membrane area (3mm x 3mm, 4mm x4 mm, 5mm x 5mm) and silicon membrane thickness around 30µm. Irradiations have been performed using high energy electron beam (6 MeV) with doses between 10kGy and 120kGy. The membrane breaking is obtained around 100kGy for 3mm x 3mm area, 80kGy for 4mm x 4mm area and 35kGy for 5mm x 5mm area. These results are consistent with analytical model of membrane breaking and validate two main points that are the hermetic sealing of the cavity at low temperature with HDPE and the strong outgasing of HDPE under irradiation [2]. where V HDPE is the HDPE volume (HDPE density equal to 0.95mg/mm 3 ).

For a typical HDPE volume (between 1 mm 3 and 10 mm 3 ), it is then possible to obtain under pressure around few bars for a large range of radiation dose (Table 2). The minimum V c volume is limited by the accuracy of the HDPE cavity filling (V c < V HDPE / 10). On the contrary, small V HDPE /V c ratio are not difficult to obtain as we have just to reduce the HDPE cavity and increase the micro-chamber surface. Pressure around few bars is useful for the full dose range as the pressure is high enough to reduce the mechanical parasitic effects and low enough to be compatible with simple micro-structure.

Test structure

Test structures have been designed in order to measure the silicon membrane deflection with an optical profiler after irradiation (Figure 5). Boss membrane has been chosen as it allows easier membrane deflection measurement and provides a quasi-uniform membrane displacement that will optimize the RF transducer sensitivity. The process used is similar to the process described in §3. The cavity depth in glass1 and the membrane thickness have been fixed respectively to 10µm and 50µm.

With the dimensions described previously, the volumes of HDPE micro-cavity (V HDPE ) and cavity chamber (V c ) are equal respectively to 8.8mm 3 and 0.7mm 3 . Figure 6 gives the expected under pressure P g generated for different ratio of HDPE cavity filling. The maximum volume of HDPE is a little bit lower than 8.8mm 3 (around 7.7mm 3 ) as the 50µm depth silicon channel reduce this potential volume. As it is also difficult to fill completely the HDPE cavity, a reasonable cavity filling between 25% and 50% is expected, that will give an overpressure between 0.3bar and 1bar for 100kGy. This overpressure leads to a membrane deflection between 1µm and 3µm that is high enough for easy measurements.

Irradiations on test structures are in progress and results will be given during the conference.

CONCLUSION

The first prototypes of miniature cavity with integrated polymer have been fabricated and the principle of hydrogen outgasing has been validated by silicon membrane breaking under over pressure. The ongoing works are focused on the polymer outgasing sensitivity using silicon membrane deflection. The next step will be the integration of hydrogen micro-cavity and MEMS pressure sensor.

1 SPECIFICATIONS

 1 Figure 1 : Sensor principle

Figure 3 :

 3 Figure 3 : Design of planar RF resonator

Figure 4 :

 4 Figure 4 : Schematics of first test structure

Figure 5 :

 5 Figure 5 : Test structure with boss silicon membrane Top : cross section drawing, Bottom : mask view

Figure 6 :

 6 Figure 6 : Expected under pressure P g generated for different ratio of HDPE cavity filling and for 100kGy

Sensor lifecycle Maximum Gamma dose Accuracy Aging effects Neutron level

	Normal NPP monitoring	2 years	20MGy	10%	Important	10 14 n/cm 2
	Accident NPP monitoring	Few days	10MGy	20%	Unimportant	0 (if no nuclear reaction)
	High energy physics facilities	10-15 years	Few MGy	10%	Important	0

Table 2 :

 2 Typical ratio volume for different dose range giving few bar under pressure

	Dose range	V HDPE / V c	P g (bars)
	10 kGy	10	1
	100 kGy	5	5
	1 MGy	0,5	5
	10 MGy	0,05	5

EWSHM 2014 -Nantes, France

ACKNOWLEDGEMENTS

This work was done in the framework of FP7-MNT-ERA.NET-DOSIMEMS project funded by EU via Foundation for Polish Science and by Midi-Pyrénées Region.


Rémy Vincent

Mikael Carmona

Olivier Michel

Jean-Louis Lacoume

PASSIVE ACOUSTIC SENSOR NETWORK LOCALIZATION; APPLICATION TO STRUCTURE GEOMETRY MONITORING

Keywords: Passive identification, Travel-time estimation, Sensor network localization, Geometry monitoring

In this article, passive acoustic sensor network localization is presented and applied to geometry structure monitoring. The methodology relies on passive travel-time estimation, which is here defined as the retrieval of an inter-sensor propagation delay, using uncontrolled ambient sources in an homogeneous acoustic propagation medium. Our approach relates to passive linear systems identification through the use of codas correlations to form estimators. We provide practical performances of such estimators and propose two approaches for geometry monitoring, both illustrated in the case of a steel beam.

INTRODUCTION

Localization of sensors within a network is a required step for monitoring the geometry of mechanical structures. This monitoring was historically achieved with informative stresses of the structures that were user-controlled and typically relies on Euclidean distances measurements between a controlled source and the sensors, e.g. GPS, optical devices, marker based techniques. In this paper, we instead work in passive context, which implies that no controlled sources are used. All sensors embedded in the investigated medium are excited by ambient uncontrolled sources only. In this context, ambient noise appears to be a useful source as long as it is sensed by all or some of the sensors. This approach relates to passive linear systems identification.

Pioneering work in passive identification was conducted in structural health monitoring [START_REF] Farrar | System Identification from Ambiant Vibration Measurements on a Bridge[END_REF], acoustic [START_REF] Weaver | Ward identities and the retrieval of Green's functions in the correlations of a diffuse field[END_REF] and seismology [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF]. In structural health monitoring, applications were released to determine structural parameters such as seismic velocity, modal damping and frequencies [START_REF] Farrar | System Identification from Ambiant Vibration Measurements on a Bridge[END_REF][START_REF] Prieto | Anelastic Earth structure from the coherency of the ambient seismic field[END_REF][START_REF] Snieder | Extracting the building response using seismic interferometry: Theory and application to the Millikan Library in Pasadena, California[END_REF]. Lately, it was even shown in [START_REF] Prieto | Impulse response of civil structures from ambient noise analysis[END_REF] that multi-story building impulse responses can be generated from ambient noise correlations only, from major earthquakes. In this article, the proposed method is different, for it is based on the passive geometry retrieval (shape, angle or curvature) using acoustic ambient noise instead of seismic one. Closest work is [START_REF] Sabra | Using Ocean Ambient Noise for Array Self-Localization and Self-Synchronization[END_REF] with underwater acoustic sensor network localization. The corner stone of the approach resides in accessing the distances between sensors, which we will tackle using indoors aerial acoustic pressure fields.

In [START_REF] Carmona | Ward identities for visco-acoustic and viscoelastic propagation media[END_REF], Green function retrieval was achieved in acoustics and a discussion was conducted on the damping model: in aerial acoustic, weak viscous damping was retained over constant damping, often applied with seismic wavefields [START_REF] Gouedard | Crosscorrelation of random fields: mathematical approach and applications[END_REF]. Our work is presently motivated by distance estimation performances assessment, which is divided in two parts. The first part is about detecting a change in the structure geometry: detection theory allows to set a threshold from the spatial resolution on the estimated inter-sensor Euclidean distances, which corresponds to a certain false alarm rate. The second part is about assessing structural parameters of the monitored structure, from the sensor network geometry, such as maximum curvature, maximum bending or Young modulus.

The article is organized as follows: in section 2, passive travel-time estimation theory is introduced. Performances of such estimators in aerial acoustics are derived and discussed. Experimental performances of distances estimation is conducted with a 7-microphones network and results are presented in section 3. This setup is extended to the two suggested approaches to passive monitoring: passive detection of changes in the sensor geometry and passive shape retrieval. The paper eventually ends summarizing pros and cons of the approach. As for the perspectives, suggestions are made to reach even better accuracy.

PASSIVE TRAVEL-TIME ESTIMATION

Passive/Active context

In an active context, travel-time estimation relies on a controlled transmitter and receiver. The signal transmitted u(t) propagates and is sensed at the receiver as r(t). r(t) is comonly designed as a delayed version of the emitted signal associated with an attenuated amplitude and corrupted by an additive measurement noise: r(t) = u(t -τ, λ ) + b(t), where τ is the delay, λ the attenuation and b the noise. The argmax of the correlation between r and u is the natural propagation delay estimator in the sense that it corresponds to the matched filter. Estimation theory then allows to derive a Cramer-Rao lower bound CRLB(τ) on the active travel time estimator, see for instance [START_REF] Ziv | Some Lower Bounds on Signal Parameter Estimation[END_REF].

In the passive context, sources are not controlled and source signals r(t) are consequently not available. The presented approach tackles the problem from a different angle: second order statistics are used to achieve the passive estimation of impulse responses between sensors in an homogeneous acoustic propagation medium, through the use of Ward identity [START_REF] Carmona | Ward identities for visco-acoustic and viscoelastic propagation media[END_REF]. Travel-times can then be extracted from those retrieved Green functions.

Propagation equations and Green functions

Studied media are homogeneous acoustic linear propagation media. Let f be some pressure source field and let s be the corresponding resulting pressure field. In this context, the propagation equation relating f and s has been derived by Stokes [START_REF] Royer | Elastic Waves in Solids[END_REF]:

∂ 2 ∂t 2 + η 2 ∇ ∂ ∂t + v 2 ∇ f (t, x) = s(t, x) ( 1
)
where η is the viscous-damping coefficient and v the wave celerity. As noticed and shown several times in the literature [START_REF] Carmona | Ward identities for visco-acoustic and viscoelastic propagation media[END_REF][START_REF] Gouedard | Crosscorrelation of random fields: mathematical approach and applications[END_REF][START_REF] Weaver | Ward identities and the retrieval of Green's functions in the correlations of a diffuse field[END_REF], dissipation has proven necessary in passive estimation theory. The Green function G(t, x| y) of a linear propagation medium is its impulse response. It stands for the field sensed at x at time t, when a pulse was triggered from y at time t = 0 s. The Green function relates the source and the corresponding response field as a filtering process:

f (t, p) = G(t -t ′ , p -x)s(t ′ , x)dt ′ d x = (G ⊗ T,S s)(t, p) (2)
where ⊗ T,S stands for the generalized convolution in both time and space. When the medium is homogeneous with a low damping coefficient and in the high-frequency regime [START_REF] Wright | New Directions in Linear Acoustics and Vibration, Quantum Chaos, Random Matrix theory and Complexity[END_REF], the Green function writes as a sum of delays with corresponding attenuated amplitudes:

G(t, p| x) = ∑ n≤0 a n δ (t -t n ) ( 3
)
where δ is the Dirac distribution, a n are the damping coefficients and t n are the delays. The first delay is the time of arrival:

t 0 = r v ( 4
)
where r is the Euclidean distance between p and x.

Green correlation

In the passive context, sources are unknown stochastic processes. While propagating, their wavefronts are shaped by the medium; filtering induces information on the medium which is sensed in the field.

To extract this information, the presented approach relies on Ward identity for which the framework deals with 2 nd order statistics of the sensed fields. Let γ f (τ, x, y) be the cross-correlation function of a stochastic and stationary field f at two locations x and y:

γ f (τ, x, y) = E [ f (u, x) f (u + τ, y)] (5)
Combining eq. ( 2) and ( 5) gives rise to a generalized interference formula, relating the auto-correlation function of field f generated by a source at s, to the signals sensed at two locations x and y:

γ f (τ, x, y) = C ⊗ T,S γ s (τ, x, y) (6)
where:

C = G ⊗ T,S G - (7)
and where G -(t, x| y) = G(-t, y| x) is the time reversed Green function. C is called Green correlation and in equation ( 6), it is the auto-correlation function of the pressure field generated by a spatiotemporal white source (i.e. γ s (τ, x, y) = δ (τ)δ ( xy)). Passive identification is then achievable provided sources are white in time and space, as in ( 6) information on the medium only remains.

In practice, ambient sources are used, which are barely white even within finite spectral bandwidth. However, seismologists [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF] have experimentally shown that the cross-correlations of codas of seismic events resulted in a good approximation of the cross-correlation of a propagated white noise.

This results is still true for acoustic propagation media. In Fig. 1, an instance of acoustic coda is displayed; it is a chunk of the field generated by an uncontrolled impulse source. Note that rigorous definition of the coda is difficult and that we thus base its extraction quite similarly to [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF], with a powerbased rule. Furthermore, a model of the coda was proposed in [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF] which allows to locate information in the phase of the coda. As a consequence, one-bit coda signals will be legitimately processed in the passive identification protocol. Note however that 1-bit correlations are distorded versions of the theoretical correlations. Corresponding correction is found in [START_REF] Pincinbono | Sur certains problèmes concernant la détection des signaux faibles[END_REF], see eq. ( 8):

γ s (τ) = sin π 2 γ s(τ ) (8)
where for signal s, γ s is its theoretical auto-correlation function and γ s its normalized 1-bit quantized version. In practice, it is clear that 1-bit quantification ease practical deployment in terms of data transmission, data storage and computation time. Note that 1-bit quantification furthermore preservs the argmaximum position, which is a fundamental property to our purpose.

Ward identity and travel-time retrieval

The Green correlation is now related to the Green function with the use of a Ward identity [START_REF] Weaver | Ward identities and the retrieval of Green's functions in the correlations of a diffuse field[END_REF]. In the acoustic framework and for a low viscous damping, the Ward identity writes [START_REF] Carmona | Ward identities for visco-acoustic and viscoelastic propagation media[END_REF]:

∂ 3 ∂t 3 γ G (t, x, y) = 1 η γ s (t, s) ⊗ T,S
Odd G(t, x, y) (

where the third time derivative and parameter η are dependent of the viscous-damping model. From Ward identity [START_REF] Prieto | Anelastic Earth structure from the coherency of the ambient seismic field[END_REF], one can observe that Green function can be retrieved from an estimation of the Green correlation. Because the air is weakly viscous, the notion of travel-time is preserved The red dot informs that an acoustic event is detected (hand-clap). The black dot specifies the start of the part of the signal that will be used for identification. The green dot ends the codas as the signal to noise ratio reaches some threshold. Thresholds on the power of the signal were used to process codas extractions. Each coda is then 1-bit quantized. Subcodas are formed from them, with a sliding window of given length and overlap. These sub-codas are used to perform estimation for each pair of sensors.

through the reconstruction of the Green function expressed in (3). As a consequence, an estimator of the travel-time will be the maximum argument of the estimated Green function. This is illustrated in Fig. 2, where a reconstructed Green correlation and retrieved Green function are displayed, between two microphones with sampling frequency 51kHz. In this example, sensors were 48cm away from each other. Amplitudes are normalized. 0.25 seconds of codas were used to compute 1-bit Green correlations and process identification.

The dot stands for the result of the passive travel-time estimation: 1.41ms, which leads to an estimated distance of 48.36cm.

PERFORMANCES AND APPLICATIONS TO STRUCTURAL GEOMETRY MONITORING

Resolution, convergence and variance

Performances of passive travel-time estimation will now be discussed in an air acoustic propagation medium. Sensors are omnidirectional microphones sampled at 51kHz.

Sound celerity is 343m/s when temperature is 20 • C. The spatial resolution thus is 6.2mm at 51kHz. It goes down to 6mm at 10 • C and up to 6.3mm at 35 • C. Note that the influence of altitude is even more negligible.

Codas are noises of interest. They are to be opposed to actual noises with respect to estimation; they can be described as parts of signals that are only seen by either one of the sensors. However, although cross-correlations of these latter are expected to interfer with travel-time extraction, it is only true at very small propagation times. The decay of their correlation function is indeed generally strong.

Acoustic codas are temporally white with a bandwidth that, in our case, is larger than the sampling rate. It is a basic setup in which, Shanon-Nyquist sampling requisite is hardly met as increasing sampling frequencies can raise technological issues. Note that this estimator is created from the sampling of the continuous estimator of the travel-time which has a bias. This latter is expected to be upper-bounded by the spatial resolution. If not, one can assess that misdetection occured. The estimator can be computed so that its variance reaches the order of magnitude of the spatial resolution. However, the numerical computation leads to a variance that has some properties:

• Correlations are approximated, with residues that have a variance that decreases proportionnaly to the number of points they were computed on; • The Green correlation is an average of N correlations of the type above, which leads to a corresponding variance that decreases proportionaly with N. • Eventually, the Green function is estimated from the 3rd-time differenciation of the Green correlation. One can show such numerical differenciation multiplies the variance by 20.

As shown in Fig. 2, we experimentally noticed that 250 averaged correlations to form the Green correlation, from typically 0.25sec of 1-bit codas, with sampling frequency f s = 51kHz, lead to 2 likely values for n (others are very unlikely events and will consequently be discarded). Because of this, it means that detection fails H 0 or succeeds H 1 , such that: The continous travel-time estimator has an actual variance that is of the same magnitude than the resolution ε. The actual digital estimator for the travel-time takes only a couple of discreet values. In the figure, the area of the shaded patch for instance stands for the probability that sample n is picked as travel-time estimate. Note that actual τ is shifted from the actual travel-time τ by a possible bias of the estimator µ.

H 0 : µ ≥ 1/ f s misdetection H 1 : µ ≤ 1/ f s detection and σ ≤ 1 6 f s (10)

Extension to a sensor network: from distances to positions

Let us now raise the reasoning to a sensor network of n sensors. For each pair of sensors, intersensor distance ought to be estimated. This leads to a maximum count of (n -1)(n -2)/2 intersensor distances accessible. These distances can be concatenated in a symmetric matrix D, where D i j is the estimated distance between sensor i and sensor j. Because these distances were estimated with performances described above, some additive estimation residue matrix B is used to express the estimation:

D = D(X) + B ( 11
)
where D(X) is the actual distances matrix, function of positions in matrix X. Also, note that from a practical point of view, it is realistic to account for the fact that distances may remain unestimated or badly done so. This is accounted for by using adjacency matrix A where A i j is a coefficient in between 0 and 1 which basically relates how confident the estimation is. It is for instance 0 if the distance couldn't be estimated and 1 if estimation protocol succeeded. Eventually, a model of the estimated distance can write using the Hadamar product ⊗:

D = A ⊗ (D(X) + B) (12)
In equation [START_REF] Sabra | Using Ocean Ambient Noise for Array Self-Localization and Self-Synchronization[END_REF], assessing the sensor positions X is an inverse problem. Various algorithms are available, one of them being Torgerson's algorithm [START_REF] Torgerson | Theory & Methods of Scaling[END_REF]. Note that in our case, the adjacency matrix A has coefficients that either are 1 (detection) or 0 (misdetection). When matrix A displays unrecovered distances, Kruskal's algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF] can be used to achieve positions estimations. Relative positioning of the sensor network can become asbolute with the use of anchors. Eventually, the sensor network geometry grants access to parameters such as curvatures, angles or bendings. Estimation errors in the distances will then have an impact on the estimated parameters, whose performances study is very particular to the applications.

We illustrate the argumentation with the example of a beam. The beam is 3m long, doubly supported at 1m from each extremities. It is made of steel, has a Young modulus E = 210e6N/m and quadratic moment I = 3.9 × 10 -8 m 4 . Fig. 4 displays the experimental setup. Sampling frequency is 51kHz with corresponding spatial resolution 6.2mm. Events are handclaps and codas are extracted as introduced in fig. 1.

        (13)
Practically speaking, acoustic codas successfully allow to retrieve distances at scales of a 'standard' room, from a centimeter to a few meters. In this study case, the beam dimensions allowed to retrieve all but 1 distances within a centimetric precision. It is not clear to us why the distance estimation between sensor 1 and 5 consistently failed.

Applications in bending detection and retrieval

Detection

In the first approach, extremities of the beam are loaded and the goal is to characterize the detectability of a change in the geometry of the beam. In this setup, sensors at each extremity of the beam are the ones that are most influenced by a change in the geometry from loads. As a consequence we will investigate a change in the distance between those two sensors and identify detectability of a bend. Fig. 5 displays the evolution of the equivalent change in the bending given 10 kg loads and a varying Young modulus. At a sampling rate of 51kHz, spatial resolution is 6.2mm. At given loads equal to 10kg, sensors would become this quantity closer when the Young modulus gets below 0.39E 0 . Corresponding maximum bending is located in the middle of the beam and is 4.3mm. A bend in this beam thus becomes detectable after this threshold. The shaded area in the figure is the area where a bend is not detectable as it would be subresolution. The next detectable bends have corresponding bending that are millimetric.

Evaluation

For the second approach, sensors are located in the medium with no control. The approach had now better rely on reconstructing the geometry of the beam. Fitting splines from retrieved sensor positions is an instance on how to achieve it in 2 dimensions. Fitting surfaces would be the extension in 3 dimensions. We now provide an example of geometry reconstruction, see fig. 6. The inverse problem that relates distances to positions proves to be very sensitive to the variance of distances estimation. Given a 6.2mm spatial resolution and a variance of the same magnitude, we display the result of meter-wise geometry retrieval. Variance in the positions estimation is centimetric, which prevents to use this configuration to assess the maximum bending of the beam from 4. This variance will however allow relevant applications at larger scales, e.g. for the monitoring of a vault. Any multi-dimensional scaling approach used for this inverse problem should be characterized through the evolution law that relates the variance in the positions σ 2 X to the variance in the distances σ 2 D . Although it may be described by very specific parameters, which depends on the study case, a general remark is that σ 2 X ≥ σ 2 D . Geometrical parameters to be estimated from the reconstructed geometry should then be of a much bigger order of magnitude in order for the estimation to be make sense.

CONCLUSIONS AND PERSPECTIVES

Performances of passive travel-time estimation have been presented in the air acoustic framework. Experiments have shown that the standard deviation of the estimator is lower than the sampling frequency. Thus, with a sampling at 51.2 kHz, accuracy is better than the centimeter. This performance allows to retrieve the geometry of a sensor network distributed at larger scales.

However, calculus derived in the context of an instrumented steel beam show that accuracy needs to be improved to detect changes and to retrieve the shape (or the bending) with an acceptable precision. Thus, two extensions are under investigation. The first one is related to the theoretical study of Green correlation retrieval with codas; even if the passive identification protocol works, all parameters are not entirely understood. Extraction of codas is probably the most probant one.

Secondly, experimentations with a higher sampling frequency will naturally improved the accuracy of the travel-time retrieval and will allow to discover the real limit of the estimation for the moment limited by technology.

Figure 1 :

 1 Figure 1 : Instance of acoustic codas

Figure 2 :

 2 Figure 2 : Retrieved Green function and travel-time

Figure 3 :

 3 Figure 3 : Travel-time estimator

Figure 4 :

 4 Figure 4 : Beam setup Below are distance estimation results for the 7 sensors. In (13), inter-sensor distances are displayed, corrupted by the estimation residue (in red). All lengths are expressed in centimeters. Passive travel-time estimation was processed 100 times and the displayed estimation residues result from the average of the discreet estimated values, weighted by their empirical occurance probability.         c1 c2 c3 c4 c5 c6 c7 c1 0 60+0.4 140-0.1 220-0.1 280-98.5 20+1.3 10+0.6 c2 0 80+1.6 160+0.2 220+0.5 40-0.6 50-0.3 c3 0 80+0.3 140+0.6 120+0.1 130+0.7 c4 0 60+0.7 200+0.4 210-1.1 c5 0 260-0.7 270-1.4 c6 0 10-0.1

Figure 5 :

 5 Figure 5 : bending and detection

Figure 6 :

 6 Figure 6 : bending and detection

EWSHM 2014 -Nantes, France


M.-A Torres-Arredondo

Julián Sierra-Pérez

Erion Zenuni

Guénaël Cabanes

José Rodellar

Alfredo Güemes

C.-P Fritzen

A PATTERN RECOGNITION APPROACH FOR DAMAGE DETECTION AND TEMPERATURE COMPENSATION IN ACOUSTO-ULTRASONICS

Keywords: Damage Detection, Acousto-Ultrasonics, Signal Processing, Pattern Recognition, Wavelet Transform, Principal Component Analysis, Self-Organizing Maps, Temperature Compensation

The global trends in the construction of modern structures require the integration of sensors together with data recording and analysis modules so that their integrity can be continuously monitored for safe-life, economic and ecological reasons. This process of measuring and analysing the data from a distributed sensor network all over a structural system in order to quantify its condition is known as structural health monitoring (SHM). Guided ultrasonic wave-based techniques are increasingly being adapted and used in several SHM systems which benefit from built-in transduction, large inspection ranges, and high sensitivity to small flaws. However, for reliable health monitoring, much information regarding the innate characteristics of the sources and their propagation is essential. Moreover, any SHM system which is expected to transition to field operation must take into account the influence of environmental and operational changes which cause modifications in the stiffness and damping of the structure and consequently modify its dynamic behaviour. On that account, special attention is paid in this paper to the development of an efficient SHM methodology where robust signal processing and pattern recognition techniques are integrated for the correct interpretation of complex ultrasonic waves within the context of damage detection and identification. The methodology is based on an acousto-ultrasonics technique where the discrete wavelet transform is evaluated for feature extraction and selection, linear principal component analysis for data-driven modelling and selforganizing maps for a two-level clustering under the principle of local density. At the end, the methodology is experimentally demonstrated and results show that all the damages were detectable and identifiable.

INTRODUCTION

The well-known properties of guided ultrasonic waves have led to a burst of studies using Lamb waves for the detection and analysis of defects in structures. Nevertheless, in reality, structures are susceptible to varying environmental and operational conditions which affect the measured signals.

It is well known that temperature as well as damage can have similar effects on the dynamic behaviour of a structure [START_REF] Sohn | Effects of environmental and operational variability on structural health monitoring[END_REF]. These environmental and operational changes of the system can often mask slight changes in the structural dynamic responses caused by damage and covert the damage detection step into a very complex process. Therefore, it is very important to understand the impact of these changing conditions and take them into account. For example, Fritzen et al. [START_REF] Fritzen | Elimination of temperature effects on damage detection within a smart structure concept[END_REF] modified an existing subspace-based identification method with temperature compensation for damage diagnosis. Moll et al. [START_REF] Moll | Compensation of Environmental Influences for Damage Detection using Classification Techniques[END_REF] studied the compensation of environmental influences using a simulation model and a laboratory structure. Buethe et al. [START_REF] Buethe | Applications of Self-Organizing Maps in Structural Health Monitoring[END_REF] evaluated self-organizing maps in order to distinguish between environmental changes and damage of within a structure so that false alarm minimization could be accomplished under changing environmental conditions. Torres and Fritzen [START_REF] Torres Arredondo | Ultrasonic guided wave dispersive characteristics in composite structures under variable temperature and operational conditions[END_REF] presented theoretical developments, numerical and experimental results in order to analyze the effects of all the aforementioned sources of variability on wave propagation velocities, directionality and attenuation. This contribution extends previous works from the authors in the field of damage detection and classification [START_REF] Torres Arredondo | Data-Driven Multivariate Algorithms for Damage Detection and Identification: Evaluation and Comparison[END_REF][START_REF] Torres Arredondo | Acoustic Emission Testing and Acousto-Ultrasonics for Structural Health Monitoring[END_REF]. On that account, this paper is concerned with the experimental validation of a data-driven structural health monitoring methodology where a damage detection and classification scheme based on an acousto-ultrasonic (AU) approach is followed. The layout of this paper is as follows. In the first part, an analysis of the effects of variable temperature and operation conditions on wave propagation is shown. The second part introduces the required background for understanding the proposed algorithms and methodology for damage detection. The third part presents the analysis together with the discussion of the results. Finally, concluding remarks are given in the last part.

GUIDED WAVE SENSITIVITY ANALYSIS

Several studies regarding the effects of temperature variability on the measured dynamic responses of structures have shown that temperature variation may change the material properties of a structure [START_REF] Qiu | Modeling for Suppression of Moisture/Temperature Induced Dimensional Changes in Fibrous Composite Structures[END_REF]. Additionally, factors such as material age effect, moisture content and structure operation affect significantly the wave propagation characteristics in the material. In order to depict the changes introduced by these effects, a CFRP plate made of 4 equal layers with a total thickness of 1.7mm and stacking of [0° 90° 90° 0°] was studied. Temperature tests were conducted in a temperature-controlled oven. During the test runs the temperature was raised stepwise up to T=60 ± 2°C. Figure 1(a) shows the structure with dimensions 200mm×250mm. Nine piezoelectric transducers PIC151 from PI Ceramics were attached to the surface of the structure. The excitation voltage signal (at P 5 ) was a 12V Hann-windowed toneburst with a carrier frequency of 30kHz with 5 cycles. Figure 1(b) shows that the shape of the signal changes significantly due to the wet surface influence. The influence of the variation in temperature also causes a change of structural dynamics. The dynamic response signal of sensor number two and three decreased monotonically in peak-topeak magnitude with increasing temperature (Figure 1(c)). Interestingly, the inverse effect was seen in sensor number one. The reason of these time-shifts is both thermal expansion and changes in wave velocities with temperature. The attenuation of Lamb wave can be regarded to both wave dispersion as a result of frequency dependent phase velocities and attenuation loss due to frequency/temperature dependent material damping.

ALGORITHMS AND METHODOLOGY FOR DAMAGE DETECTION

The system is based on a distributed array of permanently attached piezoelectric transducers where pairs of transducers are used in pitch-catch configuration. The dynamic responses collected from an actuation step are stored and then pre-processed by the discrete wavelet transform (DWT), as a

(a) (b) (c)
feature extraction technique, in order to calculate coefficients representing valuable time and frequency information from these responses. This procedure is repeated for all actuation steps and for all environmental conditions over a temperature range of these changing conditions. The DWT analysis is performed by means of a fast, pyramidal algorithm related to a two-channel subband coding scheme using a special class of filters called quadrature mirror filters as proposed by Mallat (see Figure 2) [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. The optimum number of level decompositions is determined based on a minimum-entropy decomposition algorithm. The family of Daubechies wavelets ('db8') was chosen for this study. The approximations coefficients are taken here since they represent the interesting dynamics of the recorded waveforms and the detail coefficients will be considered as highfrequency noise. Afterwards, the coefficients for each actuation step are fused following unfolding procedures (multiway) as it is done in multivariate statistical procedures for monitoring the progress of batch processes. This sensor data fusion procedure is depicted in Figure 3. Finally, the fused coefficients are used to build a data-driven model based on Principal Component Analysis (PCA) [START_REF] Tibaduiza | A Study of Two Unsupervised Data Driven Statistical Methodologies for Detecting and Classifying Damages in Structural Health Monitoring[END_REF]. When new structural responses are available (from an unknown state), DWT approximation coefficients are extracted and projected into the respective models by each actuation step. Finally, by retaining a certain number of principal components, squared prediction error (SPE) measures for all the actuation steps are calculated and used as input feature vectors for a self-organizing map (SOM) for the detection and identification tasks [START_REF] Torres Arredondo | Multivariate data-driven modelling and pattern recognition for damage detection and identification for acoustic emission and acousto-ultrasonics[END_REF]. In order to calculate the number of clusters inside the data and provide a way for damage classification, a simultaneous two-level clustering approach using SOM is evaluated [START_REF] Cabanes | Learning the Number of Clusters in Self Organizing Map, in Self-Organizing Maps[END_REF]. Within this approach the structure of the data and its segmentation is learnt at the same time by using both distance and density information. The clustering algorithm assumes that a cluster is a dense region of objects surrounded by a region of low density. This approach is very effective when the clusters are irregular or intertwined, and when noise and outliers are present. The proposed clustering algorithm, called the DS2L-SOM algorithm (Densitybased Simultaneous Two-Level -SOM) divides automatically a given dataset into a collection of subsets (clusters representing the pristine and damaged structural states) and the number of clusters is determined automatically during the learning process, i.e., no a priori hypothesis for the number of clusters is required [START_REF] Cabanes | Enriched topological learning for cluster detection and visualization[END_REF].

Figure 3. Unfolding procedure for sensor data fusion.

EXPERIMENTAL SETUP AND RESULTS

A fully independent experiment was performed in order to evaluate the practical performance of the proposed methodology. The experiment used pairs of transducers operating in pitch-catch mode.

The input signal to the actuators was generated using the arbitrary signal generation capability of a combined signal generator and oscilloscope instrument manufactured by TiePie Engineering, Holland. The time histories were digitized at a sampling frequency of 50 MHz and transferred to a portable PC for post-processing. The structure is a simplified aircraft composite skin panel made of carbon fibre reinforced plastic (CFRP) depicted in Figure 4.

The overall size of the plate is approximately 500×500×1.9mm and its weight is about 1.125 kg. The stringers are 36mm high and 2.5mm thick. The plate and the stringers consist of 9 plies.

Damage on the multilayered composite plate was simulated by placing magnets with different masses at the same position (see Figure 4) on the structure as artificial damage in four increasing steps of 0.024 kg. Temperature was varied from 35°C to 75°C in steps of 10°C for a total of five temperature levels. The excitation voltage signal is a 12V Hanning windowed cosine train signal with 5 cycles and carrier frequency of 50kHz. In a first step, all the collected baseline signals were processed by means of the DWT (a 8th level was found to be the optimum in all experiments), fused and presented to the DS2L-SOM algorithm in order to depict the changes introduced by the changing temperature. To find the optimal map size, a control run is repeated by hanging the map size for all the experiments here. As it can be seen from Figure 5(a) and 5(b), all the temperature levels were properly separated and distinguished. The U-Matrix, showing the average distance of a cell to its neighboring cells, allowed depicting the difference between the different formed clusters. Additionally, the baseline models are built with all the data coming from the different temperature conditions. A review of the variances retained in the components was performed in order to define the optimal number of components required for building the PCA models from the pristine structural condition. It was found that the first three hundred components included around 95% of the total variance into the reduced models for each actuation step. As it is depicted in Figure 6(a) and (b), all the damaged states were not only separated from the pristine condition but also between them. Nevertheless, due to the big amount of data used for building the baseline models, a huge number of retained components (almost all) are required for obtaining an accurate approximate model. In order to depict how this effect could be mitigated, new models are built just with baseline measurements lying in the range where damage records were taken, i.e. around 35°C. In this case, it was found that the first thirty components included around 95% of the total variance into the reduced models for each actuation step. Figure 7(a) and (b) show, that here again, the proposed methodology was capable of separating all the damage cases from the pristine condition with the additional advantage of using a much more compact and reduce model in comparison to the previous example. These results pose the requirement of improving the proposed methodology by following an approach for optimal baseline selection and compensation of the recorded signals.

CONCLUSIONS

The goal of this study was to illustrate and develop a methodology to counterbalance the effects of environmental sources of variability on the performance of baseline data-driven models within the context of structural health monitoring systems and damage detection algorithms. The proposed methodology proved being effective in the detecting and classifying the different structural states under different temperature conditions. The approach includes the combination of Discrete Wavelet Transform, Multi-Way Principal Component Analysis, Squared Prediction Error measures and Self-Organizing Maps. Future work will involve the improvement of the proposed methodology by including more advanced techniques for temperature compensation.

Figure 1 .

 1 Figure 1. EOC Influences: (a) Experimental Setup, (b) Wet Surface Effect, (c) Temperature Effects on Wave Directionality.

Figure 2 .

 2 Figure 2. Discrete wavelet transform decomposition tree.

Figure 4 .

 4 Figure 4. Experimental Setup: Stiffened Composite Panel.

Figure 5 .

 5 Figure 5. All baselines presented to the DS2L-SOM algorithm: (a) U-Matrix and (b) DS2L-SOM cluster map.This result is of special attention in baseline-based methods where the detection and characterization of damage is performed normally by means of metric indices by comparison of two dynamic response signatures what can lead to trigger a false alarm for damage detection just because of a change in environmental conditions. For the second experiment, four damage evolution steps were simulated using a local mass increase as it was defined before (all damages at 35°C).

Figure 6 .

 6 Figure 6. Damage detection with baseline models built with all the temperature range (35°C-75°C): (a) U-Matrix and (b) DS2L-SOM cluster map.

Figure 7 .

 7 Figure 7. Damage detection with baseline models built with the same temperature level (only 35°C): (a) U-Matrix and (b) DS2L-SOM cluster map.

EWSHM 2014 -Nantes, France

ACKNOWLEDGEMENT

The authors would like to express their gratitude to the Research School on Multi Modal Sensor Systems for Environmental Exploration (MOSES) and the Centre for Sensor Systems (ZESS) for sponsoring the research presented herein. Additionally, the authors would like to thank Mr. Raghavendra Singh from the University of Siegen for his help during the experiments.