[figure]

Non destructive testing using non linear vibroacoustic

B. Rousselet, U.N.S.A. Laboratoire J.A. Dieudonné U.M.R. C.N.R.S. 6621, Parc Valrose, F 06108 Nice, Cédex 2, email : br@math.unice.fr G. Vanderborck Thales Underwater Systems, Département acoustique 06903 Sophia-Antipolis CEDEX and Laboratoire J.A. Dieudonné U.M.R. C.N.R.S. 6621

Copyright (C) february 2005 B. Rousselet, G. Vanderborck Verbatim copying and distribution is permitted in any medium, provided this notice is preserved.

1

[figure]

1 Introduction

Several recent experimental studies show that it is possible to detect defects in a structure by considering its vibro- acoustic response to an external actuation.

1.1 Some previous papers

On this topic there is a vast literature in applied physics. We recall some papers related to the use of the frequency response for non destructive testing; in particular generation of higher harmonics, cross-modulation of a high frequency by a low frequency:

[figure]

• In Ekimov-Didenkulov-Kasakov (1999), [2], the authors report experiments with torsional waves in a rod with a crack: they use HF torsional wave (20kHz) and a LF flexural wave (12 Hz).

• In Zaitsev-Sas (1999), [9], the authors report experiments with plate vibration submitted to LF (20-60Hz) vibration by a shaker and HF (15-30 kHz) oscillations by a piezo-actuator. They notice that weak modulation side-lobes are observed for the undamaged sample but drastic increase in nonlinear vibro-acoustic of the damaged sample. Some theoretical explanations are provided.

• Other results may be found in Sedunov-Tsionsky-Donskoy(2002) [3],Sutin-Donskoy (1998), [1], Moussatov-Castagnede-Gusev(2002), [5] ...

• GDR 2501 (Etude de la propagation ultrasonore en milieux inhomognes en vue du controle non destructif)

3

[figure]

In Vanderborck-Lagier-Groby (2003) [8], ”a vibro-acoustic method, based on frequency modulation, is developed in order to detect defects on aluminum and concrete beams. Flexural waves are generated at two very separated frequencies by the way of two piezoelectric transducers. The low one corresponds to the first resonance f_{m}, the second one to a high non modal frequency f_{p}. The nonlinear response, due to the defects inside the structure, is detected by non-zero flexural waves at f_{p} ± nf_{m} frequencies.

see Vanderborck-Lagier(2004) ⟼ beam experimentation

[figure]

Very recent experiments

• have been performed on a real bridge by G. Vanderborck with four prestressed cables: two undamaged cables, a damaged cable and a safe one but damaged at the anchor;

• these experiments have been performed in the frame of the European program “Promoting competitive and sustainable growth ”of 15/12/99.

• The cables are roughly 100 m long, 4 tones weight, 15cm in diameters.

• The experiments have proved the presence of the damaged cable but also the safe one damaged at the anchor.

• Routine experimental checking with the lower eigenfrequencies had only proved only the presence of the very damaged cable by comparison with data collected 15 years ago.

• See Vanderborck-Lagier(2004) [10] for a presentation of the results of the experiment with a new post processing graphic presentation of experimental results.

5

[figure]

Difficulties of the experiments:

• non linearities of the shakers (including piezoelectric actuators)

• Natural non linearities: supports, links of complex multi structures as air planes, bridges etc

[figure]

Orientation

We intend to present simple spring mass models, simple bar and beam models with damage and use asymptotic expansions and numerical methods to try to get results which show some similarity with the experiments of [8]. Asymptotic expansions have been used for at least a century and for example has been used recently for numerical approximation of bifurcation of structures in PotierFerry-Cochelin and coworkers (1993) [4].

The key idea is to look at the solution in the frequency domain for the experiments and consequently for the numerics.

In a paper to be submitted (Lagier-Vandeborck) [7] several types of nonlinearities of defects are considered: contact elasticity, threshold contact model, nonlinear filling material. This last case will be considered for bar models: it may happen in case of corrosion: the voided crack is filled by a new dusty material: then the elastic crack response is related to the elastic properties of the filler. In this case it seems reasonable to consider a nonlinear elastic relation for the filler.

7

[figure]

2 Background of Fourier transform

2.1 Basic formulas

2.1.1 Fourier serie

For a detailed presentation, see for example Gasquet-Witomski [11] and for an engineering view point Lathy [6]; for a function f of period T, its expansion in fourier serie is:

f(t) = ∑_{n=−∞}^{+∞} c_{n}e^{\frac{2πint}{T}} with c_{n} = \frac{1}{T} ∫_{0}^{T} f(t)e^{\frac{−2πint}{T}} dt  (2.1)

Discrete Fourier transform: D.F.T. which may be computed quickly by the algorihm of F.F.T. To a sequence (y_{k}), k = 0, …, N − 1, is associated an other sequence (Y_{n}), n = 0, …, N − 1 with the formulas:

Y_{n} = \frac{T}{NT} ∑_{k=0}^{N−1} y_{k}e^{\frac{−2πin kT}{T N}} y_{k} = ∑_{n=0}^{N−1} Y_{n}e^{\frac{2πink}{N}}  (2.2)

[figure]

An approximation of Fourier coefficients (caution to indixes) may be obtained with :

c_{n}^{N} = \frac{T}{NT} ∑_{k=0}^{N−1} f(\frac{kT}{N} )e^{\frac{2πin kT}{T N}} or with the DFT using  (2.3)

