Fabien Treyssède
email: fabien.treyssede@lcpc.fr

Laurent Laguerre

Investigation of elastic modes propagating in multi-wire helical waveguides

Keywords: waveguide, propagation, helical, strand, finite element, experiment

Elastic guided waves have some potential for non-destructive inspection of civil engineering multi-wire steel cables. However, wave propagation inside such structures is not yet fully understood. This paper investigates multi-wire helical waveguides with special attention to the common seven-wire strand configuration (one straight core surrounded by one layer of six helical wires).

A helical coordinate system is first proposed. Though non-orthogonal, this system preserves translational invariance along the helix centreline to explicitly perform a spatial Fourier transform. Then, it is shown that for the analysis of multi-wire helical strands a twisting system -which is a special case of helical systems -is translationally invariant. A semi-analytical finite element method is developed reducing the problem on the cross-section only.

A straightforward computation of energy velocity is proposed. Dispersion curves for a single straight wire and a helical wire are first computed to verify the adequacy of the twisting system. Finally the seven-wire strand is analysed using simplified contact conditions. Theoretical dispersion curves are compared to low-frequency magnetostrictive measurements.

Introduction

Elastic guided waves are attractive tools for non-destructive testing and structural health monitoring because of their potential ability to inspect, in a single measurement, the whole width of the structure over an appreciable length along the guiding direction. Theoretical modelling effort is generally required to enhance physical insight into the complex behaviour of guided waves, which are multimodal and dispersive. Though cylindrical elastic waveguides have extensively been studied since the early works of Pochhammer in 1876 and Chree in 1889 [START_REF] Meeker | Guided wave propagation in elongated cylinders and plates[END_REF][START_REF] Zemanek | An experimental and theoretical investigation of elastic wave propagation in a cylinder[END_REF], less attention has been paid to helical and multi-wire waveguides. A typical application in civil engineering concerns steel multi-wire cables, which are widely employed in bridges (suspended, cable-stayed and prestressed) and anchored retaining wall constructions. These civil engineering structures suffer from ageing and degradation due to corrosion and fatigue.

The basic element of these cables is usually a seven-wire strand, consisting of one straight cylindrical core wire surrounded by one layer of six helical wires. The simplest and usual way to consider wave propagation in a sevenwire strand is to approximate the latter by an infinite cylinder of equivalent diameter. Various experimental studies have been conducted in the literature for the cylindrical bar and the seven wire strand, trying to find similarities between both behaviours that could be interpreted from Pochhammer-Chree solutions, either for the bare waveguide case [START_REF] Kwun | Effects of tensile loading on the properties of elastic-wave propagation in a strand[END_REF][START_REF] Laguerre | Low-frequency ultrasound reflectometry device based on magnetoelastic transducers for the non destructive evaluation of steel rods and cables[END_REF][START_REF] Laguerre | Magnetostrictive pulse-echo device for non-destructive evaluation of steel cylindrical materials using guide waves[END_REF][START_REF] Rizzo | Load measurement and health monitoring in cable stays via guided wave magnetostrictive ultrasonics[END_REF][START_REF] Rizzo | Wave propagation in multi-wire strands by wavelet-based laser ultrasound[END_REF] or the embedded one [START_REF] Pavlakovic | Inspection of tendons in post-tensioned concrete using guided ultrasonic waves[END_REF][START_REF] Beard | Ultrasonic guided waves for inspection of grouted tendons and bolts[END_REF]. Even if the non-destructive inspection strategy of a seven-wire strand with guided ultrasonic waves has gained from this approach, the authors point out the limits of Pochhammer-Chree solutions for the accurate interpretation of experimental data. Indeed, the theoretical understanding of guided ultrasonic waves in multi-wire strands is still challenging because of the complexity of this structure, due to the helical geometry of peripheral wires, the inter-wire coupling and contact effects, the presence of applied loads and concrete embedment (if any).

The aim of this paper is to propose a numerical method to study elastic wave propagation along multi-wire helical waveguides in order to help the interpretation of guided waves in these structures. According to the complexity of the problem, a computational approach is preferred to a purely mathematical approach. As a first step, inter-wire contact conditions will be simplified, and prestress and embedment will not be considered.

In order to deal with complex geometries, some of the most popular and efficient numerical techniques involve finite element (FE) methods. The socalled semi-analytical finite element (SAFE) method is a first popular approach to study uniform straight waveguides of arbitrary cross-section -see for instance Gavric [START_REF] Gavric | Computation of propagative waves in free rail using a finite element technique[END_REF], Hayashi et al. [START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF], Damljanovic and Weaver [START_REF] Damljanovic | Propagating and evanescent elastic waves in cylindrical waveguides of arbitrary cross-section[END_REF], Bartoli et al. [13]. This method assumes an axial dependence of the form e iks (k and s are the wavenumber and distance along the waveguide axis respectively), hence reducing the problem from three to two dimensions (only the cross-section is meshed). The SAFE method has also been extended to curved waveguides. Demma et al. [START_REF] Demma | The effect of bends on the propagation of guided waves in pipes[END_REF], Finnveden and Fraggstedt [START_REF] Finneveden | Waveguide finite element for curved structures[END_REF] investigated toroidal waveguides. Onipede and Dong [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF] extended SAFE methods to study uniformly pretwisted waveguides along a straight axis.

A second approach is based on the theory of wave propagation in periodic structures from Floquet's principle (this approach is more general because non-uniform waveguides can be analysed). A review can be found in Ref. [START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from southampton, 1964-1995[END_REF]. Based on a general theory presented by Mead [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF], some periodic FE approaches and procedures have then been developed -see for instance Gry and Gontier [START_REF] Gry | Dynamic modeling of railway track: a periodic model based on a generalized beam formulation[END_REF], Mace et al. [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF]. The periodic FE method allows to study the single repetitive substructure alone, thanks to the application of a set of periodic boundary conditions involving a propagation constant corresponding to the eigenvalue. Periodic FE methods only need the mesh of one repetitive cell.

Recently, Treyssède extended both periodic FE [START_REF] Treyssède | Numerical investigation of elastic modes of propagation in helical waveguides[END_REF] and SAFE [START_REF] Treyssède | Elastic waves in helical waveguides[END_REF] methods to model a single helical wire. Both methods are based on a specific curvilinear coordinate system attached to the helical waveguide. Such a system is non orthogonal but remains translationally invariant along the helix centreline, which allows a Fourier analysis (or a Floquet analysis with an arbitrarily small repetitive cell). In this paper, it is shown how the proposed approach can also be readily used to study a multi-wire helical waveguide. In Sec. 2, the helical coordinate system is first recalled. A particular twisting coordinate system is then discussed for the analysis of multi-wire waveguides. A SAFE method, whose weak variational formulation is rewritten in terms of the helical coordinate system, is then presented in Sec. 3. Section 4 gives preliminary results for a single wire both for the cylindrical and helical cases. Section 5 shows some numerical results (including dispersion curves and modeshapes) for the seven-wire strand configuration in civil engineering. From longitudinally-polarized magnetostrictive measurements in the low-frequency range (<500kHz), a first experimental validation of theoretical dispersion curves is performed both for the single constitutive wires (cylindrical and helical) as well as the assembled seven-wire strand.

Translationally invariant coordinate systems

Let (x, y, s) denote a curvilinear coordinate system attached to a curved waveguide, where x and y are the cross-section coordinates and s is the axial coordinate. In order to still speak of propagation modes, the physical system must be translationally invariant. In other words, one has to assume that an exponential axial dependence e iks can be separated from all field components (which is equivalent to perform a spatial Fourier analysis in the s direction).

This assumption indeed requires to meet the following three conditions:

1. the cross-section of the waveguide does not vary along s; 2. the material properties do not vary along s;

3. the coordinate system (x, y, s) is such that s does not appear explicitly in the coefficients of the equilibrium equations.

The first condition is an obvious geometrical condition. Condition 2 will be assumed to be satisfied throughout this work. The aim of this section is to propose some coordinate systems verifying the third non-trivial condition.

For a fundamental introduction to the use of general curvilinear coordinate systems, the reader may refer to Chapter 2 of Ref. [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF] for instance.

Helical coordinate system

A helical coordinate system is now built. One starts by defining the helix centreline curve, described by the following position vector in the Cartesian orthonormal basis (e X , e Y , e Z ):

R(s) = R cos( 2π l s + θ)e X + R sin( 2π l s + θ)e Y + L l se Z (1)
where l = √ L 2 + 4π 2 R 2 is the curvilinear length of one helix step. R and L are respectively the radius of the centreline in the (X, Y ) Cartesian plane and the helix step along the Z axis (see Fig. 1). θ is the helix phase angle in the Z = 0 plane. For instance, a seven-wire strand comprises six helical waveguides with θ = (k -1)π/3, k = 1, ..., 6. The unit tangent, normal and binormal vectors to the centreline are respectively obtained from T = dR/ds and the Serret-Frenet formulae [START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathematica[END_REF]: dT/ds = -κN, dN/ds = τ B + κT and dB/ds = -τ N (B = T ∧ N). Note that N is oriented outward the curvature in this paper (this is just a matter of taste). For a helix, both the curvature κ = 4π 2 R/l 2 and torsion τ = 2πL/l 2 are constant. In the Cartesian basis, N, B and T are expressed as:

N(s) = cos( 2π l s + θ)e X + sin( 2π l s + θ)e Y , B(s) = -L l sin( 2π l s + θ)e X + L l cos( 2π l s + θ)e Y -2πR l e Z , T(s) = -2πR l sin( 2π l s + θ)e X + 2πR l cos( 2π l s + θ)e Y + L l e Z (2)
In this way, a new coordinate system can be built from the orthonormal basis (N, B, T), for which any position vector Φ = Xe X + Y e Y + Ze Z is expressed as follows: yielding the non-orthogonal covariant basis (∂Φ/∂x, ∂Φ/∂y, ∂Φ/∂s) denoted by (g 1 , g 2 , g 3 ):

Φ(x, y, s) = R(s) + xN(s) + yB(s) (3)
g 1 = N(s), g 2 = B(s), g 3 = -τ yN(s) + τ xB(s) + (1 + κx)T(s). (4)
The covariant metric tensor, defined by (g) ij = g i • g j is:

g =      1 0 -τ y 0 1 τ x -τ y τ x τ 2 (x 2 + y 2 ) + (1 + κx) 2      . ( 5
)
g does not depend on s. As a consequence, the coefficients of any partial differential operators expressed in the so-defined coordinate system are independent on s too (condition 3 is hence satisfied).

For clarity the contravariant basis (g 1 , g 2 , g 3 ), defined by g i • g j = δ i j , is given by:

g 1 = N(s) + τ y 1+κx T(s), g 2 = B(s) -τ x 1+κx T(s), g 3 = 1 1+κx T(s) (6)
yielding the following contravariant metric tensor, defined by g mn = g m • g n :

G = g -1 = 1 g      g + (τ y) 2 -τ 2 xy τ y -τ 2 xy g + (τ x) 2 -τ x τ y -τ x 1      (7)
where g = (1 + κx) 2 is the determinant of g. The Christoffel symbol of the second kind Γ k ij , defined by Γ k ij = g i,j • g k , can be readily calculated from the Serret-Frenet formulae:

Γ k 11 = Γ k 12 = Γ k 21 = Γ k 22 = 0, Γ 2 23 = Γ 2 32 = Γ 3 23 = Γ 3 32 = 0, Γ 1 33 = -κτ y 2 1+κx -κ(1 + κx) -τ 2 x, Γ 2 33 = κτ 2 xy 1+κx -τ 2 y, Γ 3 33 = -κτ y 1+κx , Γ 1 13 = Γ 1 31 = κτ y 1+κx , Γ 1 23 = Γ 1 32 = -τ, Γ 2 13 = Γ 2 31 = -κτ x 1+κx + τ, Γ 3 13 = Γ 3 31 = κ 1+κx . (8)
2.2. Particular system for a multi-wire waveguide

Consider wires of circular cross-section in the normal plane to the guiding axis. The cross-section of a single helical wire does not change along the helical axis s (condition 1 is also fulfilled). Hence for helical waveguides, the previously described coordinate system allows a wave mode analysis. The analysis of a single straight wire, corresponding to the special case κ = τ = 0, is also readily performed. However, a question arises about the choice of the invariant coordinate system to be used for a multi-wire helical waveguide, made of both straight and helical wires. Of course, the choice κ = τ = 0 cannot be applied because condition 1 would not be satisfied for helical wires, and similarly κ = 0, τ = 0 cannot be applied because it would not be satisfied for straight wires. The same problem also occurs for any additional layer of peripheral wires because each layer would have a different curvature and torsion.

The adequate system is indeed given by κ = 0 and τ = 2π/L. It corresponds to a twisting coordinate system along the Z axis (s ≡ Z) with axial periodicity L. The (x,y) plane rotates around Z but remains parallel to (X,Y ). With this choice, a central straight wire (cylinder) has an invariant circular cross-section along Z ("a twisted cylinder remains a cylinder").

Furthermore, the cross-sections of peripheral helical wires do not change ei-ther along Z. In Sec. 4, it will be checked that this coordinate system yields right results for both a cylindrical and a helical waveguides. It should also be noted that this system coincides with the one proposed in Ref. [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF] for the analysis of pretwisted waveguides. It also had already been considered

in electromagnetics [START_REF] Nicolet | Modeling of twisted optical waveguides with edge elements[END_REF][START_REF] Nicolet | Finite element analysis of helicoidal waveguides[END_REF].

However, one has to be careful with such a system because the crosssection shape of helical wires is no longer circular (the cutting plane being not normal to the helix centreline). Let us find the cross-section cut by the plane Z =0. Equation (3) yields:

         X = (R + x) cos( 2π l s + θ) -y L l sin( 2π l s + θ) Y = (R + x) sin( 2π l s + θ) + y L l cos( 2π l s + θ) Z = L l s -2πR l y (9)
so that s = 2πyR/L at Z =0. In a helical coordinate system, the surface of a helical wire with a circular cross-section of radius a is given by (x, y) = (a cos t, a sin t), with t ∈ [0; 2π] and s ∈ R. Consequently, the cross-section shape can be parametrized as follows:

   X(t) = (R + a cos t) cos(αa sin t + θ) -a L l sin t sin(αa sin t + θ) Y (t) = (R + a cos t) sin(αa sin t + θ) + a L l sin t cos(αa sin t + θ) (10)
where α = 4π 2 R/(lL). With a twisting coordinate system, this shape must be used for generating the FE mesh (see Sec. 4).

As a side remark, the particular case τ = 0 degenerates into a coordinate system attached to a toroidal waveguide of curvature κ (of no interest in this paper).

Numerical method

SAFE formulation

One assumes a linearly elastic material, small strains and displacements with a time harmonic e -iωt dependence. There is no external force for the purpose of studying propagation modes. The 3D variational formulation governing elastodynamics is given by:

Ω δǫ T σdV -ω 2 Ω ρδu T udV = 0 ( 11
)
for any kinematically admissible trial displacement field δu = [δu n δu b δu t ] T .

Components in the orthonormal Serret-Frenet basis (N, B, T) are preferred here. δǫ denotes the virtual strain vector [δǫ nn δǫ bb δǫ tt 2δǫ nb 2δǫ nt 2δǫ bt ] T and σ is the stress vector [σ nn σ bb σ tt σ nb σ nt σ bt ] T . The superscript T denotes the matrix transpose. ρ is the material density and Ω is the structural volume.

The stress-strain relationship is σ = Cǫ, where C is the matrix of material properties (whose components are also given with respect to the Serret-Frenet basis). The volume element dV is given by dV = √ gdxdyds, where g is the determinant of the metric tensor previously defined. √ g can be understood

as the Jacobian of the transformation.

The strain-displacement relation can be written as follows:

ǫ = (L xy + L s ∂/∂s)u (12)
where L xy is the operator containing all terms except the derivatives with respect to the s-axis. Now an exponential e iks is assumed for u (e -iks for δu) and separated from all field components, k being the axial wavenumber.

∂/∂s can be replaced by ±ik. The problem is hence reduced from three to two dimensions (from the volume Ω to the cross-section S of the waveguide). Then, the FE discretization of Eq. ( 11) finally leads to the following eigenvalue problem for the column vector U containing nodal displacements:

{K 1 -ω 2 M + ik(K 2 -K T 2 ) + k 2 K 3 }U = 0 (13)
with the following elementary matrices:

K e 1 = S e N eT L T xy CL xy N e √ gdS, K e 2 = S e N eT L T xy CL s N e √ gdS, K e 3 = S e N eT L T s CL s N e √ gdS, M e = S e ρN eT N e √ gdS (14)
where N e is a matrix of nodal interpolating functions of displacement on the element and dS = dxdy.

The solution of Eq. ( 13) yields the propagation modes. Because of the symmetry of K 1 , K 3 and M and using the property detA T = detA (A is any matrix), it can easily be checked that if k is an eigenvalue of ( 13), then -k is also an eigenvalue. Hence, the eigenproblem has two sets of eigensolutions (k j , U + j ) and (-k j , U - j ) (j = 1, ..., n), representing n positivegoing and n negative-going wave types (n being the number of degrees of freedom (dofs)). With no damping, eigensolutions for which k j is purely real, purely imaginary and fully complex represent propagating waves, evanescent waves and inhomogeneous waves (decaying but oscillatory) respectively. At fixed real k, the eigenproblem ( 13) is linear for finding ω 2 . This simpler approach is useful only in the case of propagating modes in undamped systems. Given ω and finding k, the eigenproblem is quadratic. It can be recast into a generalized linear eigensystem written for [U T kU T ] T in order to be solved by standard numerical solvers -see Treyssède [START_REF] Treyssède | Elastic waves in helical waveguides[END_REF] for instance or Tisseur and Meerbergen [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF] for more details.

Expression of operators in curvilinear coordinate system

The strain-displacement relationship ǫ = 1/2(∇u + ∇u T ) (here, ǫ must be understood as a second order tensor), must be written in a general nonorthogonal curvilinear coordinate system. Using product and differentiation rules of such a system (see for instance Chapter 2 of Ref. [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF]), it can be shown that the strain tensor is given by:

ǫ ij = 1/2(u i,j + u j,i ) -Γ k ij u k . ( 15
)
Subscripts i = 1, 2, 3 denote covariant components with respect to the contravariant basis (g 1 , g 2 , g 3 ). The notation (•) ,i (i = 1, 2, 3) is used for derivatives with respect to x, y and s respectively.

Note that covariant components generally have non-physical units and the contravariant basis is non-orthogonal. Physical components with respect to the orthonormal Serret-Frenet basis (N, B, T) are preferred. Then we have the following transformations:

u i = J αi u α , ǫ ij = J αi ǫ αβ J βj (16)
with:

J =      1 0 -τ y 0 1 τ x 0 0 1 + κx      (17)
where the greek subscripts α = n, b, t denote components with respect to (N, B, T). Note that g = J T J.

Finally from Eqs. ( 15)-( 17) and after calculations, the following operators L xy and L s involved in Eq. ( 12) are:

L xy = 1 1 + κx               (1 + κx)∂/∂x 0 0 0 (1 + κx)∂/∂y 0 κ 0 τ y∂/∂x -τ x∂/∂y (1 + κx)∂/∂y (1 + κx)∂/∂x 0 τ y∂/∂x -τ x∂/∂y -τ -κ + (1 + κx)∂/∂x τ τ y∂/∂x -τ x∂/∂y (1 + κx)∂/∂y               , (18)
and

L s = 1 1 + κx               0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0               . ( 19
)
The components of C must also be expressed with respect to the Serret-Frenet basis. For an isotropic material, one has:

C = E (1 + ν)(1 -2ν)               1 -ν ν ν 0 0 0 ν 1 -ν ν 0 0 0 ν ν 1 -ν 0 0 0 0 0 0 (1 -2ν)/2 0 0 0 0 0 0 (1 -2ν)/2 0 0 0 0 0 0 (1 -2ν)/2               ( 20
)
where E is the Young modulus and ν denotes the Poisson coefficient.

Computation of energy velocity

The energy velocity is an important wave property that remains appropriate even for damped media [START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF] (as opposed to the group velocity definition v g = ∂ω/∂k, which is not generally valid for damped waveguides). From the previously defined SAFE matrices, a straightforward computation of energy velocity is proposed in this subsection.

The cross-section and time averaged energy velocity in waveguides is defined as follows [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]:

v e = S P • ndS S ( Ēk + Ēp )dS (21)
where bars denote time averaging. n is the unit vector along the propagation direction (normal to the cross-section). P is the Poynting vector, E k and E p are the kinetic and potential energies, defined by:

P α = -σ αβ uβ , E k = 1 2 ρ uα uα , E p = 1 2 ǫ αβ σ αβ . (22)
Once the eigensystem ( 13) is solved, the energy velocity of any eigenmode (k, U) can indeed be directly post-processed from SAFE matrices. First, it can easily be deduced that:

S Ēk dS = ω 2 4 Re(U T * MU), S Ēp dS = 1 4 Re{U T * (K 1 + ik(K 2 -K T 2 ) + k 2 K 3 )U}. ( 23
)
The computation of the cross-section and time averaged Poynting vector requires further developments but can also be simply expressed. Noticing that n = T, one has:

P • ndS = ω 2 Im{u * α σ αt dS}. (24)
Then, it can be checked that u * α σ αt = u T * L T s σ(1 + κx) = u T * L T s C(L xy + ikL s )u(1 + κx), so that the following useful expression holds for the averaged Poynting vector:

S P • ndS = ω 2 Im{U T * (K T 2 + ikK 3 )U}. ( 25
)
Expressions ( 23) and ( 25) allow a direct computation of the energy velocity defined by Eq. ( 21).

Preliminary results

The aim of this section is to verify that a twisting system yields the same physical results as those obtained with a straight coordinate system for a cylindrical waveguide, and similarly as those obtained with a helical system for a helical waveguide. The material is assumed to be isotropic, with no material damping. A value of 0.30 is chosen for the Poisson coefficient.

One considers a waveguide with a circular cross-section of radius a. The normalized frequency is given by ωa/c s , c s = E/2ρ(1 + ν) denoting the shear wave velocity. FE computations are held at fixed real wavenumbers k.

Six-node triangles meshes are used. The helix lay angle φ is defined from: tan φ = 2πR/L.

Cylindrical waveguide

For a cylindrical waveguide, we have φ = 0 and its axis corresponds to the Cartesian Z axis. As previously mentioned, the Cartesian system (X, Y, Z), which is the special case κ = τ = 0, obviously yields a translational invariance. As clearly shown in this subheading, a cylinder can also be analysed with a twisting coordinate system (κ = 0, τ = 0), where τ can be any value. Both cases (κ = 0, τ = 0) (reference solution) and (κ = 0, τ = 0) are computed. In order to check the adequacy of the twisting system, we choose a high torsion τ a = 0.5 (high rotation rate). For clarity, Fig. 2 (left) illustrates the cross-section FE mesh used (the same for both coordinate systems). However, one must be careful when interpreting the strong differences found for wavenumbers of flexural modes F (m,n), which are non-axisymmetric.

With a Cartesian system, the wavenumbers of F (m,n) modes occur in pairs of double roots. With a twisting system, their phase velocity become distinct due to the rotation of the (x, y) plane around Z (the axisymmetry is broken), their wavenumbers being translated by ±mτ a (see Fig. 4 for clarity). This was also obtained in Ref. [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF], where further explanations can be found.

Note that a wavenumber translation from ka to ka ± mτ a does not affect the slopes of curves in Fig. 3. Hence, the group velocity v g = ∂ω/∂k should remain unchanged with the twisting system. Because no damping is consid- Fig. 3). The wavenumber of the L(0,1) mode is left unchanged (axisymmetric mode).