y_{k} = f(\frac{kT}{N} ) c_{n} ≃ c_{n}^{N} = \{ \begin{matrix}Y_{n} for 0 ≤ n < \frac{N}{2} \\ Y_{n+N} for − \frac{N}{2} ≤ n < 0 \end{matrix}  (2.4)

with the error approximation:

c_{n}^{N} = ∑_{q≠0} c_{n+qN}  (2.5)

Beware to Scilab FFT which provides: X_{n} = ∑_{k=0}^{N−1} y_{k}e^{\frac{2πin k}{N}} ie the Fourier coefficients are computed from: Y_{n} = \frac{1}{N} X_{n}.

9

[figure]

2.1.2 Fourier transform

(Ff)(ν) = f̂(ν) = ∫_{R} f(t)e^{−2πiνt}dt  (2.6)
(F^{−1}g)(t) = (F̄g)(t) = ∫_{R} g(ν)e^{+2πiνt}dν  (2.7)

f̂^{m} = (2πiν)^{m}f̂ F(−(2πit)^{(m)}f(t)) = f̂^{m}(ν)  (2.8)

(2.9)

[figure]

{ \atop χ_{[−A,A]}} = \frac{sin(2πνA)}{πν} = 2A.sinc(2πνA) { \atop χ_{[0,A]}} = e^{−iπAν} \frac{sin(πAν)}{πν} = Asinc(πAν)

(2.10)

with the sampling function (“sinus cardinal”) sinc(t) = \frac{sin(t)}{t}  (2.11)
F(e^{2πiat}) = δ_{a}, and F(e^{−2πiat}T) = τ_{a}T̂ = δ_{a} ∗ T̂  (2.12)
(F(cos(2πiat)T) = \frac{1}{2}(τ_{a}T̂ + τ_{−a}T̂) = \frac{1}{2}(δ_{a} ∗ T̂ + δ_{−a} ∗ T̂)  (2.13)
F(cos(2πiat)χ_{[−A,A]}) = A(τ_{a}.sinc(2πνA) + τ_{−a}.sinc(2πνA))  (2.14)
F(sin(2πiat)χ_{[−A,A]}) = iA(−τ_{a}sinc(2πνA) + τ_{−a}sinc(2πνA))  (2.15)
F(cos(2πiat)χ_{[0,A]}) = \frac{A}{2}(τ_{a}e^{−iπAν}sinc(πνA) + τ_{−a}e^{−iπAν}sinc(2πνA))  (2.16)

11

[figure: 12]

[figure: 13]

[figure]

2.2 Sampling

Δ_{a} = ∑_{n∈Z} δ_{na} Δ̂_{a} = ∑_{n∈Z} e^{−2πi n aν} = \frac{1}{a}Δ_{\frac{1}{a}} le “peigne”  (2.17)
also Δ_{a} = \frac{1}{a} ∑_{n∈Z} e^{2πin \frac{t}{a}}  (2.18)

a sampling of f is: afΔ_{a} = a ∑_{Z} f(na)δ_{na} with a the sampling period. Poisson formula:

∑_{Z} f(t − na) = \frac{1}{a} ∑_{Z} f̂(\frac{n}{a})e^{2πin \frac{t}{a}} (f distribution à support compact)  (2.19)
∑_{Z} ĝ(ν − \frac{n}{a}) = a ∑_{Z} g(na)e^{−2πinνa} (ĝ distribution à support compact)  (2.20)
a(f̂Δ_{a})(ν) = ∑_{n∈Z} f̂(ν − \frac{n}{a}) = a ∑_{n∈Z} f(na)e^{−2πiνna}  (2.21)

( f tempered and f̂ and f̂ with support in [−ν_{c}, ν_{c}] )

no overlap of the spectrum for a ≤ \frac{1}{2ν_{c}}.

[figure]

Shanon For

f̂ ∈ L^{2}(R) et Supp(f̂) ⊂ [−ν_{c}, ν_{c}]  (2.22)
∀a ≤ \frac{1}{2ν_{c}} f(t) = ∑_{Z} f(na)sinc(\frac{π}{a} (t − na))  (2.23)
with the sampling function sinc(t) = \frac{sin(t)}{t}  (2.24)

Use a low-pass filter before sampling.

15

[figure]

2.3 Numerical computation of Fourier transform

f̂(ν) ≃ ∫_{−T/2}^{T/2} f(t)e^{−2πiνt}dt = T c_{νT} or with ν = \frac{n}{T}  (2.25)
f̂(\frac{n}{T}) ≃ T c_{n} for \frac{−N}{2} ≤ n ≤ \frac{N}{2}  (2.26)

Fourier coefficients c_{n} are numerically computed with FFT where y_{k} = f(\frac{kT}{N}):

f̂(ν) ≃ ∑_{k} f(\frac{kT}{N} )exp(−2πiν \frac{kT}{N} ) with ν = \frac{n}{T}  (2.27)
≃ ∑_{k} f(\frac{kT}{N} )exp(−2πin\frac{k}{N}) with n = νT  (2.28)

with sampling period \frac{T}{N} ≤ \frac{1}{2ν_{c}} (no overlap of the spectrum) but ... a function of compact support in time is not of compact support in frequency...

2.4 Exemples

Gate function Commenons par un exemple classique: f = χ_{[0,b]}, sa transformée de Fourier est:

f̂(ν) = \frac{sin(πbν)}{πν}exp(−iπbν)  (2.29)

Avec b = 1/2 et 1000 points utilisés dans l’intervalle [0, 1], on trouve les transformées sur les figures ci jointe. On pourra remarquer que le maximum est correct

Cosinus Pour la fonction cos(2πt), il est bien connu que la transformée de Fourier est δ_{1} +δ_{−1}. la transformée discrète est elle mme une approximation numérique de ∫_{0}^{T} cos(2πt)exp(−2πiνt)dt et l’on trouve un pic de hauteur la moitié de l’intervalle d’intégration. On trouve dans les figures 3 et 4, les transformée de fourier discrète calculée dans [0, 1] puis [0, 20]

Exponentielle-cosinus On constate que la transformée de Fourier de exp(10^{−3}t)∗cos(2πt) est sensiblement égale à celle du cosinus tandis que

16-1

Figure 1: Norm of the Fourier transform of the gate χ_{[0,0.5]} in [0, 9] herz 16-2 Figure 2: Norm of the Fourier transform of the gate χ_{[0,0.5]} in [0, 1000] herz

16-3

Figure 3: Norm of the Fourier transform of cos(2πt) in [0, 1000] herz 16-4 Figure 4: Norm of the Fourier transform of cos(2πt) in [0, 1000] herz

16-5

celle de exp(10^{−2}t) ∗ cos(2πt) est un peu différente, voir les figures 5 et 6 somme de deux sinus

sin_p_sin.sci

16-6 Figure 5: Norm of the Fourier transform of exp(10^{−3}t) ∗ cos(2πt) in

[0, 1000] herz 16-7

Figure 6: Norm of the Fourier transform of exp(10^{−2}t) ∗ cos(2πt) in

[0, 1000] herz 16-8

Figure 7: Norm of the Fourier transform of sin(2πt) + .2sin(20πt) in

[0, 50] herz 16-9

[figure]

3 Simplest mechanical example

in which we can exhibit intermodulations.

We consider a 1 d.o.f example of a spring mass system with a non linear spring.

mÿ = −ky − k_{3}y^{3} + ϵF sin(αt) or  (3.1)
ÿ = −ω_{0}^{2}y − \frac{k_{3}}{m} y^{3} + ϵ\frac{F}{m}sin(αt) with ω_{0}^{2} = \frac{k}{m}  (3.2)
with initial conditions y(0) = ϵη_{1}, ẏ = ϵv_{1}  (3.3)

We are going to solve this equation symbolically with an asymptotic expansion with respect to ϵ: y = ϵy_{1} + ϵ^{2}y_{2} + ϵ^{3}y_{3} + …; then numerically...

[figure]

3.1 The linear case

The first term is solution of:

ÿ_{1} = −ω_{0}^{2}y_{1} + \frac{F}{m}sin(αt) with y(0) = η_{1}, ẏ = v_{1} which gives:  (3.4)
y_{1} = Ae^{iω_{0}t} + \overline{A}e^{−iω_{0}t} + De^{iαt} + \overline{D}e^{−iαt} with  (3.5)
A = \frac{η_{1}}{2} − \frac{iF α}{2mω_{0}(α^{2} − ω_{0}^{2})} − \frac{iv_{1}}{2ω_{0}} D = −\overline{D} = −\frac{F}{2im(α^{2} − ω_{0}^{2})} or  (3.6)
with ϕ = \frac{F}{m(α^{2} − ω_{0}^{2})}, y_{1} = η_{1}cos(ω_{0}t) + (\frac{v_{1}}{ω_{0}} + \frac{α}{ω_{0}} ϕ)sin(ω_{0}t) − ϕsin(αt) (3.7)

Remarque 1. If we set η_{1} = 0, v_{1} = 0, then the term of pulsation ω_{0} has magnitude \frac{α}{ω_{0}} times the magnitude of the term of pulsation α; this is not a good choice for the non linear case in which \frac{α}{ω_{0}} is of order 100; it seems good choice η_{1} = 0, v_{1} = ω_{0}(−\frac{α}{ω_{0}} + 1)\frac{F}{m(α^{2}−ω_{0}^{2})} = (−α + ω_{0})ϕ

18

[figure]

3.2 Other terms

The term y_{2} is zero but the third term satisfies:

ÿ_{3} = −ω_{0}^{2}y_{3} + \frac{k_{1}}{m} y_{1}^{3}  (3.8)

to simplify, we assume η_{1} = 0 and set ϕ = \frac{F}{m(α^{2}−ω_{0}^{2})}, ψ = \frac{v_{1}}{ω_{0}} + \frac{αϕ}{ω_{0}} so that:

y_{1}^{3} = (ψsin(ω_{0}t) − ϕsin(αt))^{3} =  (3.9)
−1/4ψ^{3}sin(3ω_{0}t) + 3/2ψ(1/4ψ^{2} + ϕ^{2})sin(ω_{0}t)  (3.10)
−3/4ψϕ^{2}sin((ω_{0} + 2α)t) − 3/4ψϕ^{2}sin((ω_{0} − 2α)t)  (3.11)
+3/4ψ^{2}ϕsin((α + 2ω_{0})t) + 3/4ψ^{2}ϕsin((α − 2ω_{0})t)+  (3.12)
−3/2(ψ^{2}ϕ + 1/2ϕ^{3})sin(αt) + 1/4ϕ^{3}sin(3αt)  (3.13)

(3.14)

If we go on in the expansion, we get terms of pulsation α + 4ω_{0}, α + 6ω_{0} etc

[figure]

3.2.1 Numerical issues

For α = 40π, ω_{0} = 2π, F = 100, ϕ = .6348445087e − 2, αϕ = .7977691380, (3.15)
v_{1} = −.7578806812, ψ = −5.366518580, 3/4ψ^{2}ϕ = −.1371241364  (3.16)
ϕ^{3} ≪ ψ^{3}  (3.17)

General tendency:

• The pick of 3ω_{0} is much larger than the pick in α ± 2ω_{0} which are the most natural picks in the experiments;

• it is delicate to find datas such that the secondary picks at α ± 2ω_{0} actually appear when the differential equation is solved numerically.

• Question: algorithm and software for detecting the secondary picks?

• then find (by optimization) datas such that the secondary picks are important: criteria for damage.

20

[figure]

Figure 8: Linear response y, ν_{ω_{0}} = 1

[figure]

Figure 9: Fourier of linear response, ν_{ω_{0}} = 1

22

[figure]

fig y1.ps etc      :
k=(2*%pi)^2; m=1; nua=10; alfa=2*%pi*nua; F=1450; lam=.2;
A=[0 1;-k/m 0]; dt=.01; tmax=6*%pi;
k1=950; sk1=" ,k1="+string(k1);
function [xdot]=f3(t,x)
xdot=A*x+ [0;1 ]*( ( -k1*([1,0]*x)^3  )+
 F*sin(alfa*t)  - lam*[0,1]*x)
endfunction
function [Jf3]=jacf3(t,y)
Jf3=A+  [0,0; ( -3*k1*([1,0]*y)^2  ) , - lam ]
endfunction

[figure]

Figure 10: zoom of non linear response y, ν_{ω_{0}} = 1

24

[figure]

Figure 11: non linear response y, ν_{ω_{0}} = 1

[figure]

Figure 12: Fourier of non linear response y, ν_{ω_{0}} = 1

26

[figure]

Figure 13: Fourier of non linear response y ν_{ω_{0}} = 1

[figure]

4 Transverse vibrations: vibrating masses on streched cables in large displacement

Work performed with Theissen (doctoral student of U. Muenster); Erasmus students N. Goris and I. Altrogge worked on this topic during their stay in UNSA (2004-2005). We consider n masses attached to horizontal springs (or cables) which are in tension T_{0}, at rest ; the tension is positive when the cable is in traction which is assumed; at rest the mass m_{i} is submited to the force T the masses are moving (vertically) transversely to the springs; we denote by uper case letters quantities in the rest position and lower case in the current configuration.

4.1 Masses in vertical displacement

Here we assume that the masses can move only vertically.

• L_{i} lenth at rest; l_{i} lenth at time t; as the masses are moving vertically: l_{i}^{2} = L_{i}^{2} + (y_{i} − y_{i−1})^{2}

• and the change of tension of the linear elastic spring due to the change of of lenth T_{i} = T_{0} + k_{i}[l_{i}(y) − L_{i}] = T_{0} + k_{i}(\sqrt{L_{i}^{2} + (y_{i} − y_{i−1})^{2}} − L_{i}). this tension is directed

[figure]

along the axis of the spring.

• Denote by θ_{i}, the angle of the spring with the horizontal axis, we have

• y_{1} = L_{1}tan(θ_{1}), y_{i} − y_{i−1} = L_{i}tan(θ_{i}) y_{n} − y_{n−1} = L_{n}tan(θ_{n}). We enforce y_{n} = 0. See the picture with two masses and 3 cables.

30

[figure]

The equation of the dynamics:

m_{i}y_{i}” = −T_{i}sin(θ_{i}) + T_{i+1}sin(θ_{i+1}) + u_{i} i = 1 … n  (4.1)

where −T_{i}sin(θ_{i}) + T_{i+1}sin(θ_{i+1}) is the vertical component of the force acting on mass i; we assume that there is no horizontal movement so the horizontal component of the force does not work. The applied load on mass i is denoted by u_{i}; it is the control to be determined.

[figure]

Set

ζ_{i} = \frac{(y_{i} − y_{i−1})}{L_{i}}, and note that sin(arctan(ζ_{i})) = \frac{ζ_{i}}{\sqrt{1 + ζ_{i}^{2}}} so that  (4.2)
T_{i}sin(θ_{i}) = (T_{0} + k_{i}(L_{i}\sqrt{1 + ζ_{i}^{2}} − L_{i})) \frac{ζ_{i}}{\sqrt{1 + ζ_{i}^{2}}} =  (4.3)
(T_{0} − k_{i}L_{i})\frac{ζ_{i}}{\sqrt{1 + ζ_{i}^{2}}} + k_{i}L_{i}ζ_{i}  (4.4)

Possible approximations:

T_{i}sin(θ_{i}) = (T_{0} − k_{i}L_{i})(ζ_{i} − \frac{1}{2}ζ_{i}^{3} + \frac{3}{8}ζ_{i}^{5} + O (ζ_{i}^{7})) + k_{i}L_{i}ζ_{i} =  (4.5)
T_{0}ζ_{i} + (T_{0} − k_{i}L_{i})(−\frac{1}{2}ζ_{i}^{3} + \frac{3}{8}ζ_{i}^{5} + O (ζ_{i}^{7}))  (4.6)

Same expansion for T_{i+1}sin(θ_{i+1}) with ζ_{i+1} = \frac{(y_{i+1}−y_{i})}{L_{i+1}}

32

[figure]

4.1.1 Linearized equation

m_{i}y′′_{i} = −T_{0} (\frac{(y_{i} − y_{i−1})}{L_{i}} + \frac{(y_{i+1} − y_{i})}{L_{i+1}}) + u_{i}

corrector equations may be obtained; details for 1 d.o.f below.

\begin{matrix}l_{i}(y) − L_{i} = \frac{(y_{i} − y_{i−1})^{2}}{2L_{i}} − \frac{(y_{i} − y_{i−1})^{4}}{8L_{i}^{3}} + O((y_{i} − y_{i−1})^{6}) & (4.7) \\ sin(θ_{i}) = sin(atan(\frac{y_{i} − y_{i−1}}{L_{i}}) = & (4.8) \\ \frac{y_{i} − y_{i−1}}{L_{i}} − \frac{(y_{i} − y_{i−1})^{3}}{2L_{i}} + \frac{3(y_{i} − y_{i−1})^{5}}{8L_{i}} + O((\frac{y_{i} − y_{i−1}}{L_{i}}))^{7} & (4.9) \end{matrix} ; 33-1

[figure]

4.2 Case with 1 d.o.f

4.2.1 Model with 1 d.o.f In this case, with y_{0} = 0, y_{2} = 0 we have

m_{1}y_{1}” = −T_{1}sin(θ_{1}) + T_{2}sin(θ_{2}) + u_{1}  (4.10)

with θ_{1} = atan(y_{1}/L_{1}), θ_{2} = −atan(y_{1}/L_{2})

m_{1}y_{1}” = −T_{1}sin(atan(\frac{y_{1}}{L_{1}} )) − T_{2}sin(atan(\frac{y_{1}}{L_{2}} )) + u_{1}  (4.11)

Linearized equation

m_{1}y_{1}” = −T_{0}(\frac{1}{L_{1}} + \frac{1}{L_{2}} )y_{1} + u_{1}  (4.12)

The numerical solution of this model may be performed without stiff hypothesis with scilab routine ode; (sin(tan) is Lipshitz) but it is not obvious to prescribe the right mechanical constants to obtain clear intermodulation peaks;

also trouble of the experiments!

[figure]

4.2.2 Approximation Here set ζ_{1} = \frac{y_{1}}{L_{1}} , ζ_{2} = −\frac{y_{1}}{L_{2}} . Start from previous approximation

−T_{1}sin(θ_{1}) + T_{2}sin(θ_{2}) =  (4.13)

T_{0}(ζ_{2} − ζ_{1}) − (T_{0} − k_{1}L_{1})(−\frac{ζ_{1}^{3}}{2} + \frac{3ζ_{1}^{5}}{8} ) + (T_{0} − k_{2}L_{2})(−\frac{ζ_{2}^{3}}{2} + \frac{3ζ_{2}^{5}}{8} ) + O(ζ_{1}^{7} + ζ_{2}^{7}),(4.14)

expand y_{1} = ϵη_{1} + ϵ^{2}η_{2} + ϵ^{3}η_{3} + O(ϵ^{4}) to get  (4.15)
−T_{1}sin(θ_{1}) + T_{2}sin(θ_{2}) =  (4.16)
−ϵT_{0}(\frac{1}{L_{1}} + \frac{1}{L_{2}} )η_{1} − ϵ^{2}T_{0}(\frac{1}{L_{1}} + \frac{1}{L_{2}} )η_{2} − ϵ^{3}T_{0}(\frac{1}{L_{1}} + \frac{1}{L_{2}} )η_{3}+  (4.17)
\frac{ϵ^{3}}{2} (\frac{T_{0} − k_{1}L_{1}}{L_{1}^{3}} + \frac{T_{0} − k_{2}L_{2}}{L_{2}^{3}} ) η_{1}^{3} + O(ϵ^{4})  (4.18)

35

[figure]

• The term in ϵ provides the linearised equation,

• the second equation provides η_{2} = 0

• and the term in ϵ^{3},

mη′′_{3} = T_{0}(\frac{1}{L_{1}} + \frac{1}{L_{2}} )η_{3} + \frac{1}{2} (\frac{T_{0} − k_{1}L_{1}}{L_{1}^{3}} + \frac{T_{0} − k_{2}L_{2}}{L_{2}^{3}} ) η_{1}^{3}  (4.19)

equation similar to what is obtained for the simplest mechanical example!

[figure]

4.2.3 A possible damage of a cable

is breakage of several fibers, this will cause decrease of rigidity k_{1} say for cable 1.

• Let us start with undamaged cables of same rigidity k. If we note L_{0}, the common length of the unstressed cables, and L their common stressed lenth, their tension is T_{0} = k(L − L_{0});

• now, after damage, k_{1} < k = k_{2}, cable 1 becomes longer and cable 2 shorter, L_{1} > L_{2}, the tension goes down to T_{00} = k_{1}(L_{1} − L_{0}) = k_{2}(L_{2} − L_{0});

• note the limit case of cable 1 broken is k_{1} = 0 so that the cable 2 gets lenth L_{0} but the system is no longer working properly!

• Before such a breakdown, if the change of tension is substantial, this causes a substantial change of the fundamental frequency; indeed, this is the routine monitoring of cable bridges!

• The nonlinear vibroacoustic testing aims at monitoring the cables before such a substantial change.

[figure]

4.2.4 Datas • L_{0} unstressed length,

• L half of the lenth of the span, or lenth of each of the stressed undamaged cables.

• k undamaged srping constant,

• from which “undamaged” tension T_{0} = k(L − L_{0}),

• L_{1} (with L_{0} < L_{1} < L) increased lenth of the damaged cable,

• from which, L_{2} = 2L − L_{1} decreased lenth of the undamged cable,

• from which “damaged”tension T_{0d} = k(L_{2} − L_{0}),

• from which spring constant of the damaged cable k_{1} = \frac{T_{od}}{L_{1}−L_{0}}

39

[figure]

Figure 14: ynl,f10,z,F12050

[figure]

Figure 15: ynl,f20,z,F16050

41

[figure]

5 A non linear string model

A model of non linear string has been introduced first by Kirchoff in 1877 and rederived by Carrier in 1945.

y_{tt} − T (∫_{0}^{l} y_{x}^{2})y_{xx} = f  (5.1)

For the classical linear string model, T is the tension of the string, assumed to be constant; in a next step, a natural asumption is:

T = T_{0} + k ∫_{0}^{l} y_{x}^{2}

it involves the linearized change of lenth as the length of the deformed string is:

l(y) = ∫_{0}^{l} \sqrt{1 + y_{x}^{2}}

[figure]

• Several mathematical studies of this type of equations have been performed recently (Medeiros(1994), Clark- Lima (1997).

• Following the lines of the discrete model, we intend to investigate a string made of two materials (safe and dameged).

• For a damaged string, k will be small on a small portion of the string:

T = T_{0} + k ∫_{0}^{d−ϵ} y_{x}^{2} + k_{d} ∫_{d−ϵ}^{d+ϵ} y_{x}^{2} + k ∫_{d+ϵ}^{l} y_{x}^{2}

43

[figure: 44]

[figure]

6 Masses free to move in a plane

Here, we assume that the masses can move freely; we denote:

\binom{X_{i}}{Y_{i}} the position at rest, \binom{x_{i}}{y_{i}} the curent position  (6.1)

• L_{i} lenth at rest; l_{i} lenth at time t; as the masses are moving freely:

l_{i}(x, y)^{2} = L_{i}^{2} + ((x_{i} − x_{i−1}) + (y_{i} − y_{i−1}))^{2} • and the change of tension of the linear elastic spring due to the change of of lenth T_{i} = T_{0,i} + k_{i}[l_{i}(x, y) − L_{i}] =. this tension is directed along the axis of the spring: T⃗_{i} = T_{i}τ⃗_{i}

[figure]

• Denote by θ_{i}, the angle of the spring with the horizontal axis, we have,

τ⃗_{i} = \binom{cosθ_{i}}{sinθ_{i}}

• y_{1} = x_{1}tan(θ_{1}), y_{i} − y_{i−1} = (x_{i} − x_{i−1})tan(θ_{i}) y_{n} − y_{n−1} = x_{n}tan(θ_{n}); but here, it is more convenient to use:

y_{i} − y_{i−1} = l_{i}(x, y)sinθ_{i}, x_{i} − x_{i−1} = l_{i}(x, y)cosθ_{i}

,

47

[figure]

Equation of the dynamics

m_{i}x_{i}” = −T_{i}cos(θ_{i}) + T_{i+1}cos(θ_{i+1}) + f_{i} i = 1 … n  (6.2)
m_{i}y_{i}” = −T_{i}sin(θ_{i}) + T_{i+1}sin(θ_{i+1}) + g_{i} i = 1 … n  (6.3)

We can express θ_{i} with respect to x_{i}, y_{i}, to obtain:

m_{i}x_{i}” = −T_{i} \frac{x_{i} − x_{i−1}}{l_{i}(x, y)} + T_{i+1} \frac{x_{i+1} − x_{i}}{l_{i+1}(x, y)} + f_{i} i = 1 … n  (6.4)
m_{i}y_{i}” = −T_{i} \frac{y_{i} − y_{i−1}}{l_{i}(x, y)} + T_{i+1} \frac{y_{i+1} − y_{i}}{l_{i+1}(x, y)} + g_{i} i = 1 … n  (6.5)

[figure]

7 Actively controled system, non destructive testing

• The case of an actively controled system is prospective; real experiments are not yet performed.

• Idea: to detect damage in real time taking adavantage of the data processed by the real time actuators used for the optimal control; real time control, research group: “Echtzeit Optimierung grosser Systeme”in Germany.

• Example of the vibrating masses: the forces u_{i} are now the control we consider the simple case of a quadratic functional:

F (u) = ∫_{0}^{t_{f}} (∑_{i} u_{i}^{2}(t)) dt

with final time conditions:

y_{i}(t_{f}) = 0, y′_{i}(t_{f}) = 0

• The initial conditions may be seen as a perturbation of the system, the active control brings to rest the system;

49

[figure]

• this process is supposed to be performed regularly during the lifetime of the system; in practice y_{i} is measured by sensors and the control u_{i} is a force performed by actuators; both devices transform electric energy in mechanical energy.

• the communication between both devices goes trough some computer

• If we are able to distinguish the response of a damaged system from an undamaged one, this opens the path of monitoring controled systems in real time as a dayly routine during their life.

• Numerical approach: to solve damaged and undamaged system and compare

• Perturbation approach, introduce a small parameter ϵ and expand the solution with respect to it; theoretical basis: the controled system should satisfy second order sufficient conditions (Malanowski, Maurer ...)

[figure]

Datas for an example of controlled 2 masses worked out by by K. Theissen (U. Muenster)

\begin{matrix}T_{0,1} = T_{0,2} = T_{0,3} = & 1 \\ k_{1} = k_{2} = k_{3} = & 5 \\ m_{1} = m_{2} = & 1 \\ L_{1} = L_{2} = L_{3} = & 1 \\ t_{f} & 100 \end{matrix}

51

[figure]

Figure 16: Frequences of u_{1}

[figure]

Figure 17: Frequences of u_{2}

53

53-1

[figure]

8 Bar models with defects

Bar models with longitudinal waves (dynamical traction and compression) are considered.

ρ\frac{∂^{2}u}{∂t^{2}} − \frac{∂n}{∂x} = f(x, t)  (8.1)

With a non linear stress-strain law:

n = E(A\frac{∂u}{∂x} + ϵχ_{[a,b]}(\frac{∂u}{∂x})^{3})  (8.2)

Also a linear law is considered with a modified equation::

n = EA\frac{∂u}{∂x}  (8.3)

53

[figure]

it may correspond to the action of a non linear spring acting on part of the bar :

ρ\frac{∂^{2}u}{∂t^{2}} − \frac{∂n}{∂x} + ϵχ_{[a,b]}u^{3} = f(x, t)  (8.4)

We could as well assume that the applied load is of order epsilon without any assumption on the nonlinearity. Assuming ϵ to be small an approximate solution is searched for with the following ”ansatz”:

u = u_{0} + ϵu_{1} + … d’où  (8.5)
u^{3} = u_{0}^{3} + 3ϵu_{0}^{2}u_{1} + …  (8.6)
\frac{∂u}{∂x}^{3} = \frac{∂u_{0}}{∂x}^{3} + 3ϵ\frac{∂u_{0}}{∂x}^{2} \frac{∂u_{1}}{∂x} + …  (8.7)

(8.8)

[figure]

From which we get for the non linear law:

n = E(A\frac{∂u_{0}}{∂x} + ϵ (A\frac{∂u_{1}}{∂x} + χ_{[a,b]}(\frac{∂u_{0}}{∂x})^{3}) + …  (8.9)

and for the linear law:

n = EA(\frac{∂u_{0}}{∂x} + ϵ\frac{∂u_{1}}{∂x}) + …  (8.10)

Using these expansions, with the non linear law, the following system is obtained:

55

[figure]

ρ\frac{∂^{2}u_{0}}{∂t^{2}} − EA\frac{∂^{2}u_{0}}{∂x^{2}} = f(x, t)  (8.11)
ρ\frac{∂^{2}u_{1}}{∂t^{2}} − EA\frac{∂^{2}u_{1}}{∂x^{2}} = −E \frac{∂}{∂x}(\frac{∂u_{0}}{∂x})^{3}χ_{a,b}  (8.12)

For the modified equation the same equation for u_{0} is found but for u_{1}:

ρ\frac{∂^{2}u_{1}}{∂t^{2}} − EA\frac{∂^{2}u_{1}}{∂x^{2}} = −(u_{0})^{3}χ_{[a,b]}  (8.13)

[figure]

Theoretical justification of the expansions:

Non liner law The situation is complex in full generality: non linear hyperbolic equations exhibit a singularity after a finite time! But: the experiments are performed during a short time interval and the Fourier transforms are computed on these time intervals! Following a suggestion of Guy Metivier we are addressing the problem during a small initial time interval in which the solution is smooth: plan to use an approximation of the equation with a fixed point method proposed in Majda. In any case we should smooth the characteristic function (the material is changing smoothly)!

Modified equation The situation is simpler; we can use a priori inequalities for this type of equation.

57

[figure]

8.1 Explicit Solution

Coefficients are assumed to be constant and we consider:

Clamped at both ends: u(x, 0) = 0 = u(x, l); Eigenfunctions are introduced :

EA\frac{∂^{2}ϕ}{∂x^{2}} = −λρϕ  (8.14)
ϕ(0) = 0 = ϕ(l)  (8.15)

we find λ_{k} = \frac{k^{2}π^{2}}{l^{2}} \frac{EA}{ρ}, on pose ω_{k} = \sqrt{λ_{k}} and the normalised eigenfunction:

ϕ_{k} = \sqrt{\frac{2}{l}}sin(\frac{kπ}{l} x).

[figure]

8.1.1 Computation of u_{0} let us consider a force of f́requency \frac{α}{2π}

f(x, t) = F cos(αt)sin(\frac{kπ}{l} x)  (8.16)

with initial velocity: \frac{∂u}{∂t}(x, 0) = 0 The solution

u_{0} = \frac{F cos(αt)}{ρ(−α^{2} + λ_{k})}sin(\frac{kπ}{l} x)  (8.17)

corresponding to an initial condition

u_{0}(x, 0) = \frac{F}{ρ(−α^{2} + λ_{k})}sin(\frac{kπ}{l} x)\frac{∂u_{0}}{∂t} (x, 0)  (8.18)

59

[figure]

For the initial condition:

u_{0}(x, 0) = a_{0}sin(\frac{kπ}{l} x)  (8.19)

the solution is:

u_{0}(x, 0) = [\frac{F}{ρ(−α^{2} + λ_{k})}(cos(αt) − cos(ω_{k}t)) + a_{0}cos(ω_{k}t)] sin(\frac{kπ}{l} x)  (8.20)

[figure]

8.1.2 Computation of u_{1} Considering the first solution with a global non linearity, we get:

u_{0}^{3} = \frac{cos(αt)^{3}}{ρ^{3}(−α^{2} + λ_{k})^{3}} sin^{3}(\frac{kπ}{l} x) =  (8.21)
\frac{cos(αt)^{3}}{ρ(−α^{2} + λ_{k})^{3}} [cos(3αt)sin(3\frac{kπx}{l}) − 3cos(3αt)sin(\frac{kπx}{l})  (8.22)
+3cos(αt)sin(\frac{3kπx}{l}) − 9cos(αt)sin(\frac{kπx}{l})]  (8.23)
\frac{∂u_{0}}{∂x}^{3} = \frac{k^{3}π^{3}}{l^{3}} u_{0}^{3}; \frac{∂}{∂x} \frac{∂u_{0}}{∂x}^{3} = \frac{k^{3}π^{3}}{l^{3}} \frac{∂u_{0}^{3}}{∂x} = (8.24)

(8.25)

solution u_{1} with frf́requency \frac{3α}{2π} or \frac{2α}{2π} for a quadratic non linearity.

61

\begin{matrix}\frac{cos(αt)^{3}}{ρ^{3}(−α^{2} + λ_{k})^{3}} \frac{k^{4}π^{4}}{l^{4}} [3cos(3αt)cos(3\frac{kπx}{l}) + 3cos(3αt)cos(\frac{kπx}{l}) & (8.26) \\ +9cos(αt)cos(\frac{3kπx}{l}) − 9cos(αt)cos(\frac{kπx}{l})] & (8.27) \end{matrix} ; 61-1

[figure]

Second case for the second pair of boundary conditions, we set:

c = \frac{F}{ρ(−α^{2} + λ_{k})} d = (−\frac{F}{ρ(−α^{2} + λ_{k})} + a_{0})  (8.28)

62

[figure]

Now we have:

u_{0} = (c cos(αt) + d cos(ω_{k}t)) sin(\frac{kπx}{l})  (8.29)

(u_{0})^{3} = [\frac{c^{3}}{4} cos(3αt) + \frac{3c}{2}(\frac{c^{2}}{2} + d^{2})cos(αt)+

\frac{3c^{2}d}{4}(cos((ω_{k} + 2α)t) + cos((ω_{k} − 2α)t))+

\frac{3cd^{2}}{4} (cos((2ω_{k} + α)t) + cos((2ω_{k} − α)t)+)

+\frac{3d}{2}(\frac{d^{2}}{2} + c^{2})cos(ω_{k}t)\frac{d^{3}}{4} cos(3ω_{k}t)]

\frac{1}{4} (3sin(\frac{kπx}{l}) − sin(\frac{3kπx}{l}))  (8.30)

(8.31)

[figure: \begin{matrix}\frac{∂}{∂x} (\frac{∂u_{0}}{∂x} )^{3} = \frac{k^{3}π^{3}}{l^{3}} \frac{∂u_{0}^{3}}{∂x} = & (8.32) \\ \frac{3k^{4}π^{4}}{4l^{4}} [\frac{c^{3}}{4} cos(3αt) + \frac{3c}{2}(\frac{c^{2}}{2} + d^{2})cos(αt)+ \\ \frac{3c^{2}d}{4}(cos((ω_{k} + 2α)t) + cos((ω_{k} − 2α)t))+ \\ \frac{3cd^{2}}{4} (cos((2ω_{k} + α)t) + cos((2ω_{k} − α)t)+) \\ +\frac{3d}{2}(\frac{d^{2}}{2} + c^{2})cos(ω_{k}t)\frac{d^{3}}{4} cos(3ω_{k}t)] \\ (cos(\frac{kπx}{l}) − cos(\frac{3kπx}{l})) & (8.33) \end{matrix}]_{64}

[figure]

We notice clearly terms of frequency \frac{α}{2π} and \frac{3α}{2π} but also cross-modulations: \frac{ω_{k}+2α}{2π} et \frac{2ω_{k}+α}{2π} and frequencies \frac{3ω_{k}}{2π} \frac{ω_{k}}{2π}. This last term provides secular terms for the corrector term u_{1}; they ought to be eliminated for example by using some renormalization technique:

t = s(1 + ϵω_{1} + …)  (8.34)

We notice that the perturbation is larger if α is close to ω_{k}. this fact is used in practice: the applied load uses two frequencies with the low one at the first resonance in [8]. Here the low frequency is excited by the initial conditions.

[figure]

9 Conclusion

• Some simple models governed by ODE pr PDE show intermodulations;

• But what is the relative level of secondary peaks for a given set of datas deserves investigations: indeed it is also the difficulty of the real experiments

• Need to include other behaviors: shocks, friction

• Need of more precise models: non linear beams including tractional, flexural, torsional effects

• Mixture of local models for the defect and global models for the undamaged structure to obtain precise results at low computational cost.

66

[figure]

References

[1] D. Donskoy A. Sutin. Nonlinear vibro-acoustic nondestructive testing technique. In Nondestructive characterisation of material, 7 Ed R.E. green. Plenum press, New York, 1998, ?

[2] V.V.Kasakov A.E.Ekimov, I.N.Didenkulov. Modulation of torsional waves in a rod with a crack. J.Acoust. Soc. AM., 3(106):1289–1291, 1999.

[3] N. Sedunov Tsionskiy D. Donskoy, A. Ekimov. Nonlinear seismo-acoustic land mine detection and discrimination. . Acoust. Soc. Am, 111:2705–2714, 2002.

[4] N. Damil M. Potier-Ferry E.H. Boutyour, B. Cochelin. Calculs non linéaires par des méthodes asymptotiques-numériques: applications aux structures élastiques. In Colloque national en calcul de structures, 11-14 mai 1993. Hermes, 1993.

[5] A. Moussatov-B. Castagnede-V. Gusev. Frequency up-conversion and frequency down-conversion of acoustic waves in damaged materials. Physics letter A, 301:281–290, 2002.

[6] B.P. Lathy. Signals, systems and communications. J. Wiley and sons, 1965.

[7] G. Vanderborck M. Lagier. Nonlinear acoustic spectroscopy for crack detection in

[figure]

solids. JASA, to be submited.

[8] M. Lagier P. Tèmin G. Vanderborck P. Dufourcq, JP. Groby. Détection vibro-acoustique non linéaire d’ endomagements dans une structure poutre. Communication au Congrès français de mécanique, septembre 2003.

[9] P. Sas V. Zaitsev. Nonlinear vibro-acoustic response of a metal sample with a discontinuity like defect as related to damage detection problems. In Proceedings of DECT 99, Las Vegas, Nevada, 1999.

[10] Lagier Michel Vanderborck Gerard. Application of non-linear ultrasonic spectroscopy to health monitoring and damage detection in structures,. 38p. In 75th Shock and Vibration Symposium, Virginia Beach (VA) USA, du 18/10/2004 au 21/10/2004.

[11] C. Gasquet-P. Witomski. Analyse de Fourier et applications. Masson, 1990.