Figure 5: Plot of normalized frequency vs. wavenumber for a cylindrical waveguide. Black:

Cartesian system (reference), gray: twisting system.

ered here, the group velocity is also equal to the energy velocity [START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF], and can be computed from Eqs. ( 21), ( 23) and ( 25). The energy velocity is plotted in Fig. 5: the same curves are obtained for the straight system and the twisting system. Consequently, both systems clearly yield equivalent results from a physical point of view.

Helical waveguide

One considers a waveguide with a helix radius R = 2a and a strong helix lay angle φ = π/4. This yields (κa = 0.25, τ a = 0.25) for the helical system and (κa = 0, τ a = 0.5) for the twisting system. The computed solution given by the helical system, for which the cross-section is circular (same mesh as Fig. 2 -left), is considered as the reference solution (see Refs. [START_REF] Treyssède | Numerical investigation of elastic modes of propagation in helical waveguides[END_REF][START_REF] Treyssède | Elastic waves in helical waveguides[END_REF]). Figure 2 (right) exhibits the cross-section mesh associated with the twisting system. This cross-section, corresponding to a Z =0 plane cut as parametrized by Eq. ( 10), is no longer circular as considered earlier.

One must also be careful when considering results obtained with the helical system: in order to transform results from a curvilinear helical s-axis to a straight Z -axis, the energy velocities must be divided by the step ratio l/L (the wavenumbers k must be multiplied by l/L). Figure 6 shows the dispersion curves computed for both the helical and twisting coordinate systems.

No difference is observed in wavenumber predictions, demonstrating the adequacy of the twisting system. Note that the differences previously observed for flexural modes inside a cylinder do not occur here because both systems rotate at the same rate along Z. Obviously, results for energy velocities are also identical.

Analysis of the seven-wire strand

Now wave modes propagating inside a typical seven-wire steel strand are studied both numerically and experimentally. The core wire radius is a=2.7mm. The helical wires have a radius equal to 0.967a and a pitch L=240mm, yielding a lay angle of φ = 7.9 21) is computed thanks to the useful formula [START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF] and ( 25).

Numerical results

The twisting coordinate system is (κa = 0, τ a = 0.0705). The FE mesh, corresponding to a Z =0 plane cut, is given by Fig. 7. Though not visible in the figure, there is no physical contact between peripheral helical wires and their cross-section is not exactly circular (even if this non-circularity has a negligible effect here). For a better accuracy of numerical results, the mesh has been strongly refined at physical contact points, yielding 6975 dofs. For simplicity, one assumes stick contact conditions between the core and outside wires (no slip, no separation and no friction are considered). This is equivalent to suppose that friction is high enough to prevent any slipping between wires, and to consider continuous displacement in every directions at the six contact points. The results are computed at fixed real wavenumbers k for the normalized frequency range [0;2]. We first compute the dispersion curves for single constitutive wires. Fig-

ure 8 shows the frequency vs. wavenumber and energy velocity vs. frequency plots for a single core wire and for a single peripheral wire. As already shown in previous papers [START_REF] Treyssède | Numerical investigation of elastic modes of propagation in helical waveguides[END_REF][START_REF] Treyssède | Elastic waves in helical waveguides[END_REF], minor differences are found between cylindrical and helical wires for small lay angles (7.9 • being a small lay angle). A very small velocity decrease can be observed because waves travel at slower velocity in the curved wire. It should be noted that the shift observed for the F (1, 2) mode is indeed mainly due to the fact that the cross-section radius of a peripheral wire is slightly different from that of the central one. An interesting point to be outlined for the helical wire is that a strong decrease of energy velocity occurs for the L(0, 1) mode at lowest frequencies. This decrease is confirmed by experiments in the next subsection. Note that this phenomenon also occurs for the torsional T (0, 1) mode (see Fig. 8).

Figure 9 gives the frequency vs. wavenumber plot and the energy velocity plot for the seven-wire strand. Compared to Fig. 8, a far more complex behaviour is observed resulting from strong interactions between wires. It clearly shows an apparent cut-off of the fastest mode (compressional-like L(0, 1) mode) around ωa/c s = 0.35, corresponding to 68kHz. This phenomenon is also confirmed by experiments in the next subsection. In fact, this apparent cut-off is due to a rapid veering of branches occurring near 0.35 -this can be more clearly observed in the normalized frequency vs.

wavenumber dispersion curves (see Fig. 9 -left). As a consequence, the upper curve of the energy velocity plot (Fig. 9 -right) is indeed composed of two distinct branches, describing the behaviour of two distinct wave modes.

Though somewhat subjective, a visual inspection of modeshapes has been performed in order to identify the evolution of both branches on a wider frequency range, as obtained in Fig. 10. The frequency vs. wavenumber plot exhibits a lower branch and an upper branch, denoted 1 and 2 respectively.

As observed in the energy velocity plot (Fig. 10 -right), the fastest mode corresponds to branch 1 at low frequencies and branch 2 at higher frequencies.

Comparing Figs. 8 and 10, it may be seen that the behaviour of branch 1 at low frequencies (1a portion) tends to the compressional L(0, 1) cylinder mode.

Its group velocity then quickly drops around the cut-off frequency of branch 2, near ωa/c s = 0.35. Branch 1 then asymptotes to the flexural F (1, 1) cylinder mode (1c portion), while branch 2 becomes the L(0, 1) cylinder mode (2c portion). It must be emphasized that both branches asymptotes to the energy velocity of the L(0, 1) mode inside a single wire (of radius a)

-see Figs. 8 and 10 -as opposed to a cylinder with an equivalent radius to the strand (of about 3a). Note that another veering happens to branch 2 around ωa/c s = 1.7, where the energy velocity suddenly drops.

Surprisingly, a similar trend can be observed in pipes [START_REF] Lowe | Defect detection in pipes using guided waves[END_REF], where the L(0,1) and L(0,2) modes play the role of branches 1 and 2: in pipes, the L(0,1) mode at low frequencies is the extensional type mode (corresponding to S0 in a plate); its group velocity drops at around the cut-off frequency of the L(0,2) mode; the L(0,1) mode then asymptotes to the flexural type mode velocity (A0 in a plate), while the L(0,2) mode becomes the extensional-like mode.

The phenomenon related to branch 1 and branch 2 could also be interpreted as curve veering [START_REF] Perkins | Comments on curve veering in eigenvalue problems[END_REF][START_REF] Pierre | Mode localization and eigenvalue loci veering in disordered structures[END_REF] (repulsion of branches, veering away from each other instead of crossing), encountered in eigenvalue problems of weakly coupled systems. Both branches swap properties around ωa/c s = 0.35 in a continuous but rapid way, and branch 2 seems to continue on the path previously followed by branch 1, which is characteristic of curve veering phenomena.

Figure 11 shows the modeshape evolution of branch 1 computed at points 1a, 1b and 1c. These points are located in Fig. 10 and correspond to the solutions (ka, ωa/c s , v e /c s ) =(0.1,0.16,1.56),(0.31,0.35,0.12) and (1.34,1.00,0.94) respectively. In this paper, modes are normalized with respect to the mass matrix M. Between points 1a and 1b (in the veering zone), it can be seen that the modeshape undergoes significant transformations. At point 1c, it becomes more complex and exhibits a combined flexural-torsional motion.

Figure 12 shows the modeshape evolution of branch 2, computed at points 2a, 2b and 2c (located in Fig. 10). These points correspond to the solutions (ka, ωa/c s , v e /c s )=(0.1,0.36,0.00), (0.2,0.37,0.35) and (0.63,1.00,1.55) respectively. Comparing the modeshapes of points 1a and 2c, one can conclude that branch 1 and branch 2 have interchanged their shapes (this is another feature of veering phenomena). Their global axial motion confirms that they are compressional-like modes. However, these modes are not exactly similar because the real parts of their axial displacement have opposite signs. Some computations have also been performed for a straight strand made of seven straight wires (results not shown for conciseness). Though dispersion curves are different from that of the helical strand, the veering phenomenon around 0.35 still exists in a straight strand. For a 7.9 • lay angle strand, we can conclude that the curvature of peripheral wires has a weak influence on this phenomenon compared to inter-wire contact effects.

Experiments

The final aim of this work is to perform a first validation of the seven-wire propagation model by comparing theoretical results to experiments in the low-frequency range for the compressional-like L(0, 1) modes. This was done on the basis of axisymmetric longitudinal guided waves measurements using magnetostrictive transducers. Explanations on this experimental device can be found in Laguerre et al. [START_REF] Laguerre | Low-frequency ultrasound reflectometry device based on magnetoelastic transducers for the non destructive evaluation of steel rods and cables[END_REF][START_REF] Laguerre | Magnetostrictive pulse-echo device for non-destructive evaluation of steel cylindrical materials using guide waves[END_REF].

The device consists of two small encircling dynamic coil transducers, one transmitter and one receiver, each being mounted coaxially with a solenoid coil for the static longitudinal polarizing magnetic field. The time pulse excitation to the transmitting transducer is a low radio-frequency burst with a frequency bandwidth equal to its central frequency. The detected time signal is then processed using a short-time Fourier transform (spectrogram method) to recover the group velocity dispersion curves. The experimental group velocity dispersion curve is then deduced from the spectrogram by dividing the known travelling distance separating two arrivals with the associated group-time delay at each frequency.

A first step is devoted to the analysis of the individual constituting wires, that are the cylindrical core wire and the peripheral helical wire respectively.

Figure 13 shows the comparison between measured and computed group velocities for the L(0, 1) mode of the single central cylindrical wire (radius of 2.7mm). For SAFE computations, note that group and energy velocities are equal because no damping was considered in this paper [START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF]. We observe a very good agreement between both curves, with the classical behaviour of the L(0, 1) mode group velocity, starting from the bar velocity estimated here at 5274m.s -1 and decreasing with frequency in this low-frequency range (less than 360kHz here).

The same procedure was applied to the measured time waveforms of the single peripheral helical wire (radius of 2.61mm, lay angle of 7.9 • ). In Fig. 14, the particular L(0, 1) low-frequency behaviour intrinsic to the helical geom- etry described in the previous subsection is very well reproduced (low frequencies propagate much slower than higher frequencies). Note that such a velocity decrease is also in accordance with experimental observations of Kwun et al. [START_REF] Kwun | Effects of tensile loading on the properties of elastic-wave propagation in a strand[END_REF].

Finally, the case of the seven-wire strand is considered. Numerical vs.

experimental comparisons were performed by superimposing the SAFE computed group time-delays at each frequency upon the measured spectrograms.

This only concerns the direct transmitter-to-receiver wave whose characteristics are not influenced by the strand-end reflections. A missing frequency band is observed in Fig. 15 in our experimental dataset. Such a missing band was initially experimentally found by Kwun et al. [START_REF] Kwun | Effects of tensile loading on the properties of elastic-wave propagation in a strand[END_REF], and referred to as 'notch frequency'. This phenomenon is well reproduced by the sevenwire strand model, assuming stick contact conditions between constitutive wires. As explained in the previous subsection, this missing frequency band is related to curve veering between two distinct wave modes. The veering As a final remark, it should be noted that the 'notch frequency' found by Kwun et al. [START_REF] Kwun | Effects of tensile loading on the properties of elastic-wave propagation in a strand[END_REF] in their experiments was around 80kHz at very low prestress, which is significantly greater than the one found in this paper. In fact, their strand had a nominal radius of 12.7mm with a pitch of 22cm. This gives a lay angle of 6.9 • : as explained in the previous subsection, the difference with our lay angle of 7.9 • can be neglected. However, their core wire radius was 2.16mm instead of 2.7mm here, and their bar velocity was slightly lower.

With a =2.16mm and E =2.1e+11Pa, the normalized frequency of 0.35 found in this paper yields a satisfying result of 83kHz. (A), (B) and (C) respectively denote the arrivals related to the electro-magnetic coupling between coils, the direct mechanical wave and the first reflection of the direct wave.

Conclusion

In this paper, a SAFE method has been proposed to analyse elastic guided wave propagation along multi-wire helical strands. It is is based on a specific twisting coordinate system that can be considered as a special case of the helical system. This system allows to preserve translational invariance along the straight axis of multi-wire waveguides made of a central straight wire and several peripheral helical wires. This has been checked by comparing reference results for a single cylindrical waveguide and for a single helical waveguide with those computed for a strongly twisting coordinate system. A useful expression has also been proposed for computing the energy velocity in a straightforward way. Dispersion inside a typical seven-wire strand has then been investigated by assuming stick contact conditions and results have been compared to experimental tests.

From a physical point of view, a far more complex behaviour than for single wires has been observed. It has been found that the fastest compressionallike mode inside the seven-wire strand looks like the L(0, 1) one for a single wire, excepted in a frequency region typically around 68kHz where a quick drop of energy velocity occurs due to curve veering between two distinct modes. This value agrees well with experimental results obtained in this paper with encircling longitudinally-polarized magnetostrictive transducers, where curve veering manifests itself as a missing frequency band whose central frequency is around 67kHz. Hence, the proposed seven-wire propagation model allows a relevant interpretation of this experimentally observed phenomenon, often referred to as 'notch frequency' in the literature. Further studies should deal with the effects of initial loading.

Figure 1 :

 1 Figure 1: Left: helix centreline curve (one step) with its associated Serret-Frenet basis and (x, y, s) helical coordinate system. Right: plane view. The point s = 0 lies in the (X, Y ) plane.

Figure 2 :

 2 Figure 2: Left: cross-section FE mesh of a cylindrical waveguide (and of a helical waveguide analysed with the helical coordinate system). Right: cross-section FE mesh of a helical waveguide R = 2a, φ = π/4 analysed with the twisting coordinate system.

Fig- ure 3

 3 exhibits the dispersion curves (normalized frequency vs. wavenumber) obtained from both systems. Wavenumbers of compressional L(0,n) and torsional T (0,n) modes are left unchanged because of their axisymmetry.

Figure 3 :

 3 Figure 3: Plot of normalized frequency vs. wavenumber for a cylindrical waveguide. Black:Cartesian system (reference), gray: twisting system.

Figure 4 :

 4 Figure 4: Wavenumber translation for some examples of F (m,n) modes (same legend as

Figure 6 :

 6 Figure 6: Plot of normalized frequency vs. wavenumber for the helical waveguide R = 2a, φ = π/4. Black: helical system (reference), gray: twisting system.

Figure 7 :

 7 Figure 7: Cross-section FE mesh for the seven-wire strand φ = 7.9 • (6975 dofs).

Figure 8 :

 8 Figure 8: Dispersion curves for single wires. Black: central wire, gray: peripheral wire. Left: plot of normalized frequency vs. wavenumber, right: plot of normalized energy velocity vs. frequency.

Figure 9 :

 9 Figure 9: Dispersion curves for the seven-wire strand. Left: plot of normalized frequency vs. wavenumber (inset: zoom on the veering region of compressional-like L(0, 1) modes), right: plot of normalized energy velocity vs. frequency.

Figure 10 :

 10 Figure 10: Branch identification of compressional-like L(0, 1) modes for the seven-wire strand (black: branch 1, gray: branch 2). Left: plot of normalized frequency vs. wavenumber. Right: plot of normalized energy velocity vs. frequency. The modeshapes corresponding to points 1a, 1b and 1c are given by Fig. 11, those corresponding to 2a, 2b and 2c are given by Fig. 12.

Figure 11 :

 11 Figure 11: From left to right: real part (top) and imaginary part (bottom) of the displacement modeshape of branch 1 computed at points 1a, 1b and 1c (located in Fig. 10). Gray meshes correspond to the undeformed cross-section. The view of imaginary part is 2D for clarity.

Figure 12 :

 12 Figure 12: From left to right: real part (top) and imaginary part (bottom) of the displacement modeshape of branch 2 computed at points 2a, 2b and 2c (located in Fig. 10). Gray meshes correspond to the undeformed cross-section. The view of imaginary part is 2D for clarity.

Figure 13 :

 13 Figure 13: Central wire (cylindrical). Normalized group velocity vs. frequency measured with the magnetostrictive device (x-mark) and computed with the SAFE model (solid line) for the L(0, 1) mode.

Figure 14 :

 14 Figure 14: Peripheral wire (helical). Normalized group velocity vs. frequency measured with the magnetostrictive device (x-mark) and computed with the SAFE model (solid line) for the L(0, 1) mode.

Figure 15 :

 15 Figure 15: Seven-wire strand. SAFE predictions (gray solid lines) of group-time delays for the compressional-like L(0, 1) modes superimposed to the experimental spectrogram.

  • . Mechanical properties are as follows: E=2.17e11Pa, ν=0.28, ρ=7800kg/m 3 . It is outlined that peripheral wires do not contact each other (this is a rather widespread design criterion for minimizing friction effects). As before, FE computations are performed at fixed real wavenumbers k and six-node triangles meshes are used. The energy velocity defined by Eq. (

Acknowledgements

The authors wish to thank the anonymous reviewers for their helpul comments.


Michèle Basseville
email: michele.basseville@irisa.fr

Frédéric Bourquin
email: bourquin@lcpc.fr

Laurent Mevel
email: laurent.mevel@inria.fr

; Houssein Nasser
email: houssein.nasser@tudor.lu

Fabien Treyssède
email: fabien.treyssede@lcpc.fr

Handling the Temperature Effect in Vibration Monitoring: Two Subspace-Based Analytical Approaches

Keywords: 3͑367͒ CE Database subject headings: Vibration, Damage, Assessment, Temperature effects, Analytical techniques, Structural dynamics

The dynamics of most civil engineering structures is affected by the ambient temperature. This raises the issue of discriminating changes in modal parameters due to damage from those due to such effects. A statistical parametric damage detection algorithm based on a null space residual associated with output-only subspace identification and a 2 test built on that residual has been designed by some of the writers. The purpose of this paper is to propose two extensions of this detection method which account for the temperature effect. The first extension uses a thermal model for deriving a temperature-adjusted null space. The second extension exploits the thermal model together with a statistical nuisance rejection technique. Both methods are illustrated on a laboratory test case within a climatic chamber.

Introduction

Detecting and localizing damages for monitoring the integrity of structural and mechanical systems are topics of growing interest due to the aging of many engineering constructions and structures and to increased safety norms. Automatic global vibration-based structural health monitoring ͑SHM͒ techniques have been recognized to be useful alternatives to visual inspections or local nondestructive evaluations performed manually ͑Pandey and Biswas 1994; [START_REF] Natke | Model-aided diagnosis of mechanical systems: Fundamentals, detection, localization, assessment[END_REF][START_REF] Farrar | Vibration-based structural damage identification[END_REF][START_REF] Catbas | Parameter estimation for multiple-input multiple-output modal analysis of large structures[END_REF][START_REF] Fritzen | Recent developments in vibration-based structural health monitoring[END_REF][START_REF] Alvandi | Assessment of vibration-based damage identification techniques[END_REF][START_REF] Gao | Experimental verification of the flexibility-based damage locating vector method[END_REF][START_REF] Zang | Structural health monitoring and damage assessment using frequency response correlation criteria[END_REF]. When processing vibration data ͑recorded with, e.g., accelerometers͒, a major focus is on extracting modal parameters ͑modal frequencies and associated damping values and mode shapes͒ which characterize the dynamics of the structure. It has been widely acknowledged that changes in frequencies bear useful information for damage detection, and information on changes in mode shapes is mandatory for performing damage localization.

It is also widely recognized that the dynamics of large civil structures, such as bridges, is significantly affected by the ambient temperature and other environmental effects, such as changes in loads, boundary conditions, moisture, wind, etc. Typically temperature variations have a major effect on the stiffness and thus on the dynamics of the structure, whereas the humidity mainly affects its mass. Moreover, changes in modal parameters due to environmental effects can be larger than changes due to damages ͑Moorty and Roeder 1992; [START_REF] Cornwell | Environmental variability of modal properties[END_REF][START_REF] Alampalli | Effects of testing, analysis, damage and environment on modal parameters[END_REF][START_REF] Rohrmann | Structural causes of temperature affected modal data of civil structures obtained by long time monitoring[END_REF]Peeters and De Roeck 2001;Kullaa 2003͒. This raises the issue of discriminating between changes in modal parameters due to damages and changes in modal parameters due to environmental effects. Handling environmental effects has recently become a major research issue for SHM. A number of methods have been designed for addressing those issues for different civil engineering application examples ͑Peeters et al. 2001; [START_REF] Golinval | A comparative study on damage detection of Z24-Bridge: One-year monitoring with varying environmental conditions[END_REF][START_REF] Kullaa | Structural health monitoring under variable environmental or operational conditions[END_REF][START_REF] Steenackers | Structural health monitoring of the Z24 Bridge in presence of environmental changes using modal analysis[END_REF]Yan et al. 2005a,b͒. Typical damage indicators used for SHM include modes, mode shape curvatures, transfer functions, autoregressive ͑AR͒ model coefficients, stiffness or flexibility matrix coefficients, etc. Those indicators generally result from the processing of measured data ͑identification͒ and depend on both the structural parameters and the environmental parameters. Most model-based approaches to the handling of the environmental parameters, and especially the temperature effect, involve a model of the effect of the temperature on the eigenfrequencies ͑modes͒ or other damage indicators. For example, an analytical model of the temperature-induced movements and associated stresses in a bridge has been proposed in Moorty and Roeder ͑1992͒, an AR model has been used for the temperature effect in Peeters and De Roeck ͑2001͒, and a linear adaptive filter has been proposed in [START_REF] Sohn | An experimental study of temperature effect on modal parameters of the Alamosa Canyon Bridge[END_REF].

The investigation performed for the large-scale Alamosa Canyon Bridge in New Mexico in Sohn et al. ͑1999͒ involves a linear adaptive filter. It has been found that the fundamental frequency has been modified by about 5% within a 1-day monitoring period, and this frequency is correlated with the temperature difference through the bridge deck. A simple linear filter, taking two time and two spatial temperature measurements as its inputs, has been used for modeling the variability of the frequencies over the day. A regression model has also been used for estimating the confidence interval of the frequencies for a new temperature profile. In Rohrmann et al. ͑2000͒, the environmental effects on the eigenfrequencies identified on the West End Bridge in Berlin, Germany, has been investigated for over a 3-year period. A thermomechanical model describing the temperature effect has been proposed. A frequency shift proportional to the mean temperature of the bridge and to the temperature gradient has been observed, and the higher order modes have been found to be more affected by the temperature effect than the lower order modes. This approach requires the measurement of the temperature over the whole bridge at any time, which is difficult to transpose when it comes to handling other environmental effects. In Peeters and De Roeck ͑2001͒, a regression analysis of the identified frequencies over the temperature has been performed for the Z24 Bridge in Switzerland. A series of damages has been realized. The relation between the frequencies and the temperature has been found to be linear for positive temperatures ͑in °C͒ and nonlinear for negative ones due to the variability of the asphalt elastic moduli, which significantly affects the bridge stiffness ͑Wood 1992͒. An AR model has been identified using the data associated with positive temperatures only. The prediction error, used as a damage indicator, succeeded in detecting the damage based on the first three frequencies.

In Kullaa ͑2003͒, factor analysis is used for deriving a relationship between the modes and the temperature. Basically factor analysis considers that the environmental effects are nonobservable and statistically independent variables, and performs a decomposition of the correlation matrix of monitored features. The damage indicator is a factor independent from the temperature ͑or the other environmental effects͒. This technique has been used for a crane vehicle, for which the five first eigenfrequencies compose the damage indicator. The method proposed in Vanlanduit et al. ͑2005͒ is of a similar kind. Two approaches are proposed in Deraemaeker et al. ͑2008͒. The first one consists in using modal filters for extracting damage indicators insensitive to the environmental parameters; the second consists in using factor analysis for eliminating the environmental effects from the damage indicator vector without measuring those effects; see also Kullaa ͑2004͒, Vanlanduit et al. ͑2005͒, and Yan et al. ͑2005a͒. The damage indicators are the eigenfrequencies, the mode shapes, and the Fourier transform of the output of modal filters used for reducing data measured with a large sensor array. The application example is a numerical model of a simulated bridge. The sensitivity of the damage indicators to noise levels, damages, and temperature gradients has been investigated for both approaches. The results from Kullaa ͑2003͒ and Deraemaeker et al. ͑2008͒ indicate that a factor analysis may help in performing damage detection without measuring the environmental parameters, provided that the dimension of the nonobservable vector is smaller than the dimension of the damage indicator vector.

A damage detection and localization method based on a subspace residual and 2 -type global and sensitivity tests has been proposed in Basseville et al. ͑2000, 2004͒ and successfully experimented on a variety of test cases ͑Mevel et al. 1999, 2003a,b;[START_REF] Peeters | Online vibration based crack detection during fatigue testing[END_REF][START_REF] Basseville | Output-only subspace-based damage detection-Application to a reticular structure[END_REF]. A subspace-based damage detection method inspired from that method has been proposed in [START_REF] Fritzen | Elimination of temperature effects on damage detection within a smart structures concept[END_REF] for handling the temperature effect. The main idea consists in saving different reference signatures identified under different temperature conditions, and in applying the subspace-based detection algorithm to new data with a reference modal parameter updated according to the actual temperature value. This technique has been applied to a composite plate within a furnace with controlled temperature. The test values have been shown to remain nearly unchanged for safe condition under different temperature levels and to be highly sensitive to structural damages. However, this approach is limited by the number of reference modal signatures available at different temperatures.

The purpose of this paper is to propose two different extensions of the subspace-based damage detection algorithm in Basseville et al. ͑2004͒, which address the above issues for handling the temperature effect, and to investigate their relevance when applied to a laboratory test case. The rationale for the two extensions is the following. The subspace-based residual involves the null space of a matrix ͑observability matrix͒ built on reference modes and mode shapes identified on the safe structure at a specific and known temperature. When processing data that will be recorded later at a different temperature that is also assumed to be known, it is mandatory to update the null space to account for possible modal changes due to temperature effects. This is achievable either by updating the stiffness matrix, thanks to a thermal model, and computing the associated modes and mode shapes, or by updating directly the modes and mode shapes using a semianalytic approach. This is the first extension. The second extension handles the temperature effect as a nuisance parameter. A model of the temperature effect developed by some of the writers is used together with a statistical nuisance rejection approach, for deriving an extended 2 test performing the rejection of the temperature effect seen as a nuisance for monitoring.

The paper is organized as follows. First, the key components of the subspace-based vibration monitoring algorithms are summarized. Second, an analytical model of the temperature effect on the modal properties is described for an Euler-Bernoulli beam, and analytical expressions for thermal sensitivities of the modal parameters are provided. Then, the handling of the temperature effect with the subspace-based approach is investigated and the two methods of temperature-adjusted subspace detection and temperature rejection, respectively, are proposed. Next, a laboratory test case is described, and numerical results obtained with these two new algorithms are reported. Finally, some conclusions are drawn.

Subspace-Based Detection

Some of the writers have advocated a subspace-based algorithm for modal monitoring. The key components of this algorithm are now recalled.

Subspace-Based Residual

It is well known ͑Juang 1994; Ewins 2000͒ that the continuoustime model of interest for vibration-based structural analysis and monitoring is of the form

MZ ¨͑t͒ + CZ ˙͑t͒ + KZ͑t͒ = v͑t͒ ͑ 1͒
where M , C , K = mass, damping, and stiffness matrices, respectively; Z͑t͒ = displacement vector at time t; and the unmeasured input force v͑t͒ϭstationary white noise. ͓In case of colored input excitation, the corresponding poles are generally highly damped and can easily be separated from structural modes. The nonstationary case is analyzed in Benveniste and Mevel ͑2007͒.͔ Moreover, the modes and mode shapes ⌿ are solutions of

det͑ 2 M + C + K͒ = 0, ͑ 2 M + C + K͒⌿ = 0 ͑2͒
It is also well known that, because of Eq. ͑1͒, vibration-based structural monitoring boils down to the problem of monitoring the eigenstructure of the state transition matrix F of a discrete-time linear dynamic system

ͭ X k+1 = FX k + V k+1 Y k = HX k ͮ ͑3͒ namely the pair ͑ , ͒ solution of det͑F -I͒ = 0, ͑F -I͒ = 0 ͑4͒
Actually, matrix F is related to the physical parameters M , C , K through

F = e L , L = ͫ 0 I -M -1 K -M -1 C ͬ where = sampling interval. The relation to the continuous-time eigenstructure ͑2͒ is = e and = def H contains the observed components of ⌿ . Let the ͑ , ͒'s be stacked into = def ͩ ⌳ vec⌽ ͪ ͑5͒
where ⌳ = vector whose elements are the 's; ⌽ = matrix whose columns are the 's; and vec= column stacking operator.

It is assumed that a reference parameter 0 is available, identified on data recorded on the reference ͑undamaged͒ system. Given a new data sample, the detection problem is to decide whether it is still well described by 0 or not. The subspace-based detection algorithm, a solution to that problem proposed in Basseville et al. ͑2000, 2004͒, builds on a residual associated with the output-only covariance-driven subspace-based identification algorithm.

This algorithm consists in computing the singular value decomposition of

H ˆp+1,q 0 = def R ˆ0 0 R ˆ1 0 ] R ˆq-1 0 R ˆ1 0 R ˆ2 0 ] R ˆq 0 ] ] ] ] R ˆp 0 R ˆp+1 0 ] R ˆp+q-1 0 ͑6͒
͑the empirical Hankel matrix filled with correlations computed on reference data corresponding to the undamaged structure state͒. Let 0 be the identified parameter. It can be characterized as follows ͑Basseville et al. 2000͒. Let O p+1 ͑͒ be the observability matrix written in the modal basis, namely

O p+1 ͑͒ = ΄ ⌽ ⌽⌬ ] ⌽⌬ p ΅ ͑7͒
where ⌬ = def diag͑⌳͒ and ⌳ and ⌽ are as in Eq. ͑5͒. Then, if the orthonormal matrix S is in the null space of the observability matrix associated with the reference condition, that is

S T O p+1 ͑ 0 ͒ = 0 ͑8͒
then it must also be in the null space of the Hankel matrix filled with correlations of the data taken from the reference model, namely

S T H ˆp+1,q 0 = 0 ͑9͒
Thus H ˆand O should have the same left null space. Note that matrix S in Eq. ͑8͒ depends implicitly on and is not unique; nevertheless it can be treated as a function S͑͒ ͑Basseville et al.

2000͒.

Let H ˆp+1,q be the empirical Hankel matrix filled with correlations of new data from the ͑possibly damaged͒ system. To check whether a damage has occurred, in other words if the new data are still well described by the reference parameter 0 , the idea is to check whether this new matrix H ˆp+1,q has the same left null space. For this purpose, a subspace-based residual is defined as

n ͑ 0 ͒ = def ͱ n vec͓S͑ 0 ͒ T H ˆp+1,q ͔ ͑ 10͒
where n = sample size.

Test for Global Modal Monitoring

Testing if = 0 holds true or not, or equivalently deciding that the residual n ͑ 0 ͒ is significantly different from zero, requires the probability distribution of n ͑ 0 ͒ which is generally unknown. A possible approach to circumvent this situation is the statistical local approach on which tests between the two hypotheses are assumed to be close to each other, namely

H 0 : = 0 and H 1 : = 0 + ␦/ ͱ n ͑11͒
where vector ␦ is unknown, but fixed. For large n, hypothesis H 1 corresponds to small deviations in . Define the mean sensitivity:

J͑ 0 ͒ = def ͉ -1/ ͱ nץ / ץE 0 n ͉͑͒ = 0 and the covariance matrix: ⌺͑ 0 ͒ = def lim n→ϱ E 0 ͑ n n T ͒.
Then, provided that ⌺͑ 0 ͒ is positive definite, and for all ␦, the residual n ͑ 0 ͒ in Eq. ͑10͒ is asymptotically Gaussian distributed under both hypotheses in Eq. ͑11͒ on ͑Basseville 1998͒, that is

n ͑ 0 ͒ --→ n→ϱ N͓J͑ 0 ͒␦,⌺͑ 0 ͔͒ ͑12͒
Thus a deviation ␦ in the system parameter is reflected into a change in the mean value of the residual n ͑ 0 ͒. Matrices J͑ 0 ͒ and ⌺͑ 0 ͒, which depend on neither n nor ␦, can be estimated prior to testing, using data on the reference system. Consistent estimates J ˆ͑ 0 ͒ and ⌺ ˆ͑ 0 ͒ of J͑ 0 ͒ and ⌺͑ 0 ͒ write as follows ͑Basseville et al. 2004; Zhang and Basseville 2003͒:

J ˆ͑ 0 ͒ = ͓I qr S͑ 0 ͒ T ͔͓H ˆp+1,q 0T O p+1 † ͑ 0 ͒ T I ͑p+1͒r ͔O p+1 Ј ͑ 0 ͒ ͑13͒ O p+1 Ј ͑ 0 ͒ = ⌳ 1 Ј ͑p͒ 1 0 ⌳ 1 ͑p͒ I r 0 0 ⌳ m Ј ͑p͒ m 0 ⌳ m ͑p͒ I r ⌳ i ͑p͒T = ͑1 i i 2 . . . i p ͒, 1 ഛ i ഛ m ⌳ i Ј ͑p͒T = ͑0 1 2 i . . . p i p-1 ͒, 1 ഛ i ഛ m ⌺ ˆ͑ 0 ͒ = 1/Kᐉ ͚ k=1 K k k T , n Ϸ Kᐉ with residual k computed using Y ͑k-1͒ᐉ+1 , . . . ,Y kᐉ ͑14͒ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 / 369
Assume that J͑ 0 ͒ is full column rank. ͓Typically, the rank of J͑͒ is equal to the dimension of , namely the number of modes plus the number of mode shapes times the number of sensors.͔ Then, deciding whether residual n ͑ 0 ͒ is significantly different from zero, stated as testing between the two hypotheses in Eq. ͑11͒, can be achieved with the aid of the generalized likelihood ratio ͑GLR͒ test. The GLR test computes the log-likelihood ratio of the distributions of the residual n ͑ 0 ͒ under each of the hypotheses in Eq. ͑11͒. Since both distributions are Gaussian, as displayed in Eq. ͑12͒, the GLR test in that case writes sup

␦ 0 ͕-͓ n ͑ 0 ͒ -J ˆ͑ 0 ͒␦͔ T ⌺ ˆ-1 ͑ 0 ͓͒ n ͑ 0 ͒ -J ˆ͑ 0 ͒␦͔͖ + n ͑ 0 ͒⌺ ˆ-1 ͑ 0 ͒ n ͑ 0 ͒
This boils down to the following 2 -test statistics:

n 2 ͑ 0 ͒ = def n T ͑ 0 ͒⌺ ˆ-1 ͑ 0 ͒J ˆ͑ 0 ͒ ϫ͓J T ͑ 0 ͒⌺ ˆ-1 ͑ 0 ͒J ˆ͑ 0 ͔͒ -1 J ˆT͑ 0 ͒⌺ ˆ-1 ͑ 0 ͒ n ͑ 0 ͒ ͑15͒
which should be compared to a threshold. In what follows, the âbove J and ⌺ and the index n when clear from context are dropped for simplicity.

The test n 2 ͑ 0 ͒ is asymptotically 2 distributed with rank͑J͒ degrees of freedom and noncentrality parameter ␦ T F͑ 0 ͒␦, where F͑ 0 ͒ is the asymptotic Fisher information on 0 contained

in n ͑ 0 ͒ F͑ 0 ͒ = def J͑ 0 ͒ T ⌺͑ 0 ͒ -1 J͑ 0 ͒ ͑ 16͒
2

Tests for Focused Monitoring

Focusing the monitoring on some subsets of components of the parameter vector is of interest for many purposes. For example, it may be desirable to monitor a specific mode and the associated mode shape . For this purpose, and without loss of generality, let , J, and F be partitioned as

= ͩ a b ͪ, J = ͑J a J b ͒, F = ͩ F aa F ab F ba F bb ͪ = ͩ J a T ⌺ -1 J a J a T ⌺ -1 J b J b T ⌺ -1 J a J b T ⌺ -1 J b ͪ ͑17͒
For deciding between ␦ a = 0 and ␦ a 0, two approaches may be used ͑Basseville 1997͒.

The sensitivity approach assumes b = 0 and consists in projecting the deviations in ␦ onto the subspace generated by the components ␦ a . This results in

˜a 2 = def ˜a T F aa -1 ˜a, where ˜a = def J a T ⌺ -1 ͑18͒
is the partial residual. Note that ˜a 2 might be sensitive to a change in b which is a limitation in the actual focusing capability of that test and that, for a b, residuals ˜a and ˜b are correlated ͑Basseville 1997͒.

The min-max approach consists in viewing parameter ␦ b as a nuisance, and rejecting it by replacing it in the likelihood with its least favorable value. This results in the test

a Ã2 = a ÃT F a Ã-1 a à ͑19͒
where the robust residual a à is defined by

a à = def ˜a -F ab F ab -1 ˜b ͑20͒
and where good candidates for performing the desired monitoring focused on a .

F a à = def F aa -F ab F bb -1 F ba .

Thermal Effect Modeling

We now introduce the key steps in modeling the temperature effect on the laboratory test case used in the application section.

Assuming an Euler-Bernoulli beam and neglecting rotating inertia, the eigenproblem governing the vibrations of a clamped planar beam, axially prestressed, is given by the following equations ͑Géradin and Rixen 1997͒:

Ά EI d 4 w͑x͒ dx 4 -N 0 d 2 w͑x͒ dx 2 -A 2 w͑x͒ = 0 ͉w͑x͉͒ x=0,L = 0; ͯ dw͑x͒ dx ͯ x=0,L = 0 · ͑21͒
Here ϭeigenpulsation and w͑x͒ϭcorresponding transversal displacement. The relation with the discrete-time mode and mode shape in Eq. ͑4͒ and continuous-time mode in Eq. ͑2͒ is

= def 1 arctan I͑͒ R͑͒ = I͑͒ w = w = def ͑22͒
In Eq. ͑21͒, E, I, , and A denote Young's modulus, the inertia momentum of the cross section, the density, and cross-sectional area, respectively ͑all assumed constant͒. N 0 ϭquasistatic axial preload given by

N 0 = EA͑x͒ ͑ 23͒ ͑x͒ = def 0 ͑x͒ -␣␦T, ␦T = def T 0 ͑x͒ -T ref ͑24͒
In Eq. ͑24͒, 0 ϭmechanical strain; T 0 ϭcurrent temperature; T ref ϭreference temperature for which there is no stress ͑corresponding to the temperature just before the tightening beam ends͒; and ␣ denotes the thermal expansion coefficient. It is important to note that N 0 remains spatially constant-even though a temperature gradient might exist along the beam-because no external axial body forces nor axial surface tractions are present during the experiments ͑in particular, gravity effects are negli-gible͒. Thus, only one measure at a given point is necessary. This is obtained by thermally compensated strain gauges, yielding a direct measure of ͑x͒ = N 0 / EA.

Updating the Modal Parameters

It can be shown that the eigenproblem in Eq. ͑21͒ yields an infinite number of solutions for ᐉ ͑ᐉ = 1 , . . . ,ϱ͒, given by the characteristic equation ͑Bokaian,

1988͒ 2␥ ᐉ -␥ ᐉ + ͓1 -cos͑␥ ᐉ -L͒cosh͑␥ ᐉ + L͔͒ + ͑␥ ᐉ +2 -␥ ᐉ -2 ͒sin͑␥ ᐉ -L͒sinh͑␥ ᐉ + L͒ = 0 ͑25͒
with the following notation:

␥ ᐉ Ϯ = ͫ ͱ A EI ᐉ 2 + ͩ N 0 2EI ͪ 2 Ϯ N 0 2EI ͬ 1/2

͑26͒

Eq. ͑25͒ must be solved numerically because no analytical solutions are available for the clamped case. The mode shapes are then given by

w ᐉ ͑x͒ = C ᐉͫ ͑cos ␥ ᐉ -x -cosh ␥ ᐉ + x͒ - cos ␥ ᐉ -L -cosh ␥ ᐉ + L ␥ ᐉ + sin ␥ ᐉ -L -␥ ᐉ -sinh ␥ ᐉ + L ͑␥ ᐉ + sin ␥ ᐉ -x -␥ ᐉ -sinh ␥ ᐉ + x͒ ͬ ͑27͒
where C ᐉ = normalization constant.

For each value of temperature shift ␦T, N 0 / EA, and thus ͑x͒ in Eq. ͑23͒, is measured; the eigenpulsations ᐉ 's are computed from Eqs. ͑25͒ and ͑26͒ and the mode shapes w ᐉ ͑x͒'s are computed from Eq. ͑27͒. The corresponding temperature-updated modal parameter T is then recovered from Eqs. ͑22͒ and ͑5͒. This is summarized as follows:

␦T → ͑x͒ → ͑ ᐉ ͒ ᐉ → ͑ ᐉ , ᐉ ͒ ᐉ=1,m → T ͑28͒
The number m of identified modal parameters ͑ ᐉ , ᐉ ͒ is limited by the number r of sensors and results from experimental to analytical mode matching.

In the study reported in this paper, the mode shapes w ᐉ ͑x͒'s are assumed to be unaffected by the temperature effect. Thus Formula ͑27͒ is not used, and the ᐉ 's in Eq. ͑28͒ are the mode shapes stacked in 0 .

Computing the Thermal Sensitivities

The variational formulation associated with the equations in Eq. ͑21͒ is written, for any complex number v, as

EI ͵ 0 L d 2 v dx 2 d 2 w dx 2 dx + N 0͵ 0 L dv dx dw dx dx -2 A ͵ 0 L vwdx = 0

͑29͒

with ͉v͉ x=0,L = ͉dv / dx͉ x=0,L = 0, where Lϭbeam length. A temperature change causes a variation ␦N 0 of the axial prestress, and thus a variation of the modes. The first variation of this formulation writes

EI ͵ 0 L d 2 v dx 2 d 2 ␦w dx 2 dx + ͵ 0 L ͩN 0 dv dx d␦w dx -2 Av␦w + ␦N 0 dv dx dw dx ͪdx = ␦͑ 2 ͒A ͵ 0 L vwdx ͑30͒
From Eq. ͑29͒, the sum of the first three terms in the left-hand side of Eq. ͑30͒ is equal to zero. Thus, writing Eq. ͑30͒ for v = w = w ᐉ , the sensitivity of the eigenfrequency ᐉ with respect to the axial prestress writes

␦͑ ᐉ 2 ͒ = ͐ 0 L ͑dw ᐉ /dx͒ 2 dx A͐ 0 L w ᐉ 2 dx ␦N 0 ͑31͒
Knowing the mode shapes w ᐉ , the expression in Eq. ͑31͒ can be computed semianalytically.

The sensitivity with respect to the temperature T is then obtained from ␦N 0 =-EA␣␦T, which results from Eqs. ͑23͒ and ͑24͒. From Eq. ͑31͒ and Basseville et al. ͑2004͒, it is deduced that the sensitivity of the ᐉth component ͑mode͒ of the parameter with respect to the temperature T writes

J T ͑ᐉ͒ = def ץ ᐉ ץ ᐉ ץ ᐉ ץ ᐉ ץ ᐉ ץT = ץ ᐉ ץ ᐉ ץ ᐉ ץ ᐉ 1 2 ᐉ ץ ᐉ 2 ץT = E␣͐ 0 L ͑dw ᐉ /dx͒ 2 dx 2 ᐉ ͐ 0 L w ᐉ 2 dx ͩ I͑͒ -R͑͒ ͪ ͑32͒
The sensitivity matrix ͑actually, a vector͒ of the modal parameter with respect to the temperature T has m components of the Form ͑32͒. The other components of J T are equal to zero since the prestress effects upon the mode shapes are negligible, at least for the temperatures considered during the experiment ͑see 3͒, and are thus neglected in what follows. Therefore J T writes

J T = ⎝ ⎜ ⎜ ⎜ ⎛ Ά΄ ] E␣ ͵ 0 L ͑dw ᐉ /dx͒ 2 dx 2 ᐉ ͵ 0 L w ᐉ 2 dx I͑͒ ] ΅ · m elements Ά΄ ] 0 ] ΅· r ϫ m elements Ά΄ ] - E␣ ͵ 0 L ͑dw ᐉ /dx͒ 2 dx 2 ᐉ ͵ 0 L w ᐉ 2 dx R͑͒ ΅ · m elements Ά΄ ] 0 ] ΅· r ϫ m elements ⎠ ⎟ ⎟ ⎟ ⎞

͑33͒

Two Analytical Approaches to Temperature Handling

In this section, two analytical approaches are proposed for handling the temperature effect in vibration-based SHM within the framework of the subspace-based detection algorithm ͓Eqs. ͑10͒ and ͑13͒-͑15͔͒. The first approach updates the null space S in Eq. ͑8͒ to account for possible modal changes due to temperature

JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 / 371
effects. The second approach handles the temperature as a nuisance parameter.

Temperature-Adjusted Null Space Detection

Having recorded both a new data sample ͑Y k ͒ k=1,. . .,n from the ͑possibly damaged͒ system and the corresponding temperature T, the idea is to update the modal parameter using Eq. ͑28͒, then to compute O p+1 ͑ T ͒ from Eq. ͑7͒ and the adjusted null space S͑ T ͒ such that

S T ͑ T ͒O p+1 ͑ T ͒ = 0 ͑34͒
Then, filling the empirical Hankel matrix ͓Eq. ͑6͔͒ with the correlation matrices of the new data, the residual is computed as

n ͑ T ͒ = def ͱ n vec͓S T ͑ T ͒H ˆp+1,q ͔ ͑ 35͒
and the proposed temperature-adjusted damage detection test n 2 ͑ T ͒ is computed as in Eq. ͑15͒.

Rejecting the Temperature Effect Seen as a Nuisance

The method used in this paper consists in using the material presented above for focused monitoring with the modal parameter and the temperature T playing the role of a and b , respectively. More precisely, assuming small deviations in both and T, the mean value of the residual n now writes

E 1 ͑ n ͒ = J͑ 0 ͒␦ + J T ␦T ͑36͒
where J͑ 0 ͒ is as in Eq. ͑12͒, J T is defined as

J T = def J͑ 0 ͉͒J T ͉ 0 ͑37͒
and where the matrix ͉J T ͉ 0 contains the sensitivities of the modes with respect to the temperature, computed at 0 from Eq. ͑33͒. Because of the expression for J T , computing the 2 test in Eq. ͑19͒ associated with the robust residual ͓Eq. ͑20͔͒ requires, in addition to J͑ 0 ͒, the computation of ͑ ͉J T ͉ 0 ͒.

Application

The global test ͓Eq. ͑15͔͒ and the two new algorithms of the previous section have been experimented on a vertical clamped beam equipped with accelerometers and installed within a climatic chamber, a laboratory test case provided by Laboratoire Central des Ponts et Chaussées.

Description of the Test Case

The experimental device is depicted in Fig. 1. A vertical beam is clamped at both ends on a workbench made of four vertical thick columns and two horizontal decks. This workbench is made of steel, whereas the beam is made of aluminum. The whole apparatus is set inside a climatic chamber with controlled ambient temperature. Because steel and aluminum have different thermal expansion coefficients ͑11.7ϫ 10 -6 and 23.4ϫ 10 -6 K -1 , respec-tively͒, a temperature change naturally induces a significant axial prestress, constant along the beam in the absence of external forces.

The test beam is instrumented with four accelerometers, located at nodes of the fifth flexural mode in order to avoid nodes of Modes 1-4 ͓a truncation up to the fourth mode has been performed in the analysis reported in Basseville et al. ͑2006͔͒. A pair of aluminum strain gauges with thermal compensation is also bonded on the beam. The half sum of both strains gives a direct measure of the actual axial prestress ͑divided by Young's modu-lus͒, without requiring any temperature measurement. This measurement is of primary importance for the method using a temperature-adjusted null space. Some temperature sensors have also been used in order to check that gauge measurements are closely related to thermal variations.

Tests are carried out inside the climatic chamber, first by stabilizing the ambient temperature for 1 h and then cooling down for 17 h with a slope of -1 °C / h. The beam is acoustically excited by a loudspeaker with a white noise input. Strains and temperature measurements are saved every second. Acceleration measurements are automatically triggered every 30 min for 600 s with a sampling frequency of 1,280 Hz, which is sufficient for Modes 1-4 ͑the fourth frequency is below 500 Hz͒.

An example of experimental results is displayed in Fig. 2. It clearly shows that prestress and thus frequencies increase as temperature decreases with time ͑note that prestress is taken positive when tensile͒. At the end of the test, the first frequency has been increased by about 16%. Modes 2, 3, and 4 ͑not shown here͒ have been increased by about 8, 5, and 3%, respectively. Axial prestress has thus a stronger effect for lower frequencies, which is consistent with known results ͑Bokaian 1988͒. As displayed in Fig. 3, the effect of the axial prestress on the first four mode shapes is hardly visible, and thus negligible, at least for the temperatures considered during the experiment. Hence the assumption made above on mode shapes unaffected by the prestress is valid in the present case. As for the damping coefficients, it is well known that they are subject to a larger uncertainty than the frequencies ͑Gersch 1974͒. It may happen that this uncertainty encompasses the temperature effects, but no knowledge is available about that.

Finally, two sets of experiments are performed, for the safe and damaged cases, respectively. In each set, 37 temperature scenarios are undertaken, and each scenario is repeated 10 times ͑one scenario every 30 min͒. The effect of a horizontal clamped spring attached to the beam plays the role of a local nondestructive damage. Indeed, the presence of the spring slightly modifies the frequencies, as a damage would do, causing an additional offset that is nearly constant throughout the experiment. This damage may be tuned by increasing the stiffness and/or the spring attachment point. In this experiment, the spring is located at L / 5 and its stiffness is about k = 4,000 N • m -1 . This choice has been designed from a finite-element model in order to give slight deviations ͑less than 1%͒ from the no-spring ͑safe͒ case. Note that the finite-element model used is simplified as a beam-spring system ͑the workbench is not accounted for͒. Experimental deviations of Modes 1-4 are about +0.8, +0.4, -0.1, and -0.2%, respectively. Though weak, the deviations for Modes 3 and 4 are slightly negative. This is probably due to the mechanical coupling of the beam with the spring bracket, which has a non-negligible flexibility compared to the workbench ͑this bracket is connected to the workbench and is not accounted for either in the finiteelement model͒. Fig. 4 shows how the first four frequencies are affected by the temperature scenarios, under both safe ͑blue͒ and damaged ͑red͒ conditions. In that theory, the prebending effects have been neglected upon dynamics, which is acceptable for experiments because cooling has been chosen instead of heating in order to avoid any prebuckling bending.

Implementation Issues

Some comments are in order on the selection of the key parameters and the computation of the key matrices for each of the two proposed methods. The reference parameter vector 0 , the null space S͑ 0 ͒ in Eq. ͑8͒, and the key matrices J ˆ͑ 0 ͒ in Eq. ͑13͒ and ⌺ ˆ͑ 0 ͒ in Eq. ͑14͒ for the global test ͓Eqs. ͑10͒ and ͑15͔͒ are estimated on the 10 realizations of the first scenario corresponding to a thermal constraint approximately equal to 20 m / m. The computations in Eqs. ͑28͒ and ͑33͒, for the modal updating and sensitivities, respectively, have been made with the following values of beam thickness and length: h = 0.01 m and L =1 m ͑depth: b = 0.03 m͒. Young's modulus and density are E = 6.95e +10 Pa and = 2,700 kg• m -3 . Those characteristics have been experimentally checked with static bending tests and also by comparing theoretical and experimental free eigenfrequencies.

Temperature-Adjusted Subspace Detection

Based on the temperature measurement for the current scenario, the adjusted parameter T in Eq. ͑28͒ can be computed using the analytical model ͓Eqs. ͑25͒-͑27͔͒. Since the mode shapes are assumed constant, this provides us with the evolution with respect to ͑x͒ of the frequencies only.

In the experiment, the clamped boundary condition-obtained with tightening jaws-is not perfect, which explains a negative offset in the experimental frequencies compared with theory ͑less than 3%͒. As displayed in Fig. 4, the computed and identified first frequencies are different, under both safe and damaged conditions, respectively. Comparing the black plots with respect to the corresponding blue ones, the frequencies in Eq. ͑28͒ appear slightly biased, and the accuracy of T can only be assumed in the vicinity of 0 . For this reason, based on an identified reference parameter ˜0, a modified parameter vector ˜T is computed, which contains mode shape values, damping coefficients, and unbiased frequencies. Actually, assuming that the slope of the black plots in Fig. 4 is correct, a simple shift helps getting rid of the bias

͑ ∀ ͒ ¯T͑͒ = def T ͑͒ + ¯0͑ 0 ͒ -T ͑ 0 ͒ ͑ 38͒
The proposed algorithm for the temperature-adjusted test based on the residual in Eq. ͑35͒ runs as follows. Having at hand a reference modal parameter 0 for the safe structure, and having recorded an n-size sample of new measurements ͑Y k ͒, namely for every realization of every temperature scenario: • Compute the adjusted parameter vector T in Eq. ͑28͒ and update it into ˜T from Eq. ͑38͒;

• Compute the observability matrix O͑ ˜T͒ from Eq. ͑7͒ and the null space S͑ ˜T͒ from Eq. ͑34͒; • Compute the Hankel matrix in Eq. ͑6͒ and the residual in Eq. ͑35͒, namely n ͑ ˜T͒ = def ͱ nvec͓S T ͑ ˜T͒H ˆ͔;

• Compute the Jacobian matrix J ˆ͑ ˜T͒ as in Eq. ͑13͒ and the covariance matrix ⌺ ˆ͑ ˜T͒ as in Eq. ͑14͒; and • Compute the 2 test:

n 2 ͑ ˜T͒ = def n T ͑ ˜T͒⌺ ˆ-1 ͑ ˜T͒J ˆ͑ ˜T͒ ϫ͓J ˆT͑ ˜T͒⌺ ˆ-1 ͑ ˜T͒J ˆ͑ ˜T͔͒ -1 J ˆT͑ ˜T͒⌺ ˆ-1 ͑ ˜T͒ n ͑ ˜T͒ ͑39͒

Temperature Rejection

For each value of temperature shift ␦T, N 0 / EA, and thus ͑x͒ in Eq. ͑23͒, is measured, and the ␦͑ ᐉ 2 ͒'s are computed from Eq. ͑31͒. The proposed algorithm runs as follows. Having at hand a long data sample for the safe structure and a reference modal parameter vector 0 identified on that data sample: • Compute the observability matrix O͑ 0 ͒ from Eq. ͑7͒ and the null space S͑ 0 ͒ from Eq. ͑8͒; • Compute the sensitivities with respect to 0 and T, namely J ˆ͑ 0 ͒ in Eq. ͑13͒, ͉J T ͉ 0 from Eq. ͑33͒, and J T in Eq. ͑37͒;

• Compute the covariance ⌺ ˆ͑ 0 ͒ in Eq. ͑14͒, and Fisher matrix in Eq. ͑16͒ for the entire parameter vector: and F = def

Temperature-Adjusted Subspace Detection

In Fig. 6, the values of the test ͓Eq. ͑39͔͒ are displayed. Whereas in Fig. 5 ͑left͒, it is not possible to discriminate between the damaged and the undamaged values for the test ͓Eq. ͑15͔͒; in Fig. 6 ͑left͒ it is possible to set a threshold for separating the test values in the damaged and safe scenarios. Fig. 6 ͑right͒ provides us with the same information about the dispersion of the test values, which are much higher in the damaged case, as theoretically predicted. For the plots in Fig. 6 ͑right͒, ˜0 has been estimated on scenarios corresponding to the middle value ͑x͒ = 71.497 m / m. It should be noticed that the choice of ˜0 for computing ˜T influences the behavior of the detection test.

Temperature Rejection

The values of the min-max 2 test ͓Eq. ͑40͔͒ are displayed in Fig. 7 for all the scenarios. Fig. 8 displays a zoom near the upper values of the reference scenario 0 corresponding to a prestrain of ͑x͒ -␣ = 20.344 m / m. A better visualization of the comparative capabilities and limitations of the min-max 2 test ͓Eq. ͑40͔͒ with respect to the global 2 test ͓Eq. ͑15͔͒ is achieved in Fig. 9. In Fig. 9 ͑left͒, it appears that the min-max test is smaller than the global one, as expected, and thus less sensitive to small increases in ͑x͒. In Fig. 9 ͑right͒, it is seen that damage detection can be achieved for larger values of when using the min-max test than when using the global one. Actually, the range for ͑x͒ where the tests values under safe conditions can be separated from the val-ues under damaged conditions extends from 40 m / m for the global test to 65 m / m for the min-max test. Fig. 9 ͑right͒ also shows that the nuisance rejection is only local, as predicted by the small deviation hypothesis, and fails when ͑x͒ too much deviates from the nominal value.

Comparison with Global Empirical Null Space Detection

The parametric algorithm proposed in this paper has been compared to a nonparametric approach to temperature handling which consists in using an empirical null space ͑Balmès et al. 2008͒. The following comments are in order. The two methods differ in two ways. In this paper, a parametric kernel is generated from a model, whereas for the nonparametric approach the kernel is learned from a collection of reference scenarios. The parametric kernel is expected to yield lower values for the safe scenario. The drawback is that one must know the value of ͑x͒ for each scenario, but there is no need for a training period. The second difference between the two methods is that the residual covariance is here computed for each scenario, and only once for the nonparametric method. The present approach is thus computationally more expensive. But it is also more robust with respect to changes in the excitation.

Conclusions

The problem of handling the temperature effect in vibration-based monitoring of civil engineering structures has been addressed. Two analytical approaches have been proposed, based on a subspace-based damage detection algorithm previously proposed by some of the writers and on a model of the temperature effect on the modal parameters. The first approach is based on a temperature-adjusted null space plugged into that damage detection algorithm. In the second approach, the temperature is seen as a nuisance effect, and the thermal model is used together with a statistical rejection approach to nuisance parameter handling. The relevance of the proposed algorithms has been illustrated on a laboratory test case within a climatic chamber.

The model of the temperature effect on the modal parameters considered in this article involves only an axial prestress, which is relevant for the application example. For more complex structures, thermal prestress in all directions should be taken into ac- count and the model should be extended to a three-dimensional temperature field. In such a case, a relevant parameterization for the thermal model ͓here Eqs. ͑23͒ and ͑24͔͒ should be selected, and the problem of handling the thermal parameters as nuisance parameters becomes more involved. Another issue for future research is the validation on structures in-operation. 2006 for the temperature-adjusted null space method, and at the 4th World Conference on Structural Control and Monitoring, San Diego, July 2006 for the temperature rejection method.

Notation

The following symbols are used in this paper:

A ϭ cross-sectional area; C ϭ damping matrix; E ϭ Young's modulus; E ϭ expectation operator; F ϭ Fisher information matrix; F, H ϭ 2m ϫ 2m state transition and r ϫ 2m observation matrices; f ϭ frequency ͑discrete-time model͒; H 0 ϭ null hypothesis for ; H 1 ϭ alternative hypothesis for ; H p,q ϭ Hankel matrix, with p block rows and q block columns; I ϭ inertia momentum of the beam cross section; J ϭ sensitivity matrix of with respect to ␦; K, M ϭ stiffness and mass matrices; L ϭ beam length; N 0 ϭ quasi-static axial preload ͑in newton͒; O p ϭ observability matrix, with p block rows; R i ϭ lag i correlation matrix of Y; S ϭ ͑p +1͒ • r orthogonal matrix; T ϭ temperature; V k , X k ϭ discrete-time state noise and state vector;

w ϭ transversal displacement associated with ; Y k ϭ discrete-time measurement vector; Z͑t͒ ϭ vector collecting the displacements of the degrees of freedom of the structure; ␣ ϭ thermal expansion coefficient; ⌬, ⌳ ϭ diagonal matrix and vector, respectively, filled with the 's; ␦ ϭ deviation in ; ͑x͒ ϭ strain; n , ˜a, a à ϭ residual, partial residual, and robust residual vectors; , 0 ϭ modal parameter vector and reference value; , ϭ eigenvalue of F and mode; ϭ beam density; ϭ sampling frequency; ⌺ ϭ covariance matrix of n ; ͑t͒ ϭ continuous-time state noise;

, ϭ eigenvector of F associated with and its observed components; n 2 , ˜a 2 , a Ã2 ϭ chi-square tests associated with n , ˜a, a * , respectively; ⌽ ϭ matrix whose columns are the 's; ⌿ ϭ mode shape associated with mode ; ϭ eigenpulsation ln ; and ˆϭ A ˆ: estimated value of a scalar, vector, or matrix A.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Workbench and instrumented test beam

Fig. 3 .

 3 Fig. 3. Temperature effect on the first four mode shapes

Fig. 4 .

 4 Fig. 4. Temperature effect on the four frequencies in both safe and damaged cases

F͑ 0 ,Fig. 6 .

 06 Fig. 6. 2 test ͓Eq. ͑39͔͒ with temperature-adjusted null space

Fig. 9 .

 9 Fig. 9. Comparing Fig. 5 ͑right͒ and Fig. 7 ͑right͒

  Note that a à is the residual from the regression of the partial residual ˜a in ␦ a on the partial residual ˜b in ␦ b ͑nuisance͒. Under both hypotheses in Eq. ͑11͒, test statistics a Ã2 is 2 distributed with l a degrees of freedom and noncentrality parameter ␦ a T F a à ␦ a under ␦ a 0 and for all ␦ b . This property makes the residual a à and the associated 2 -test a

	Ã2

/ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

/ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

/ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

/ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

/ JOURNAL OF ENGINEERING MECHANICS © ASCE / MARCH 2010 Downloaded 16 Feb 2010 to 137.121.21.30. Redistribution subject to ASCE license or copyright; see http://pubs.asce.org/copyright

Acknowledgments

This work was carried out within the Project No. 2003-54 CON-STRUCTIF, coordinated by IRISA, in the framework of the French Computer and Security program ͑ACI S&I͒. Preliminary versions of parts of this paper were presented at the 3rd European Workshop on Structural Health Monitoring, Granada, S., July