Docteur En

Andrey Trifonov

Président Pr

Emmanuel Moulin

Dr Cédric Payan

Mdc Hdr

Pr Koen

Van Den Abeele

Dr Sigrun Membre

Hirsekorn

Dr Anissa Meziane

Dr Steven Delrue

Dr Membre Invité

Vladislav Aleshin

Olivier Bou Matar -Lacaze

CENTRALE LILLE THESE Présentée en vue d'obtenir le grade de

non linéaire pour des applications aéronautiques / Nonlinear acoustic nondestructive testing for aeronautical applications

technique is in its extremely high sensitivity owing to the fact that weak changes in sample parameters are accumulated and finally greatly amplified during the formation of the coda wave i.e. multiple scattered late arriving part of the reverberation signal. The other technique has a complimentary strength and offers a possibility of a remote detection. More precisely, the processes of signal generation and detection are totally remote while the pumping wave excitation is still performed in a contact manner. The developed techniques are tested on samples with artificially fabricated defects at known locations. The performance of each method is accessed and the potential for obtaining robust nonlinear images is demonstrated.

The second (theoretical) part of the work is concerned with a theoretical description of contact acoustical nonlinearity and its use for creating of a numerical toolbox capable of simulating wave propagation in complex structures containing internal contacts. A physical model describing the tangential shift of two contacting bodies in the presence of friction has been proposed. Its result is an analytical computer-assisted solution for hysteretic relationships between normal and tangential contact displacements and loads. The contact model and derived load-displacement relationships are used as boundary conditions posed at the internal boundaries (contact surfaces) in a finite element wave propagation model programmed via existing software (COMSOL). This work has been carried out in the framework of the ALAMSA European project.

Résumé en français

Contexte

Ce travail de thèse est une contribution au développement des méthodes d'acoustique non linéaire pour le contrôle non destructif et l'imagerie de défauts de type contact dans les solides.

Ces techniques reposent sur le fait que la présence dans un solide de contacts internes (des défauts) modifie considérablement la propagation acoustique en créant des non linéarités acoustique de contact qui à leur tout génèrent des signatures non linéaires spécifiques. Ces signatures peuvent être mesurées en utilisant des techniques appropriées et ainsi utilisées afin de détecter et caractériser des défauts.

L'étude présenté dans ce manuscrit a été réalisé dans le cadre du projet Européen ALAMSA ("A Life-cycle Autonomous Modular System for Aircraft Material State Evaluation and Restoring System").

Depuis plusieurs années il est connu que les méthodes de contrôle non destructif (CND) par acoustique non linéaire permettent une détection fiable et précise de l'endommagement dans divers matériaux et structures. Le CND par acoustique non linéaire utilise maintenant toute une gamme de techniques développées pour des applications spécifiques en faisant face à des exigences particulières. Dans notre cas, ces exigences portent essentiellement sur deux aspects : la méthode doit être suffisamment fiable et sensible afin de détecter de faibles non linéarités, et en même temps capable de fonctionner en environnement industriel. D'un autre côté, une étude purement expérimentale en CND est fréquemment insuffisante. En fait, les signatures non linéaires ne caractérisent généralement pas directement le matériau et les défauts qu'il contient. Il est souvent nécessaire de fournir un modèle permettant d'interpréter le lien entre les propriétés du matériau mesurables par CND non linéaire et les paramètres du défaut (une fissure par exemple). Ce modèle doit pouvoir s'appliquer à des structures de géométrie complexe et doit reposer sur la physique des contacts internes afin de garantir sa fiabilité.

Les objectifs de cette thèse étaient donc :

Non linéarité acoustique de contact et contrôle non destructif

Le manuscrit débute par rapide état de l'art des différents types de non linéarités mécaniques, à savoir les non linéarités géométriques, les non linéarités de l'équation d'état du matériau Les techniques de CND par acoustique non linéaire développées récemment dans le cadre du projet ALAMSA sont ensuite présentées : La tomographie d'ondes guidées utilisant l'algorithme RAPID, le principe de résonance localisée de défaut (LDR) permettant d'accroitre considérablement l'amplitude de réponse du défaut, et utilisé conjointement avec la vibrométrie laser à balayage, la thermographie ultrasonore, et la shearographie.

Interférométrie non linéaire d'ondes de coda

Une nouvelle technique permettant d'obtenir simultanément une information globale sur l'état de santé d'un échantillon et de localiser un défaut a été développée. Cette technique combine l'interférométrie de coda non linéaire (CWI) et le principe de focalisation par retournement temporel (TR).

L'interférométrie d'ondes de coda est une technique permettant de mesurer des décalages temporels à l'aide des variations de phase et d'amplitude d'ondes de coda, qui sont les queues des signaux reçus lorsque les ondes se propagent dans un milieu multi-diffusant ou réverbérant. Dans cette technique le milieu se comporte comme un interféromètre « naturel ».

Ainsi, des variations extrêmement faibles ou localisées des propriétés du matériau ou de la géométrie de l'échantillon produisent des variations mesurables des ondes de coda. Afin de La technique de localisation combinant la technique CWI et le retournement temporel est ensuite testée sur une plaque de verre rectangulaire (30 x 50 x 1 cm) contenant une fissure induite par un impact réalisé à l'aide d'un pistolet pneumatique (Fig. 6). Le protocole de mesure est décrit sur la Fig. 7a. La mesure CWI est synchronisée avec le processus de focalisation par la cavité chaotique afin de mesurer les variations de  et de Kd dans une fenêtre temporelle où l'onde de coda est en interaction avec l'onde de pompe.

Probe emitter Pump emitter Probe receiver Fig. 6. Plaque de verre rectangulaire (30 x 50 x 1 cm) contenant une fissure induite par un impact réalisé à l'aide d'un pistolet pneumatique.

de pompe n'est pas focalisée (en vert), le cas où elle est focalisée sur le défaut (en noir), ou en dehors du défaut (en bleu).

No changes are visible in the stretching parameter when the pump wave is focused or not on the defect. Conversely, Kd's changes demonstrate the possibility to localize the defect using CWI. Des mesures ont aussi été réalisées sur un échantillon de GFRP de 4 mm d'épaisseur contenant une délamination artificielle plus réaliste (sans aucune inclusion de Teflon ou de Kapton). Le résultat de la mesure non linéaire obtenu en excitant l'échantillon à une fréquence LDR voisine de 56 kHz est comparé avec une image obtenue par une mesure par ultrasons aériens (300-700 kHz) en transmission sur la Fig. 10.

Modèles de contact avec friction

Les modèles de contact jouent un rôle primordial dans la description des propriétés de matériaux contenant des contacts internes. La littérature sur le sujet débute en 1880 lorsque H.

Hertz ) publie la solution, maintenant classique, de deux sphères élastiques

      da b T dN     End of current step Take new N=N n -N p , T=T n -T p , a=a(N n )-a(N p ) Check if N>0 Check if |T|/(N)<1 No   1 sgn T T N      2 1 T T T          sgn , n n D T Q N         2 p N Q T D d          1 sgn b T a        2 p N Q da b D d dN         1 2 b b b      No   1 n p N p N T D d       2 1 T T T          2 sgn , n n D T Q N         2 n N Q T D d          1 n p N p N da b D d dN          2 p N Q da b D d dN         1 2 b b b     

Yes

20

Deux exemples de résultat obtenus à l'aide du modèle pour l'évolution du déplacement tangentiel en fonction de la force tangentiel pour deux chargements temporels différents présentés en insert sont présentés sur la Fig. 13.

-2 -1 0 1 2 3 -2 -1 0 1 -2 0 2 4 6 -2 -1 0 1 2 3 -2 -1 0 1 2 -2 0 2 4 (a) b/(a 0 ) (b) b/(a 0 ) T/(N 0 ) T/(N 0 ) T/(N0) T/(N0) N/N0 N/N0

Modélisation de la propagation d'ondes élastiques dans un matériau contenant des fissures

Toute simulation de la propagation d'ondes élastiques dans un matériau contenant des fissures nécessite un modèle réaliste de fissure. Un tel modèle doit être formulé pour un petit volume de matériau incluant la fissure, c'est-à-dire à une échelle mésoscopique (petite devant la longueur d'onde mais grande devant les aspérités de la fissure) à laquelle des conditions aux limites sur les surfaces de la fissure doivent être définies. Le modèle de fissure doit donc fournir les conditions aux limites représentant le lien entre les forces et les déplacements à la surface de la fissure. D'un point de vue numérique, le modèle de fissure doit échanger des données de type force-déplacement avec un module de mécanique du solide qui résout les équations de l'élasticité dans le volume du matériau. Le modèle de fissure et le module de mécanique du solide peuvent tous deux être pilotés par les forces ou par les déplacements.

Afin d'obtenir dans le cas où la friction est prise en compte un échange de données explicite nous avons choisi d'utiliser un modèle de fissure piloté par les déplacement (voir Fig. 14).

CRACK MODEL

Normal

SOLID MECHANICS

Input

, ab Output , NT , NT , ab   N N a    , T T b a  Tangential

Introduction and objectives

Context

Generally, this PhD thesis work is related to nondestructive testing based on contact acoustical nonlinearity. The principle underlying this research is the fact that the presence of internal contacts (defects) in solids considerably modifies acoustic propagation and results in contact acoustical nonlinearity which, in turn, generates a number of detectable nonlinear signatures. These signatures can be accurately measured by using appropriate techniques and then used for interpreting the material properties and characterizing the defects.

The study described here was performed in the framework of a European Project entitled "A

Life-cycle Autonomous Modular System for Aircraft Material State Evaluation and Restoring

System" (ALAMSA). The four-years project (2012-2016) funded by the Framework 7

Programme of the European Union brought together the efforts of ten European academic, research and industrial partners. Its main objective was the creation of an innovative selfrestoring system for aircraft materials. Physically, the self-healing principle is based on the integration into composite materials of breakable fibers filled with chemical reagents. Any event that results in breaking the fibers liberates the reagents whose chemical reaction creates a rigid agent that solidifies the damaged material.

For many years it has been known that nonlinear acoustic non-destructive testing (NDT) is capable of robust and precise detection of damage in various materials and structures and therefore suggests an opportunity to test the final efficiency of a self-healing process. The nonlinear acoustic NDT now uses a whole range of techniques, each developed for specific applications in order to meet some particular requirements. In our case, these requirements include at least two essential aspects: the method should be sufficiently robust and sensitive in order to detect weak nonlinearities and, at the same time, it should be suitable for using in real field but not laboratory conditions. In fact, these requirements are related to each other; indeed, the necessity to remotely test a complex structure with a number of geometric features often makes the measurable signatures weak even when the actual damage is strong.

Generally, an attempt to apply the developed techniques in real field greatly increases the requirements to robustness and sensitivity. Ideally, the goal should not only be to develop a damage detection technique but an imaging method capable not only of detecting damage located somewhere but of its localization as well.

On the other hand, purely experimental investigations in NDT are frequently not sufficient.

The matter is that a nonlinear signature by itself does not characterize the materials and the damage in there directly. In any case, some interpretation that links the measurable properties and the actual damage parameters is of interest. Certainly, such an interpretation should be applicable for structures of complex geometries and, in addition, should be based on physics of internal contacts otherwise its efficiency is not guaranteed.

Objectives

Summarizing these desired requirements we formulate the following objectives of the present study:

Experimental

Geometric, material and contact nonlinearities

Section 1 contains a brief review of mechanical nonlinearity types, namely geometric, material, and contact. The geometric nonlinearity is related to a change of a distance between two points under loading. The origin of material nonlinearity is in the inharmonicity of interatomic potentials. Contact nonlinearity appears when the material contains internal contacts such as cracks, delaminations, imperfect intergranular boundaries etc.

Geometric nonlinearity

Geometric nonlinearity appears when one attempts to link a change of distance between two points in a strained material with the displacement vector. This link is given by expression (1):

22 '2 ik i k dl dl dx dx   , (1)
where dl and ' dl are distances between close points in unstrained and strained bodies,

respectively, ) 3 , 2 , 1 (  i x i
are Cartesian coordinates, and ik  is the strain tensor calculated as 1 2

jj ik ik k i i k uu uu x x x x              , ( 2
)
through the displacement vector i u .

Usually in mechanics and acoustics all strains are small, i.e. changes in the sample size are much less than the size itself. This means that all derivatives of the kind / ik ux  are much less than 1. Hence, the product of the derivatives in Eq. ( 2) can be neglected and

1 () 2 ik ik ki uu xx     . ( 3
)
The neglected term represents the geometric nonlinearity. An attempt to keep it produces nonlinear equations of solid mechanics . Account for geometric nonlinearity is essential only for specific systems, e.g. bending of thin bars or thin plates, cable-stayed structures, etc.

Material nonlinearity

The origin of the material nonlinearity is related to the potential energy of interactions between atoms forming the lattice. In a perfectly linear material, the interatomic potential would have a parabolic (harmonic) shape. At the same time, a more realistic approximation is the Lennard-Jones potential [Len-24] describing the interaction between a pair of neutral atoms or molecules:

12 6 0 ( ) 2 mm LJ rr U r U rr                  , (4)
illustrated in Fig. 1. Here, 0

U is depth of the potential well, m r is the distance at which the potential reaches its minimum. It is straightforward to show that this potential demonstrates a weak inharnonicity that starts to be more important when r deviates considerably from the rest radius m r :

  2 3 4 0 1 36 252 1113 m m m LJ r r r r r r U r U r r r                              
.

(5)

The presence of cubic and fourth-order terms is the reason of material nonlinearity. Indeed,

the value 2 (( ) / ) m r r r 
is related to strain while the derivative / LJ dU dr characterizes stress, therefore the non-quadratic terms in Eq. ( 5) produce a nonlinear stress-strain relationship.    ) invariants in which the sum over indexes is meant. Hence in the isotropic case

2 2 2 3 2 3 3 3 ik ll ik il kl ll ik ll K A C WB                   , (6)
where K and  are the bulk and shear moduli, respectively. This fact can be confirmed by differentiation of the first two components of W over its arguments ik  which produces Hooke's law for isotropic materials:

1 2 3 ik ll ik ik ik ll K              , (7)
where ik  is the stress tensor and ik  is the Kronecker delta.

In Eq. ( 6) three nonlinear constants A, B and C have been introduced. Thus the full constitutive model of material nonlinearity is formulated using five constants only, and the stress-strain relationship is given. Description of anisotropic materials requires more material constants, but the stress-strain relation is still written as a closed-form expression. As it will be shown in the next section, the contact nonlinearity case is generally more complex.

Contact nonlinearity

Contact nonlinearity is the third class of mechanical nonlinearities considered here. The former two classes are related to uniform materials containing no defects in their structure. At the same time, pure, uniform and regular materials are exceptional in common life. 99% of the time mankind produces, treats and uses materials which have impurities, irregularities, inclusions, defects etc., which are inherent properties of their microstructure. Internal defects can be roughly categorized in 3 types: 1D dislocations, 2D internal contacts and 3D pores, voids, etc. Amongst these, the second type is the most essential in terms of material performance, since the presence of internal contacts manifests itself in the most drastic way.

Indeed, influence of dislocations is negligible if we speak about seismology or building constructions, whereas pores and voids usually contribute to the most interesting material properties much less than cracks and contacts do (e.g., failure loads, acoustic and static nonlinearities, sound attenuation, etc). This makes solids with internal contacts to be an extremely important class of materials, and justifies the fact that an accurate description of their mechanical properties is critical.

Most typical examples of materials with internal contacts are unconsolidated granular materials in which the only physical link between the constituents (grains) is through internal contacts, and consolidated materials in which there exist a solid matrix whose properties are modified by the presence of contacts. In the latter class, two groups can be distinguished. In some materials, internal contacts are present as an inherent part of their structure (consolidated grainy materials such as geomaterials or building construction materials).

Generally speaking, all solids that are not single crystals can be regarded as materials with inherent random structure at a mesoscopic scale, i.e., a scale which significantly exceeds the atomic size but is still small compared to macroscopic dimensions. Finally, there are solids in which internal contacts appear as defects (cracks, delaminations, etc.). Studies for materials of this class form a basis for theories underlying nonlinear NDT techniques.

Contact nonlinearities appear even when the material itself is perfectly linear. An obvious reason for that is the fact that the contact can be open or closed. In the former case the faces do not interact while in the latter one there is an interaction. This effect alone results in bimodality i.e. dependence of elastic moduli in a material on the states of contacts. Another nonlinear mechanics is related to the contact geometry. In most cases, contacting faces have some profiles, including regular shapes or random topographies such as roughness. The simplest example is Hertzian spheres. Even for perfectly elastic materials and for spheres always staying at contact (no bimodality effect), the force-displacement relationship in such a system is nonlinear. The issue is that higher displacements involve deformations of deeper layers of the material, whereas in the linear case the 'amount' of strained material stays the same.

In the next sections we consider existing models for contact nonlinearity in more detail.

Existing models for contact nonlinearity

For modeling the nonlinear elastic behavior of materials, there exist a vast set of models. Here we concentrate only on theories capable of producing nonlinear stress-strain relationships, which is of primary importance for building up numerical models imitating nonlinear wave propagation, for the final purpose of creating a numerical tool for nonlinear NDT applications.

Besides, there are many models that predict other characteristics, such as wave dispersion and attenuation, slow dependence of parameters on time, modification of linear elastic properties in the presence of damage, etc. In addition to the classical Landau theory, there exist two classes of so-called "nonclassical" models: phenomenological, in which some desired behavior is directly postulated as a simple or complex stress-strain relation, and physical, when it is attempted to take into account the physical comportment of real internal contacts.

Most of the models target a 1D case only and are primarily suitable for simple wave propagation geometries.

Phenomenological models

Clapping or bimodal model

The simplest phenomenological model addresses a single contact (crack) perpendicular to the direction of plane longitudinal wave propagation. If the incident wave has a stress amplitude that exceeds the static stress of the originally closed interface, it opens the crack which results in a change of stiffness of the whole material. The simplest way of modeling this situation is to assume two different stiffnesses for positive and negative strains:

  1 C CH C         , (8)

Sliding friction model

Whilst the previous model addresses the case of normal wave incidence, the sliding friction model is related to the tangential wave-to-interface interaction. Consider the non-bounded interface between two friction-coupled surfaces subjected to an oscillating tangential traction (shear wave scattering) strong enough to cause their sliding. Then suppose that gross sliding of the interfaces occurs when the shear wave stress ε exceeds certain a value ε 1 . Then, the tangential stiffness, which has a value of C in the stick phase, drops to zero in the sliding phase, so that (9)

(see ). It is essential to note that, whereas in the clapping model both odd and even harmonics are generated, the sliding friction model is symmetric to positive and negative thus generating odd harmonics only.

These two basic mechanisms of contact nonlinearity are often present in a mixed form. On the elementary level, they produce a plausible explanation of the effects encountered. Moreover, if the geometry of the sample makes one of these two modes preferable, this can be seen in the presence/absence of even harmonics.

Preisach model

The In this formalism, the Preisach density represents the portrait of the system. By varying the density, the responses of a wide range of hysteretic systems can be imitated. The method allows one to obtain the response y(t) of a hysteretic system for any signal x(t) that makes it particularly suitable for complex acoustical excitations. There are two identifying criterions for the Preisach system (see [May-85]); using these criterions it is possible to check if a particular hysteretic model represents the Preisach system or not.

In addition, there exist procedures (see [Guy-97], [Ale-08]) for reconstructing the Preisach density from a sole particular hysteretic curve. The methods work well for high quasi-static loads where up-going and down-going branches of a hysteresis loop are clearly distinguished.

However, many researchers use the Preisach formalism for acoustical simulations as well Vn assuming a constant density or accepting some model density forms. The advantage of the approach is in its flexibility and possibility to imitate a wide class of hysteretic dependencies.

The disadvantage is related to its phenomenological character; the theory becomes physical only when the Preisach density and other model characteristics are linked to physical and geometrical parameters of the material and damage. Assuming various forms for the Preisach density, it is possible to deduce (see [ ) from the Preisach formalism a number of particular hysteretic models, such as the hysteretic quadratic nonlinearity, Nazarov, Davidenkov, and Granato-Lücke models.

Hysteretic quadratic nonlinearity

This model is deduced from the Preisach formalism by assuming and harmonic excitation with the amplitude ε, with the resulting expression for stiffness:

(10)

where  and  are material constants. Thanks to its simplicity, the model was used by a number of researchers to analyze the impact of hysteresis on the wave shapes (see [Ale-04]), to simulate resonance experiments, and to reveal theoretical amplitude dependencies for the nonlinear resonant frequency shift and the higher harmonics in resonance (see [Van-00]), etc.

Nazarov model

The model (see [Naz-03]) is applicable for the harmonic excitation of amplitude ε; it has three parameters and describes the up-doing and down-going branches separately:

            2 2 1 2 1 2 2 1 2 2 42 42 C C                                                             (11)
The model is a modification of the model that uses hysteretic quadratic nonlinearity; it can be obtained from the Preisach formalism by setting: (12) where ε o and ε c are strain-related arguments introduced instead of formal arguments x o and x c .

Davidenkov model

The adapted ) assumes the following expressions for the stress-strain curves:

                1 1 2 2 mm m mm m C m C m                                           (13)
it can be retrieved from the Preisach formalism by setting

     2 ,1 m o c o c Cm           . ( 14
)

Granato-Lücke model

The Granato-Lücke model (see [Gra-56], [Naz-09]) originates from physics of dislocations which also involves hysteretic behavior. The stress-strain curves according to the (adapted)

Granato-Lücke model contain four branches which can be expressed as follows:

1 1 2 1 2 3 1 4 3 4 ( ) 0, 0 ( ( ) ) 0, 0 ( , ) ( ) 0, 0 ( ( ) ) 0, 0 m mm m mm if m if mm C if m if mm                                                                (15)
Here the Preisach space is not constant (quadratic hysteretic, Nazarov) and does not depend only on the perpendicular coordinate

  /2 oc      (Davidenkov), but equals       2 1 12 2 1 34 0 ,0 0 m m o o c m m m o c c c o m C if C if otherwise                                (16)
Since the density equals 0 in a whole quadrant ε c <0, ε o >0, some internal loops within the large hysteresis loop can be actually non-hysteretic (reversible i.e. no difference between up-going and down-going branches). The original Granato-Lücke model does not have this property since it is not intended for non-harmonic signals at all.

Compared to the generic Preisach formalism, these models have the advantage of having less free parameters. However, they are still purely phenomenological and have an ad hoc character.

Vector and tensor Preisach models

In the vector Preisach model, it is suggested (see [May-88]) that any direction of the output vector generates its own scalar Preisach space that depends on two switching values (x o , x c ).

The elementary hysteretic units are switched by the projection of the input vector on the direction of the output vector. As such, the Preisach density additionally depends on the polar angle in 2D case, =(x o ,x c ,), and on the spherical angles, =(x o ,x c ,,), in the 3D case.

This also means that the distribution of open and closed elements can be different for each angle.

The tensorial model (see [Hel-01]) uses the Kelvin notations and representation of the stress and strain tensors as 6-component vectors. These vectors are linked with a 6x6 matrix that can be diagonalized thus presenting 6 eigenvectors. It is suggested that these 6 eigenvectors are separately multiplied by the stress vector which results in the appearance of 6 scalar quantities.

Then, 6 independent Preisach spaces are introduced with the above 6 scalar inputs; the corresponding outputs are considered as weighting coefficients for the 6 eigenvectors in a linear combination that represents the global output of the model.

The tensor Preisach model actually implied some additional simplifications compared to the vector model. They come from the imposed consideration of 6 eigenvectors in the Kelvin notations that brings into the model only 6 independent Preisach spaces, while in the vector model the number of Preisach spaces is continual. In any case, the number of free parameters is huge which makes the identification problem (retrieving material's parameters from data) hardly solvable.

Preisach-Arrhenius model

This model (see [Tor-02], [Gus-05]) is one of the different possible extensions of the Preisach formalism in which energy barriers separating the open and closed states of the hysteretic elements are introduced. These barriers are suggested to be very weak (comparable to kT). This modification produces a number of dynamic effects, such as dispersion, spontaneous transitions between the two metastable states, etc. It can be used when the scalar modeling of the stress-strain relationship is sufficient, and dynamic effects and relaxation are necessary to be taken into account.

In conclusion it can be noted that the models of the phenomenological class are suitable only for imitating contact nonlinearity (in a precise or qualitative way), but neither for explaining it nor linking it to the materials' parameters. They can be very flexible, but for the price of high number of free parameters and associated difficulties in solving the identification problem. At the same time, the simplified versions can reproduce the nonlinear effects in a qualitative manner only.

Physical models

Physical models use an approach which fundamentally differs from the pure phenomenological considerations. In this case it is assumed that there exists one or more physical mechanisms responsible for the nonlinear effects. Accepting a certain geometry (even greatly simplified) of the material and of the defects, we then select an appropriate model for the basic physical mechanism (friction, adhesion, collective transformations of dislocations, etc.) hypothetically underlying the nonlinear behavior. Below we exemplify several models of this class.

Lawn and Marshall model

In this model (for the original version, see [Law-98]; the modified model taking into account additional effects is presented by [Ale-07a]), an elastic continuum contains a large number of diversely oriented plane cracks. For each of them, the Coulomb friction law is postulated that assumes either stick state or gross sliding of the crack surfaces. The choice of the state depends on the normal and shear stresses transmitted to the crack faces from distant boundaries of the sample where external loading conditions are posed. Then, all displacements of the cracks faces (only tangential in the original model of 1998 and both tangential and normal in the modified one of 2007) are summed up with the relevant orientational weighting coefficients to produce the total strain tensor. It is seen that a large number of cracks having the postulated behavior give rise to hysteretic responses of the Preisach type, with internal loops for partial loading-unloading-reloading, etc. In other words, the model explains the macroscopic stress-strain hysteresis as a collective movement of individual cracks. Moreover, experimental results for complex quasi-static loading can be fitted by matching a very limited number of model parameters. This makes us believe that the suggested mechanism of internal friction is primarily responsible for the hysteretic constitutive behavior of micro-cracked solids. The model has not been tested for weak acoustical excitation, since at this level it most likely provokes only partial slip in micro-contacts, which is neglected by the adopted friction model. In addition, small characteristic displacements of the cracks' faces can bring into play the adhesion effects.

Micro-potential model

This theory (see [Ale-05]) represents another attempt to find a physical mechanism for hysteresis in stress-strain. The model is essentially based on two assumptions: i) adhesion implies a non-zero pull-off force for surfaces in contact; ii) in the unstrained state there is a final opening between crack faces. Here friction is not considered. The application of the two hypotheses results in the appearance of a double-well potential for some of the internal cracks in the material. Each of them can be found in one of two states: "adhesive" i.e. when the crack faces are trapped in a potential well corresponding to the non-zero pull-off force, and the open state when the cracks' faces do not touch each other. This mechanism makes it possible to use the Preisach formalism in which a hysteretic element corresponds to a single crack.

Variability in geometric parameters, such as the crack size, and the characteristic roughness of its faces, results in different switching parameters for each hysteretic element. The model predicts theoretical hysteretic stress-strain relationships that qualitatively describe the experimental behavior of micro-cracked solids under weak acoustic excitation.

Adhesion hysteresis model

In this adhesion-based model (see [Sha-94], [Ale-07b]), the stress-strain hysteresis is deduced from hysteresis in adhesion. Indeed, introduction of the surface energy and related pull-off force leads to the formation and rupture of adhesive necks connecting two asperities belonging to opposite crack faces. These necks are formed when two asperities that initially were in contact separate. The presence of a neck means that the asperities interact with some attraction force. When identical asperities, initially distant, are approaching each other, a neck will not form until the real atomic contact is reached, so that at the same distance between the asperities they do not interact. So the presence of the interaction force depends on the movement direction (approaching or separating) that actually corresponds to hysteresis at the micro-level. Again, it is possible to transfer this model to the Preisach description; a hysteretic element is represented in this case by a pair of contacting asperities. In this model, the decomposition of the continuous roughness into a set of individual asperities presents an additional problem.

Again, the model offers the opportunity to qualitatively describe nonlinear acoustic experiments in media with simple geometries.

Soft-ratchet model

This model (see ) is intended to simulate acoustical experiments on the short and long time scale and, in particular, to explain the so-called slow dynamics effects (see [Ost-01]). Slow dynamics is an effect involving various time scales; it consists in a drop in the elastic moduli just after intense straining by a powerful acoustic wave or after an impact.

Then the modulus drop slowly (logarithmically) disappears so that the modulus is gradually approaching but not reaching the original value. The effect is associated in the model with the asymmetry of creation and breaking of adhesive bonds at the micro-level. It is assumed that there is a multitude of ways to break a bound since an individual atom can belong to one or another asperity after breaking. At the same time, there is only one way to reunite two initially separated asperities. This means that the process of breaking the bound is faster than the healing process, and therefore during the intense straining the modulus is most likely to drop.

Pack of Hertz-Mindlin spheres

In this model (see [Nih-00]) for an isotropic material with internal cracks (e.g. sandstones), the real complex material's geometry is replaced by a pack of spheres interacting according to the Hertz-Mindlin laws (see [Min-53]). Known force displacement relationships for elastic spheres in contact with friction were used to calculate the response of a face-centered cubic sphere pack subjected to uniaxial strain loading ε zz . The spheres were considered as unconsolidated, so that, in order to model a solid material, the pack was confined from the sides, ε xx =ε yy =0. This greatly simplified system provides a hysteretic stress-strain response  zz vs ε zz , similar to the experimental hysteretic curves for rocks. There are a number of experimentally observed features correctly predicted by the model, such as a shift of the first hysteresis loop relatively to the subsequent loops for a periodic excitation, the independence of dynamic moduli on strain. Thus, the model has a strong potential in describing microcracked materials, although faces the difficulty in taking into account the material's consolidation.

Recently developed nonlinear acoustic NDT methods

In the last few decades, a strong interest on nondestructive testing (NDT) methods based on nonlinear elastic effects in solids has grown, driven by the request from industry for sensitive quantification and localization of micro-structural damage. Researchers have developed innovative techniques to interrogate the micromechanical behavior of materials and its effect on wave propagation by investigating the amplitude dependence of macroscopically observable properties, as speed of sound, attenuation and so on. As excellent reviews of these developments of nonlinear acoustic NDT techniques are available in the literature [Van-00b],

[Ost-01], [Guy-07], we will concentrate in this part on nonlinear NDT techniques recently developed in the framework of the ALAMSA European project.

Guided wave tomography using RAPID algorithm

Ultrasonic Guided Wave Tomography (GWT) is a non-destructive inspection method capable . RAPID GWT uses a sparse array consisting of N ultrasonic transducers arranged in a way that it covers a designated region of interest where the defect is to be detected, as depicted in Fig. 3(a). Every element in the array can act like a transmitter and a receiver. The measurement procedure is as follows.

First, element i is switched to the transmitting mode and sends out a predefined waveform and the other elements act like receivers recording the propagating guided waves. Then element i is switched back to the receiving mode and element i+1 acts as a transmitter instead. This procedure is repeated either for half of the sparse array elements (reciprocity taken in account) or for the complete array. The whole dataset has to be acquired in two different states. Originally in the conventional RAPID, it was the intact state (baseline) and the damaged state. In the nonlinear version of the method, the whole dataset is acquired using two different amplitudes of excitation. The Scale Subtraction Method (SSM), based on the subtraction of both images obtained using RAPID

(a) (b)
algorithm after scaling, is then used to extract only the nonlinear contribution to the signals. A typical result of the SSM baseline-free nonlinear RAPID GWT imaging is depicted in Fig.

3(b).

Resonant scanning laser vibrometry

In order to considerably increase the sensitivity of nonlinear acoustic imaging methods, a resonant ultrasonic activation of defects via Local Defect Resonance (LDR) can be used. The LDR concept is based on the fact that the inclusion of a defect leads to a local decrease in stiffness for a certain mass of the material in this area, which should manifest in a particular characteristic frequency of the defect. As soon as the driving ultrasonic frequency matches the LDR frequency of the defect, the amplitude of local defect vibrations increases dramatically.

Since LDR is as an efficient resonance "amplifier" of the local vibrations, it contributes appreciably to nonlinearity of the defect which exhibits transition to nonlinear regime even at moderate ultrasonic excitation levels. Under LDR conditions, the nonlinear frequency components (higher harmonics, frequency mixing and subharmonics) dominate in the vibration spectrum of the damaged area.

In the Resonant Scanning Laser Vibrometry (RESLV) technique, the nonlinear source (defect) is visualized via laser detection of nonlinear vibrations: the average vibration spectrum of the specimen is measured with a laser and after the Fourier transformation the nonlinear frequencies are mapped over the specimen surface.

The RESLV experimental setup (Fig. 4) combines piezoelectric excitation of ultrasound at the LDR frequency with laser detection of nonlinear frequency components in the defect area. A direct way to experimentally reveal the LDR frequency is to measure an individual contribution of each point of the specimen in its overall frequency response. For this purpose, an ultrasonic excitation by a wideband piezoelectric transducer is required to probe and indicate all possible resonances in every point of the specimen. The origin of each maximum is then verified by imaging the vibration pattern in the specimen at the corresponding frequency. A strong enhancement of the vibration amplitude observed locally in the defect area is identified as a fundamental defect resonance. Such procedure was implemented in all resonant methods developed prior to nonlinear imaging.

Resonant thermosonics

The use of LDR, which strongly intensifies local vibrations and results in an efficient higher harmonics generation, is beneficial for enhancing the defect thermal response and therefore the efficiency of ultrasonic thermography in Resonant Thermosonics (RETS) technique.

The RETS nonlinear imaging technique combines the resonant ultrasonic excitation with infrared imaging. Prior to RETS measurements, a scanning laser vibrometer is used to identify the LDR frequency of the defect in question.

Scan area

Defect(crack) Scan tool

Laser interferometer Excitation

Resonant shearosonics

A combined effect of the LDR and the defect-induced nonlinearity also increases the sensitivity and contrast of interferometric optical imaging such as Resonant Shearosonics (RESH). Shearography is a speckle-interferometrical method that uses images of speckle patterns at the sample surface. The pattern correlates with the shape of the object surface and can therefore be used to determine the object deformation within a fraction of the laser wavelength. The deformation produced by an external excitation (pressure, heating, ultrasound) is particularly strong in the positions of near-surface defects.

In the RESH technique, the specimen is excited at the LDR frequency in a continuous ultrasonic generation mode by using piezo-transducers. An optical sensor with continuouswave laser and CMOS digital camera is used for RESH imaging (see Fig. 8). The speckle pattern of a vibration object is integrated over the frame period and compared with a reference speckle. Due to the LDR-induced local resonant increase in defect vibrations and its nonlinearity the RESH sensitivity and contrast of imaging is substantially enhanced.

Coda wave interferometry

Coda waves are the tails of signals received in samples where the waves repeatedly sample the same region in space due to multiple scattering produced by inhomogeneities in the material (such as holes, inclusions, boundaries between different materials and so on) or due to the reverberation in samples of finite size (see Fig. 1).

Coda Wave Interferometry (CWI) is a technique for monitoring time-lapse changes based on the phase and amplitude information of coda waves. In this technique the medium works as a natural interferometer. Thus, even extremely weak or local changes in material properties or sample geometry produce measurable deviations of the coda signals. With the source and receiver unchanged during the experiment, the CWI analysis can detect perturbations of the medium by comparing coda waves. The system resembles a billiard where a small change in the initial conditions results in a considerable final change in the configuration due to multiple interactions between the ball and the billiard and the effect of an accumulation of errors.

Doublet technique and Snieder's model

Changes in the waveforms, as shown in Fig. 2, can be quantified by computing the timewindowed cross-correlation coefficient defined by: , 22

( ) ( ) () ( ) ( ) c c c cc ip cc tT i p s tT t T s t T t T t T t T u t u t t dt Rt u t dt u t dt         , ( 1
)
where t s is the time lag for the correlation, u p is the perturbed wave field, u i is the reference wave field, and t c is the center of a time window of length 2T. Accepting some assumptions about wave propagation conditions, it is possible to analytically obtain the time delay δt at which a maximum of correlation , c tT R is attained for different types of distortions.

Following Snieder (see [Sni-06]), we express the wave field displacement at a given location as the sum of the displacements of waves propagating along all possible trajectories t r : ( ) ( )

i tr tr tr u t A S t t   , ( 2
)
where t denotes time, S(t) is the source signal, and A tr and t tr are the amplitude and the travel time, respectively, of the wave following trajectory t r . Equation (2) does not specify how exactly the waves are scattered and just corresponds to the superposition principle. The trajectories include the direct propagation as well as single and multiple scattering including those produced by the nonlinear scatterers, whereas the propagation process itself is assumed to be linear. Accepting that typical wavelengths In what follows we assume that the materials quality factor is high enough so that Eq. ( 3) is valid.

Substituting Eq. ( 2) and (3) into the expression for the time-windowed cross-correlation coefficient given by Eq. ( 1) and using Ladder's approximation [ (the cross terms with different paths (tr  tr') are incoherent and thus negligible), which is valid for 21 fT    where f  is the width of the coda wave spectrum, one arrives at the following expression: For small time shifts in comparison to the dominant period of the wave, a second order expansion can be made:

(6)
where 2  is the mean-squared frequency of the coda waves that arrive in the time window:

(7) Thus, we obtain: where t c is the center time of the employed time window. Thus, the travel time perturbation depends on the arrival time of the coda wave, but not on the particular path followed, and

  2 2 , ( , ) 1 ( ) 1 2 c c t T s tr s tT R t t     , (8)
                          . ( 12
( , ) () c tr c tT tt       . ( 15
)
In a short time window we have thus 2 0

   and () c t t     . ( 16
)
The use of the normalized correlation coefficient given in Eq. ( 1

  

, this method is insufficiently precise since various paths become undistinguishable. In that situation, the use of the stretching method is preferable.

Random displacement of scatterers

Considering scatterers that move independently in three dimensions with a root-mean-square displacement  between two measurements, Snieder given by:

2 2 * 2 L t l    , ( 17
)
where l * is the transport mean free path, linked to l e by [Sch-97]:

* 1 1 1 ea l l l , (18)
where l a is the absorption length. When the velocity is constant then the variance in the travel time is related to the variance in the path length by:

L      . ( 19
)
Thus,

2 2 * 2t l      .
In this case, the maximum of the cross correlation provides pertinent information, which is linked to the mean displacement of the scatterers by:

22 , * max ( ) 1 c s t T s t Rt l    , ( 20
)

Doublet technique

From Eq. ( 1) two parameters can be extracted. The first one is t the time shift that corresponds to the maximum of the correlation coefficient. This parameter is sensitive to velocity perturbations as shown in Eq. ( 15). The second one is the decorrelation coefficient which indicates the level of decorrelation (from 0% to 100%) between two coda waves. The second parameter is sensitive to geometrical or structural changes of the scattering medium, e.g. random displacement of scatterers (see Eq. ( 20)).

In the doublet technique, the entire coda is decomposed in non-overlapping windows of width T and center t i c . For each window the time shift that corresponds to the maximum of the 55 correlation coefficient t i is calculated. Then, from the slope of the curve t i versus t i c we obtain the relative velocity perturbation   /   (see Eq. ( 16)).

Stretching technique

The stretching method consists in calculating the correlation coefficient between the perturbed coda and the reference coda interpolated at times t (1+α), where /

  

( (1 )) ( ) c c c cc ip cc tT ip tT S tT t T t T t T t T u t u t dt R u t dt u t dt             . ( (1 )) ( ) ()
) 21
The stretching method is especially suitable when the perturbed coda approximately behaves as a time stretched/compressed version of a reference coda. Indeed, this technique directly includes the fact that / const    , and thus that the induced phase shift changes linearly with time:

 

1 t t t t         . ( 22
)
Thus, the relative variation of the velocity is determined by finding the value  max that maximizes , ()

c S tT R  : max     .
In this technique the entire coda wave is treated at once. Thus, we need to interpolate the reference coda wave u i at times t(1+ i ) for different values of 

      , max max 33 2 , max 1 ( ) 62 2 ( ) (2 ) c c S tT S tT c c c R f RMS R f t T t T             , ( 23
)
where f c and f are the center frequency and the bandwidth of the coda wave spectrum, respectively. RMS  is the apparent fluctuation of the stretching parameter induced by noise. In this expression, 1 f  is the amount of time for one bit of information to be delivered, and corresponds roughly to the time of the initial source. The quantity

      T t t T    
represents the available time where coda waves are significantly larger than the noise. This RMS  value can be measured using , () value is higher than the RMS  value, then the measurement can be considered as meaningful.

It is important to note that the stretching method requires stretching of the measured coda waves from t = 0 up to t = α(t c +T). Thus, it is essential to precisely synchronize all the acquisitions with the source signal emission.

As for the doublet technique, two parameters can be extracted from Eq. ( 21) in the stretching technique. The first one is the stretching factor ε that corresponds to the maximum of the correlation coefficient for all values of α. The second one is the decorrelation coefficient

, 1 ( ) c S tT Kd R  
that indicates the level of decorrelation (from 0% to 100%) between two coda waves.

Choice of the time window used in the CWI technique

The precision of CWI techniques is highly dependent on the choice of the time window used for the analysis. Coda waves are long duration signals that probe cumulative effects. Thus, various parts of the coda wave contain different levels of information about the defect. The first part of the coda is just the ballistic signal coming directly from the source, and thus is to be disregarded as it contains very few information on the defects. Aspects related to the Signal to Noise Ratio (SNR) in the last part of coda trail have to be carefully treated in order not to interpret the noise as a defect.

We will describe now certain rules that help to choose correctly the time window used for a proper CWI analysis.

Doublet technique

1. The central time of the window t c needs to be higher than the mean free path l * :

* 4 c l t   is
generally sufficient.

2. The length of the window T needs to be sufficiently large in order to verify the coda wave interferometry criterion:

21 fT    . Thus 1 2 T f   .

t c needs to verify

c tT  to be sure that t max =τ tr is nearly constant in the window T. 4. The signal to noise ratio needs to be sufficiently high in the considered time window T.

This criterion generally induces a higher limit in time in the choice of T .

Payan et al. [Pay-09]

have also indicated that, in order to obtain a good evolution of the correlation coefficient, T needs to be at least of the order of 10 periods of the lowest frequency component of the coda wave.

Stretching technique

In the stretching method the chosen window T is generally quite long and the criterions 2

( 1 2 T f  
) and 5 ( 10 Thus, the choice of t 1 is linked to the reverberations of the sample in the considered frequency range. The choice of t 2 can be increased if the quality of the measurement (obtained SNR) is improved. One possibility is to make an averaging over a high number of measurements.

59

Determination of the signal to noise ratio

To determine correctly the time window used in CWI techniques we need to automatically calculate the SNR for each measurement. The method used is described schematically in Fig. 5. The coda wave is divided in parts of length t. Each part is Fourier transformed using a

Hanning window to reduce Gibbs oscillations in the obtained spectrum. We then calculate the intensity of the incoherent part of the coda wave ˆ( , ) . Finally, the signal to noise ratio is obtained by dividing ˆ( , ) cc I t f by the value of Î obtained before the beginning of the coda. In all our experiments, t 2 will be the time at which the SNR becomes lower than 50dB in order to ensure a high level of precision for the determination of the stretching factor ε and the decorrelation coefficient Kd (see Eq. ( 23)). underwater acoustics, due to its ability to provide spatial and temporal focusing of an ultrasonic wave. Time-reversal invariance in acoustics means that for every burst of sound s(r,t) emitted from a source and then reflected, refracted, or scattered by heterogeneities of the propagation medium, there exists a set of waves s(r,-t) that precisely retrace all these complex paths and converge at the original source, as if time was going backwards. This invariance is satisfied by the equation in non-attenuating media. The TR process leads to a spatial focusing and a temporal compression.

Spatial focusing means that the time-reversed field focuses back exactly at the source.

Temporal compression means that the time-reversed signal at the source is similar to the signal previously emitted by the source. In other words, the result of a TR process is that waves recorded on the boundary are focused back in space and time on the acoustic source, or on the scattering targets inside the region that were acting as sources. For the classical linear TR process, the returned signal focuses on the direct wave source position and not on the

Principle of one-channel time reversal

In

Chaotic cavity transducer

In order to improve the quality of the TR process, we use a so-called chaotic cavity transducer, consisting of a combination of a PZT ceramic glued to a cavity of chaotic shape on the hardware side with the time reversal principle on the software side. An applied source signal to the PZT ceramic generates a wave propagating in the cavity and the sample medium. Each time when the propagating wave in the cavity arrives at the boundary between the cavity and the sample, some part of the incident energy is reflected and continues to engender multiple reflections on the other boundaries of the cavity, whereas the other part of the energy is transmitted into the sample, as shown in Fig. 7. One of the advantages of chaotic cavity transducers is their ease of manufacture and their low cost. Two examples of manufactured chaotic cavity transducers are displayed in Fig. 8.

LabVIEW data analysis tools for one-channel time reversal focusing

A program using LabVIEW software, commercialized by National Instruments, has been developed for driving one-channel time reversal focusing experiment (see Fig. 9) and data analysis. The choice of LabVIEW, that is a graphical programming platform, was made because this software is ideal for any measurement or control system in a wide range of applications in a relatively short amount of time.

Sweep compression

Experimentally, due to the low energy in the pulse, the use of a sinusoidal pulse signal for impulse response measurements does not provide a strong received signal, leading to a poor signal to noise ratio.

  2 0 ( ) cos(2 2 ), 22 e B T T s t A t f t t t T       , (27)
whereas a rectangular-shaped window leads to Fresnel ripples at the frequency band edges.

1-bit time reversal

A 1-bit process can be used in order to increase the amplitude of the retro-focalized wave

[Mon-01].
During classical time reversal or inverse filter experiments, both the instantaneous phase and amplitude information of the received signal are sent back. For a 1-bit method, the time reversed or inverse filtered signal will be set to 1  depending on the sign of signal: 1  amplitude is set if the sign is positive, otherwise, 1

 amplitude is set. This means that the instantaneous phase information present in the zero crossing of the signal is time reversed while the instantaneous amplitude information is ignored.

Inverse filtering (Zverev)

The spatio-temporal inverse filter approach used by Tanter et al. The classical inverse filtering method is based on the inversion of the propagation operator H.

In practice this inversion is performed in the frequency domain for every frequency component of the received signal spectrum R(ω) [Zve-04]:

    1 () E H R      . ( 28
)
The inverse filtration can be applied only to those parts of the signal spectrum at which the modulus of the function H(ω) differs from zero. The parts at which this modulus is equal to zero should be excluded. To avoid errors caused by these zeros, it is sufficient to set to zero H -1 (ω) if the modulus of H(ω) takes a value below some threshold level. In practice a value of one hundredth of the maximal value of H(ω) in the considered frequency band for the threshold is a good compromise.

1-bit inverse filtering

The 1-bit inverse filtering process is similar to the one described for 1-bit time reversal. Here, the 1-bit process is performed after application of the inverse filtering described in the previous subsection, instead of use of the time reversal signal processing.

Inverse filtering (Quieffin)

In the one-channel TR system, the inverse filtering approach consists in the inversion of the eigenmode energy [Qui-04]. The procedure performs an inversion of the energy of the eigenmodes and constructs the re-emitted signal as a linear combination of all the eigenmodes of the cavity, weighted by this inversion. Doing so, the focusing process takes advantage of all the modes including those with the weakest energy, which are poorly exploited in the time reversal focusing process. Using this concept of eigenmodes, the inverse filter transfer function can be calculated as:

      1 2 H H H       , ( 29
)
where H + (ω) is the Hermitian transform (transposition and conjugation) of the propagation operator H, and

 

H  is the mean value energy of the eigenmode at frequency ω.

3.4.6. 1-bit inverse filtering (Quieffin) This signal processing method is similar to the 1-bit inverse filtering but with the use of the

Principle of nonlinear coda wave interferometry

In the nonlinear version of the CWI technique, small variations of the time delay ε and the envelope Kd of the coda wave are induced by nonlinear interactions between the coda wave and a low frequency pump wave with increasing amplitude. These variations of ε and Kd are indicators of the defect's presence. Thus, in this technique, in addition to the HF coda wave, a LF frequency pump signal should be generated. The LF component of the measured coda should be then filtered out. Therefore, in the absence of defects the filtered HF coda is not affected by the LF pump signal and remains unchanged. This is not so when nonlinear interactions are induced by the presence of defects. In this case, due to contact acoustical nonlinearity, the LF signal will generate some modulation of the HF coda via the frequency mixing process. This modulation, i.e. the change in the HF response, is measured by the CWI technique. This means that we have to measure the optimal time stretching factor ε and the decorrelation correlation coefficient (Kd) for a pair of coda waves produced without and with a LF (10-100 kHz typically) high amplitude pump. We use the stretching technique to determine these small variations. The measured variations increase with the quantity and density of the nonlinear scatterers, i.e. the defects. In materials without any nonlinear scatterers, i.e. in intact samples, no nonlinear interactions are observed. The LF pump wave is emitted in the sample using a chaotic cavity transducer with the use of the time reversal concept. Time reversal gives a means to concentrate the pump wave in a given region of the sample, thus providing an opportunity to perform a localized CWI analysis, as shown in Fig.

12. The possibility to delay the CWI probing relatively to the pumping is also introduced to enable the study of slow dynamics effects.

Experimental setup

Experimental setup description

The experimental setup illustrated in Fig. 13 contains two components (CWI and pumping).

Pumping component

To emit the pump wave, an AWG710 arbitrary waveform generator connected to a rüel Kjaer 2713 Type Amplifier is used. When a focusing of the pump wave is needed, i.e. when we want to localize the defect combining the nonlinear CWI technique and the TR process, the acquisition of the signal for the TR procedure is performed by a LeCroy 64 Xi

Oscilloscope. The pump signal can be emitted continuously or with a finite duration T PUMP . In this case, it is emitted synchronously with the probe signal after each T TRIG time interval, including a possibility to delay the probe relative to the pump by a time T DELAY (see Fig. 14).

The pump wave is emitted by a piezodisk glued to a metallic chaotic cavity in a form of truncated 40 mm diameter cylinder. This chaotic cavity transducer is a part of the TR procedure used for damage localization. To acquire the signal at the focal point needed for the time reversal process, one can use a laser vibrometer, an ultrasonic pinducer (needle transducer), or a small diameter piezoelectric disc. We used 10.5 mm diameter radial polarized piezoelectric disks KEPO FT-10.5 due to their simplicity and low cost. They have been glued to a glass plate at the location of the Time (ms) defect and in three arbitrary locations schematically shown in Fig. 13 by point's p1, p2, p3 and p4. The resonance frequency of these piezoelectric disks is 8 kHz. As we operate in a frequency range far away from this resonance frequency, their transfer function is nearly constant in the used frequency range. The range of pump frequency is 50-300 kHz. We use frequency swept signals in order to cover several resonance frequencies of the sample and to avoid the existence of "dark" regions due to resonance patterns. The pump can be injected continuously (and thus asynchronously) regarding to the probe signal or it can have finite duration T pump and be emitted synchronously with the probe signal after each T trig time interval.

In the latter case there is a possibility to delay the probe relative to the pump by the time T delay .

These synchronization parameters are managed by a pulse generator (we use a marker channel of a AWG710B generator). The pump wave is emitted by a piezoelectric disk of 40 mm diameter and 15 mm thickness glued to a metallic chaotic cavity in a form of truncated 40 mm diameter cylinder. A detailed description of the focusing procedure with time reversal and chaotic cavity transducers is given in section 3 of this chapter.

Choice of the signal processing method used in chaotic cavity transducers

For the experiments reported below, a glass plate with a localized defect was used, as shown in Fig. 15 Time reversal focusing procedure may encounter some problems in a medium of regular geometric shape (non-ergodic cavity). Indeed, in such a sample, regular resonance patterns can appear that worsen energy focusing. This is not the case when the chaotic cavity transducer is used. In Table 1 the ratio of RMS values and amplitudes between the current focusing point and the average value in two other points (on the plate) are presented. We used two options to concentrate the ultrasound field in a given point: the Quieffin inverse filter due to its ability to obtain good shape of the recompressed signal and good contrast in peak amplitude, and the time reversal technique with normalization of the envelope (similar to 1-bit time reversal) due to its ability to create a focused field of high power. We can see in Table 1 that the measurement produced at the defect provides a quite strong contrast in RMS and amplitude, for both the Quieffin filter and the time reversal procedure.

The RMS contrast for the Quieffin filter is even stronger since at p 1 and p 2 the measured RMS is negative. For the time reversal procedure, all indicators are positive in all cases but at p 2 the RMS is weak. The situation remains the same in both frequency ranges of 35 -150 kHz and of 120 -220 kHz.

In the TR process, one important parameter is the length of the time window to be processed (time reversed or inverse filtered). We measured the evolution of the energy contrast (ratio between the energy at the focal point and the one at 4 cm away from the focal point) as a function of this delay (length of the processed time window) used in the TR process. A 900 μs sweep in the frequency range 50-250 kHz is used to excite the chaotic cavity transducers. The results obtained for two signal processing methods, Quieffin inverse filter (red line) and TR (black line), are displayed in Fig. 16. The maximum energy contrast of 2.2 measured with Quieffin inverse filter is higher than the one measured when using time reversal (1.9). In both cases, the maximal value of the energy contrast is obtained for a length of the processed time window longer than 4ms. Certainly, this length depends on the length of the received coda wave and thus on the reverberating properties of the tested sample. Nevertheless, this study shows that the minimum length of the processed time window does not depend on the used signal processing technique for time reversing the received signal.

These results clearly indicate that the Quieffin inverse filter gives better results than a simple TR technique. Thus, this inverse filter is used in all the further experiments.

Choice of the frequency range of the sweep used in chaotic cavity transducers

We measure the focal spot on a glass plate with a VP-0.5 pinducer for three different frequency ranges: 15-50 kHz, 50-125 kHz, and 150-250 kHz. The results are displayed in Fig. 17. As expected, the size of the focusing spot decreases (from 3 cm in diameter to 1 cm) with increasing frequency. Nevertheless, as the pumping frequency has to be decoupled from the one used for the probing coda wave (see Fig. 15(c)) it cannot be increased above 200-250 kHz.

Indeed, in order to obtain a sufficiently long coda wave, the frequency range of the pumping limited by attenuation in the sample is typically lower than 1 MHz. It is important to note that this limits the resolution of the method to about 1cm.

Comparison with other source types

The optimized chaotic cavity transducer is compared with other acoustical sources (three PZT In conclusion, our study has demonstrated that the optimized chaotic cavity transducer can be effectively used for pump energy focusing in the 50-100 kHz frequency range, as required for the application of the time reversal coda wave mixing technique.

Measurements and results

Measurements on a thermally shocked glass plate

The nonlinear coda wave mixing technique was first applied to a thermally shocked glass plate of irregular form (approximately 20 x 15 cm) shown in Fig. 19. It contains numerous uniformly distributed cracks at the surface and in depth. For global inspection of the sample no focusing is needed so that the use of the chaotic cavity for pumping is not necessary. In this case as a pump emitter, we used a rectangular piezoceramic plate of dimensions 45 x 15 x 1.5 mm 3 with a resonance frequency of approximately 77 kHz directly glued on the lateral side of the sample (Fig. 19). Similarly to the frequency shift observed in nonlinear resonance experiments [ which is often considered as an indicator of micro-damage in materials, we expect a velocity variation dependent on the acoustic excitation amplitude which can be correlated with the elastic nonlinearities of the sample. In order to detect the amplitude dependency, the following protocol of amplitude A pump was applied (see Fig. 20(a)):

1. First, a reference wave with no pump is recorded. Each coda wave recording is performed with 64 sweeps. The procedure takes approximately 1 second. This record is used as reference for all further measurements.

2. A second measurement is performed without pump and the stretching analysis of ε and Kd is applied. The result of the measurement is used in order to estimate the drift of the coda wave and the precision of measurement.

3. During the next 10 steps, the pump wave amplitude A pump is increased step-by-step from 0 to maximum (150 V or 160 V depending on the case), with a constant increment at each step a CWI measurement is made and analyzed.

4. Next, the amplitude is raised to its maximum level and held at that same level for the five subsequent steps.

5. Afterwards, A pump is decreased directly from maximum back to 0 and then remains at this level for the five final steps. This procedure is executed in order to estimate a possible slow dynamics effect.

6. Each sequence of measurement (1 -5) is performed 10 times, then the mean value and error are calculated.

Fig. 20. (a) Measurement protocol used during the test of the highly cracked glass sample with nonlinear CWI. Results of CWI analysis in continuous mode: evolutions of (b) stretching parameter ε, and (c) Kd. A comparison is given for an intact glass sample.

As illustrated in Fig. 20, for the damaged specimen, when A pump voltage is increased from 0 to 160 V, the propagation velocity of the probing coda waves decreases by 0.41% (Fig. 20(b)), and the waveform shape is modified (Kd exhibits a 3.2% increase in Fig. 20(c)). When A pump is decreased back to 0 after an excitation to 160V, ε and Kd return to values very close to the initial ones obtained at the beginning of the test (i.e. with no pump). We can observe a slight drift of the ε parameter, possibly due to some amount of heat generated by the pump piezoceramic plate. In addition, a very small and slow dynamic relaxation effect after conditioning at 160V is observed. Such relaxation effects, that occur at time scales much larger than the pump acoustic wave period, represent a contribution to the nonlinear slow dynamic behavior Experiments were carried out using a constant amplitude pump with a finite duration and with a variable probe delay (see Fig. 14). A mono-frequency pump of 76.88 kHz (resonance frequency of the PZT ceramic used for an excitation) was used with a fixed duration of 800 µs.

The time delay of the probe was varied from 0 to 2ms. Pump repetition rate T trig is 10ms and pump amplitude is 150 V. As shown for the evolution of ε and Kd as a function of T delay in Fig. 21, when period of the probe wave become higher then period of the pump wave the decay of the curves is not instantaneous. The nature of this residual effect could be linked to slow dynamics, a well-known effect in nonlinear acoustics. We have also studied the influence of the pump duration on these slow dynamic effects appearing in nonlinear coda wave mixing. To do this, we measured for different duration of pumping, i.e. for different values of T pump (from half period of excitation (6.5 µs) to 800 µs), the evolution of the measured stretching parameters ε and decorrelation coefficient Kd. The results obtained for Kd are shown in Fig. 22.

Measurements on a rectangular glass plate with impact damage

The nonlinear coda wave mixing technique was then applied to a rectangular glass plate of 30

x 50 x 1 cm 3 . This plate contains a crack due to a bullet hit from a pneumatic gun that mimics a realistic local defect (see Fig. 25). An intact glass plate of the same size, i.e. 30 x 50 x 1 cm 3 , was used as a reference object.

Probe emitter Pump emitter

Probe receiver Fig. 25. Rectangular glass plate of 30 x 50 x 1 cm 3 dimensions containing a crack due to a pneumatic gun bullet hit.

On this sample, the CWI measurements were performed in two modes:

1.

Continuous asynchronous mode.

Results obtained with this mode are denoted as "continuous". Pump duration is T pump = 2ms and repetition rate for probe signal is T probe = 3.3ms.

2.

Pulsed pump with delayed probe. This mode is denoted as "pulsed". In this mode, duration of pump is T pump = 800µs and repetition rate for probe signal is T probe = 10ms.

The delay T delay between the pump and probing waves can be varied. In the pulsed mode, these values return to the initial value while in the continuous regime, after switching off the pump, the absolute values of ε and Kd decrease but don't return to the initial values. This is presumably due to the heating of the sample induced by the continuous excitation of the pump transducer. As for the highly cracked sample, no changes were measured for the intact sample.

Results of the CWI analysis in continuous mode

The developed system enables monitoring of ε and Kd changes in real time, as shown in Fig. 27, for the intact and locally damaged glass plates. Here, the stretching analysis is performed in continuous asynchronous mode, where each measurement takes about 0.1s. Here, pump varies in binary manner between 0 and max value. The pump is switched on and off manually at different times and for different durations in order to characterize particular features of nonlinear coda mixing under different conditions.

The overall variation of the stretching parameters is much higher for the damaged sample than for the intact sample: 1% vs 0.2% for ε and 3% vs 0.4% for Kd. Periodical fluctuations of CWI parameters seen in Fig. 27 for the intact sample are typical for temperature fluctuations.

In the damaged sample, for short pump durations one can see a quick growth and drop of the curves at the moments of pump switching on and off. For the long pump duration, after the first part of quick growing one sees a slow drift of the curves superimposed with thermal fluctuations. After switching off the pump the CWI parameters slowly relax to the initial values.

- In order to localize the defect in the glass plate with a local crack (Fig. 25), the CWI technique was combined with the TR process, as described in Fig. 28(a). The CWI measurements were synchronized with the TR process in order to perform the stretching analysis in a narrow time window applied to the part of the coda signal that is sensitive to the ultrasound pumping. Fig. 28 illustrates cases when the pump wave is not focused (green), SSM is not a spectroscopy method, since it does not consider information in a particular range of frequencies in the way band pass filters or high/low pass filters do. Instead, using SSM, the nonlinear information is obtained as the result of an operation in the time domain. We chose to use the SSM analysis as a post-processing tool since it can be used with a broad band excitation and in not only limited to continuous wave and long burst excitation signals as harmonic and intermodulation filtering methods. Moreover, one of the primary advantages of SSM is that it yields a high signal-to-noise ratio. This is due to the fact that the more traditional filtering methods extract nonlinear signatures only at the sub-or higher harmonics content, while SSM manages to extract additional nonlinear signatures that also affect the fundamental frequency range.

Local defect resonance excitation

The efficiency of the NACU technique can be enhanced by choosing the excitation frequency equal to one of resonance frequencies of the sample. However, defects introduce certain changes into the structure of resonant frequencies in a sample. In particular, new frequencies can appear which were absent in the intact sample. This phenomenon called Local Defect

Resonance can be used to drastically enhance the sensitivity of the method or to lower input power, as first described by Solodov in [Sol-11], [Sol-13]. In fact, the acoustical energy distribution in a sample excited exactly at the LDR frequency has a pronounced maximum in the vicinity of the defects. In other words, the method provides an efficient energy pumping of the defect and not to the whole sample. As the excitation amplitude increases, the LDR-"amplifier" (up to 20-40 dB) exhibits a transition to the nonlinear regime with higher harmonic generation only in the defects area.

To find an LDR frequency, a certain procedure has to be followed. It consists of frequency scans and visualization of the standing wave patterns at the sample's surface. When the excitation frequency matches the LDR frequency, the vibration amplitude of the defect significantly increases. A strong local maximum in the standing wave pattern indicates the presence of a defect in the vicinity of this position. For further analysis, the sample should be excited exactly at this frequency since it allows one to selectively "insonify" the damage.

It is important to note that the choice of the SSM method for nonlinear components extraction is preferable in comparison to more classical filtering techniques, as the LDR frequency is not known a priori for a given sample.

Experimental setup

The setup developed for damage detection using the NACU technique with SSM postprocessing is shown in The list of needed equipment is as follows:

 two arbitrary waveform generators capable of sending the following waveforms: sine burst, chirp/sweep, pulse. The frequency range is 10 kHz up to 10 MHz, with maximal amplitude of a few volts (depending on the amplifier used).

 Amplifier should be able to operate in the frequency range 10 kHz -10 MHz with a maximum output up to 200 V to a PZT ceramic (a KEPO piezodisc for example).

 A pair of air coupled ultrasonic transducers (focused if possible in order to improve the imaging resolution) suitable for a given frequency range.

 A reception amplifier dedicated to air coupled ultrasonic transducers.

 An automatic motion stage.

 An oscilloscope connected to a computer (via GPIB for example) to acquire the data.

Measurements and results

Below we describe the results obtained via the NACU technique combined with the SSM post-processing on samples fabricated by partners of the European ALAMSA project.

Measurements on a CFRP laminate with a 35x35 mm 2 delamination at half thickness

The first sample is a 300 x 300 mm 2 pre-preg (pre-impregnated composite fibers) CFRP panel with a thickness of 3 mm and containing a 35x35 mm 2 artificial delamination. The CFRP is made of 12 plies with alternating orientations (0°/45°). The delamination was introduced at half thickness into the sample during the manufacturing process by the insertion of a small sheet of Teflon.

The plate was vibrated with a KEPO PZT ceramic glued on its bottom left corner. The PZT ceramic was excited with a sweep from 20 to 220 kHz at two different amplitudes, 20 and 60

V. A scan of the response of the plate was measured with the wide-band Microacoustics air coupled ultrasonic transducer placed 5 mm from the CFRP panel.

Fig. 2. (a) Mean spectrum of the plate vibrations at the defect (A in , in red) and outside of it (A out , in blue). (b) Ratio A in /A out .

Using the localization of the delamination made with a through transmission air-coupled ultrasonic image (see Fig. 4a) we calculated the mean spectrum of the vibration of the plate on (Fig. 2a in red) and outside (Fig. 2a in blue) the defect. The ratio A in /A out of these mean spectra shown in Fig. 2b indicates the frequency band in which the response at the defect position is higher than the one measured elsewhere. This ratio is chosen as a parameter for the NACU imaging at frequencies between 120 to 220 kHz. The image obtained by subtracting three times the results obtained at 20 V from the one measured at 60 V is displayed in Fig. 3 to the through transmission NACU image, we can see (Fig. 4) that, despite that the boundaries of the defect are less pronounced, the contrast is better. Now, if we use only a limited frequency range, corresponding for example to only one of the peak of the spectra displayed in Fig. 2b, then only a small part of the defect is highlighted. This is shown in Fig. 5 for a selected frequency band ranging from 148 to 152 kHz.

Measurements on a CFRP laminate with a 20x20 mm 2 delamination at 1/4th thickness

The second sample is a 300 x 300 mm 2 prepreg CFRP panel with a thickness of 3 mm and containing a 20x20 mm 2 artificial delamination. The CFRP is made of 12 plies with alternating orientations (0°/45°). The delamination was introduced at 1/4th thickness into the sample during the manufacturing process by the insertion of a small sheet of Teflon.

As for the previous composite sample, the CFRP laminate plate with delamination of 20x20 mm 2 was vibrated with a KEPO PZT ceramic glued on its bottom left corner. The PZT ceramic was excited with a sweep from 20 to 220 kHz at two different amplitudes, 20 and 60

V. The scale subtraction NACU image, obtained by subtracting three times the results obtained at 20 V to those measured at 60 V, is shown in Fig. 6. In the first measurement set, the air-coupled ultrasonic transducer scans the plate from the surface closest to the in-depth delamination position. Again, the image presents a high contrast allowing a clear determination of the defect position. The contrast is higher than the one obtained when we use only the data measured with excitation amplitude 60 V, shown in Fig. 7. Now, making exactly the same measurements but from the other side of the plate, thus farthest from the in-depth delamination position, no increase in vibration amplitude appears, as shown in Fig. 8, even if we scan the whole frequency range. This confirms the fact that the detected LDR frequency corresponds to vibrations of the thin layer of plate between the delamination and the surface of the sample. Thus, when the delamination is deeper inside the sample, it hardly induces a detectable local resonance. Nevertheless, in our experiments, it appears that, when we select the frequency band 100-200 kHz, the delamination can be localized as a part of the plate with a decrease of amplitude of vibration, as shown in Fig. 9.

Measurements on a GRFP sample

The third sample is a 4 mm thickness GFRP plate containing an artificial delamination with a 6 mm radius (see Fig. 12(a)). The GFRP is manufactured by an ALAMSA project partner from plain weave (0°/90°) plies. The delamination was introduced into the sample during the manufacturing process. It's a realistic delamination without any kind of inclusion such as Teflon or Kapton.

We first made a linear through transmission air-coupled ultrasonic image of this GRFP sample. We excited the focused airborne transducer with a sweep signal in the frequency range 300-700 kHz, and received the transmitted signal with a wideband non-focused receiver.

As shown in Fig. 10, contrary to the Teflon inserts in the CFRP panels, the realistic delamination is not visible on the obtained image with this linear ultrasonic technique. Though less visible than in Fig. 11, the delamination is the place where spots with maximum amplitude appear.

To improve the detection of the defect, we also subtracted the result obtained at 60 V from a through transmission measurement made with a sweep in the range 300-700 kHz without additional low frequency vibrations. In this experiment we observe nonlinear mixing between the low frequency pump and the airborne acoustic wave crossing the sample. The obtained image is show in Fig. 14; we see that nonlinear processes between by the low frequency pumping (at the LDR frequency) and the airborne acoustic wave crossing the sample take place mainly at the delamination position.

Conclusions

The scale subtraction NACU imaging technique has been developed and used for visualizing defects in composite plates. It was demonstrated that the sensitivity of the technique could be drastically improved by using a vibration frequency range that includes several local defect resonances. Indeed, if we use only a limited frequency range, corresponding for example to only one LDR, then only a small part of the defect is highlighted.

In the future, we plan to combine this scale subtraction NACU technique with the classical through transmission air-coupled ultrasonic imaging techniques. In this case, the sample is placed between a pair of air-coupled ultrasonic transducers. These transducers can be both focused (optimal for the spatial resolution), or only the emitter is focused and a plane receiver is placed at the focus point of the emitter. Vibrating simultaneously the sample with a PZT ceramic glued on it, we can measure the changes induced on the transmission of the aircoupled ultrasonic waves (by nonlinear attenuation or/and nonlinear interaction). As in the previous described methods the sensitivity can be drastically improved by using a vibration frequency equal to the LDR frequency.

Brief history

As we have seen, contact models play a key role for a successful properties description of materials with internal contacts. Theoretical results essentially depend on factors taken into account or neglected by a particular contact model underlying the description at the material (macroscopic) level. This chapter is devoted to the fundamental mechanical problem of frictional contact between two solids having certain surface profiles.

The Our objective is to essentially modernize the classical results by developing more universal approaches applicable for a wide range of excitation protocols and contact geometries.

Geometric extensions

The geometric aspect of the problem can be successfully dealt with by using the reduced defined in the analogous way. For constant loading i.e. when the system is loaded first by a constant normal force and then by a constant tangential force, the notations {N, T } and {a, b} for the vector force and displacement are accepted. The general problem consists in establishing a link between {N, T x , T y } and {a, b x , b y } for both force-driven and displacementdriven systems. For 2D loading the force and displacement vectors stay in one plane so that the notations {N, T } and {a, b} are more suitable.

The objective of this section is to understand how to take into account various contact geometries. This can be made by substituting a given contact geometry with a pair of axisymmetric profiles with the same normal reaction law. Such a possibility is guaranteed by the REFP at least as a good approximation.

The REFP is an important theorem of contact mechanics which states that, for constant loading and for a wide range of contact geometries, the tangential force and displacement can be expressed through the normal force and displacement. This principle is illustrated in Fig. 1 (b) for axisymmetric bodies. Consider two situations: first one where the system is loaded only by normal force Q, and second one where both force components, N and T, are applied (N>Q). The force values are chosen in such a way that the stick zone in the second case coincides with the contact zone in the first case. Then the tangential force and displacement in the second situation are given by

        NQ T N Q b a N a N            , ( 1
)
where the dependency of the normal displacement on the normal force a=a(N) is considered as known, and  is a material constant that depends only on Poisson's ratio ,  

2 21       . ( 2
)
An important feature of Eq. ( 1) is that it does not contain any geometry-related characteristics.

Thus a simple consequence of the REFP is a statement that, for two contact systems with the same normal response, the tangential responses are also identical (see Fig. 1(c)).

Consequently, a contact between surfaces of almost arbitrary topography can be replaced by an equivalent axisymmetric system. The related restrictions that limit the class of the In Fig. 1 the word "arbitrary" should be understood in that context. The idea of geometric extensions for contact problems widely used is the method of dimensionality reduction (MDR,

[Pop-13], [Pop-15]).

The existing solution for an axisymmetric system with a circular contact area A of radius c, A=c 2 , has the following form ([Gal-61], [Sne-65], [Jäg-95]):

  2 22 0 2 () 1 c z E N c ca d c             , (3)
  22 0 () c cz a c d c        , (4)
        2 2 2 2 2 0 , 1 cr z p pz p E dr c dp r p r                , ( 5
)
where  is the radial coordinate, z() is a function describing the gap between the two bodies, and z'() and z''() are its first and second derivatives, respectively (see Fig. 1(a)). Eqs. ( 3) and ( 4) can be used to define a parametric expression for the normal force-displacement relationship a=a(N) or N=N(a).

It is straightforward to verify that Eqs. ( 3)-( 5) yield the classical Hertz solution [Lan-93] for two equal spheres when

  22 z R R    
is assumed, with R, the radius of the spheres. In this derivation, only the lowest term in the expansion over the small parameter c/R is kept.

It can be numerically demonstrated [Poh-15] that the precision of the REFP becomes worse in the case of a very strong asymmetry in the contact shape i.e. the situation when one dimension (say, x) considerably exceeds another one (y) . For moving along x, the effective  will deviate from the effective  for movements along y-axis.

Another, more detailed form of the REFP for an arbitrary 3D relief with possible multiple contacts is

        ,, r A r S r      (6)       T N A N S   , (7)
      b a A a S   , (8)
and the value of  is close to Eq. ( 2).

The consideration of this assumption certainly deserves some comprehensive comments.

(i) First of all, the above hypothesis can only be made under the standard set of assumptions typical for contact theories of this kind: the bodies are elastic, the normal directions for all individual contacts are aligned and stay aligned when sheared, the local Coulomb's law is postulated, and torsion is absent.

(ii) Strictly speaking, the approach discussed here is only valid for elastically similar materials with the Dundur's constant

        2 1 1 2 2 1 1 2 11 11 GG GG          
equal to 0 (here G 1,2 are shear moduli, (iii) The expressions and notations used in Eqs. ( 6)-( 8) in which the stick area S can substitute the contact area A reflect the acceptance of the assumption illustrated in Fig. 2. Suppose that for a certain applied normal load N=N 1 the contact zone has the configuration shown on the left-hand side figure, with a total area A 1 that, in the absence of tangential action, coincides with the zone of stick, S 1 =A 1 . Upon increasing the normal force to a larger value (N 2 ), the contact zone increases, too, as shown in the figure on the right, and some individual contact spots may merge. If we then apply a tangential force such that the stick area shrinks and forms an area S 2 that equals A 1 , the configuration of the stick zones at the right will perfectly repeat the configuration of the contact zone for N=N 1 at the left, including the number and shapes of all individual contacts.

The general form of the REFP (Eqs. 6-8) assumes a version more suitable for axisymmetric profiles:

        ,, cs         (9) 107       T N c N s   , (10)
      b a c a s   , (11)
Fig. 2. Illustration of an important assumption related to the reduced friction principle for arbitrary multiple contacts. The left set shows a configuration of contact zones for a specific normal force. Upon increasing the normal force and then gradually applying a tangential force, the configuration of stick zones at the right set will fully coincide with the configuration of the contact zones at the left for a certain value of the tangential force.

in which c is the contact radius, s is the stick zone radius,  is the radial coordinate. The dependencies (c,) N(c) and a(c) are provided by the known normal solution. Expressions (s,) N(s) and a(s) just mean that in the above functional dependencies representing the normal solutions c is substituted by s. Comparison of the axisymmetric representation Eqs.

(9)-( 11) and the short form of the REFP Eq. ( 1) suggests the following interpretation of the parameter Q: Q is the virtually reduced normal force at which the corresponding virtually reduced contact radius coincides with the actual stick zone radius. In other words, Q is the stick zone radius in the force-related coordinate  given by the mapping

  c Nc     (12)
in which N(c) is the normal solution, as previously.

The short form of the REFP Eq. ( 1) does not contain any parameters directly related to the contact geometry. All geometric information is contained in the dependency N=N(a) linking the normal force (N) and the normal displacement (a). This dependency is one of forms of the normal solution and is considered to be known. For axisymmetric profiles it is given by Eqs.

(3)-( 4). For rough surfaces, one of the appropriate models [Gre-66], [Per-02], [Biw-04],

N=N 1 , T=0 N=N 2 >N 1 , T0 slip area stick area stick area + slip area = contact area [Pag-10], [Pag-14] can be used. Anyway, the REFP in the form of Eq. ( 1) means that any arbitrary contact satisfying the theorem's conditions behaves in the same way in terms of displacement and forces as an equivalent single-contact axisymmetric system having the same normal reaction.

Method of memory diagrams

The method of Memory diagrams (MMD, see [Ale-15]) do not require any specific loading conditions. In 2D, it provides a two parametric output for a two-parametric input i.e.

    ,, N T a b  for a force-driven and     ,, a b N T 
for displacements-driven systems. 3D

extensions are also possible but are not discussed here. A memory diagram is a simplified form of representing a complex traction distribution by a graphical object that contains the same amount of information. For the force-driven system, it is convenient to use a forcerelated variable , instead of the radial variable , which is defined (see Eq. ( 12)) by a nonlinear scaling transformation

  N    
, where  formally substitutes c in the normal problem solution N=N(c). Aside from the nonlinear mapping  has the same practical sense as the radial coordinate . Using coordinate  does not limit the applicability of results to axisymmetric systems only since, as we have seen in Section 2, an arbitrary contact system for which the REFP is valid can be replaced by an effective axisymmetric system.

Simplest memory diagram for an initial curve

In order to introduce the MMD concept, we start with the simplest constant-loading case: we consider a contact zone which is initially free of traction, and then create a traction distribution by applying a tangential force T, assuming that the normal force is kept constant during the tangential loading. This particular situation has already been described in Section 2 (see Eq. (1) or Eqs. ( 9)-( 11)) and is illustrated for the contact of two spheres in Fig. 3 segment 0<<Q stays unchanged and equals 0. This was a deliberate assumption as the traction distribution in that zone can be easily calculated using Eq. ( 4) with known c and s (or N and Q) and does not contain any additional information. 4) is satisfied, so that slip is absent, and |()|<(c,) everywhere in the contact zone.

In both cases, the forces have been applied to a traction-free contact (absence of loading history, initial curve). Groups (c) and (b): Memory diagrams corresponding to the basic solutions: (c) corresponds to (a), (d) corresponds to (b).

The force balance equation in terms of traction distributions reads

  0 2 c dT        .
In the language of memory diagrams the balance equation becomes

  0 / N D d T      (13)
which can be easily verified, since the integral in Eq. ( 13) amounts to N-Q and therefore equals T/ according to Eq. (1). Thus, the memory diagram shown in Fig. 3 (c), on one hand, contains all elements necessary to calculate the traction distribution in Fig. 3 (a), and, on the other hand, represents an alternative form of the force balance equation.

 c+c () c (a) (c)  () (c+c,) c s (c,)  N 1 D() Q  N+N 1 D() N (b) (d)
constant during the application of T but rapidly increases so that   1/ / 1 dT dN



 . This slipfree case is referred to as "overloading" or "quick" enlargement of the contact zone in [Ale-12] and [Ale-13], or "complete overlapping" by , [Jäg-03]). The slip-free solution can be derived from the same equations Eqs. ( 9)-( 11) as in the previous case using the incremental technique [Min-53]. It consists in assuming a small increment N that enlarges the contact zone from c to a new radius c+c, and a small increment T that produces slip characterized by the radius s. If s lies between c and c+c no slip occurs. Using Eqs. ( 9)-( 11)

we write

      , c c s c c            ,   T dN c s c dc       ,   b da c s c dc       ,
where the last two equations respectively describe slip induced by the force T or by displacement b, applied to a contact of radius c+c. Solving for c-s+c, we finally obtain

      ,, 1
cc dc dc Tb c dN c da              (14)
for the traction increment, and

b da T dN     (15)
for the tangential force and displacement increments, where derivatives dc/dN, dc/da, and da/dN are given by the known normal solution. Besides, the condition representing the absence of slip

c s c c     corresponds to the validity condition   1/ / 1, 0 T N N       (16)
of this solution.

Fig . 3 (b) shows the traction increment corresponding to an exemplar case of the contact between two spheres. In the plot, the increments c and  are not small but obtained as a sum of small increments, each satisfying the validity condition of the slip-free solution so that the absolute value of the traction in Fig. 3 (b) is less than the limiting value  at each point.

In defining the related memory diagram, our intention is to keep the force-balance equation in the form of Eq. ( 13) for any current N. The result is shown in Fig. 3 (d) where the curvilinear piece of the memory function is given by

  1 N dT D dN      . ( 17
)
The absolute value of the memory function on the new segment N<<N+N is less than 1 which follows from Eq. ( 16). This means that the memory diagram does not end with a segment on which D()=1, and indicates the absence of slip in the considered case.

Similar as in the previous situation of slip, there is an unchanged part 0<<N of the memory diagram despite the fact that the traction in the corresponding contact area (0<<c) changes. It is not critical in the MMD since the traction in the whole contact area can be computed from the solution for the normal stress using Eq. ( 14), taking into account only those parameters pertaining to the important part N<<N+N of the diagram. Thus, again, the memory diagram in Fig. 3 (b) contains all necessary information to retrieve the original traction in Fig. 3 (d) and, at the same time, satisfies the force balance equation in the form Eq. (13).

Evolution of memory diagrams

In the general case of an arbitrary loading history, we define the memory diagram as an internal functional dependency D() satisfying the following three properties.

(I) The memory function D() is specified on the -interval limited by the current normal force N (0<<N), and evolves in such a way that its integral over this interval always equals T/ (Eq. 1) in accordance with the force balance equation.

(II) The absolute value of the memory function can not exceed 1. This requirement is related to the Coulomb friction law and agrees with the procedure established above in this section in which   1 D   in the changing part of the diagram, while in the other part D() is left unaffected.

(III) A possible adjustment of the memory function D() in the purpose of complying with the principle (I) is made by setting a final segment Q<<N of the memory diagram on which D   . This final segment corresponds to the presence of slip that propagates inward until, by proper setting of Q, the condition (I) is satisfied.

Fig. 4. Complete algorithm of the method of memory diagrams. The equations displayed in

the boxes correspond to Eqs. ( 18)-( 35).

In the following paragraphs, we explain how the conditions (I)-(III) determine the evolution of a memory diagram for arbitrary force increments N and T. At each instance, the updated memory diagram allows us to calculate the corresponding displacement increments thus providing the solution to the contact problem. The complete MMD algorithm for an individual step, going from the previous (subscript "p") to the next (subscript "n") values of the parameters, is shown in Fig. 4. In our notation, the increments are denoted as  (e.g.

N=N n -N p , T=T n -T p , etc). Since the normal solution to the contact problem is known, the normal displacement increment

    np a a N a N    , (18)
can be calculated directly , while b has to be determined as the result of the procedure. The algorithm involves only two binary choices and, therefore, it contains three cases in which the

Yes overloading   n p N n N T D d      da b T dN     End of current step Take new N=N n -N p , T=T n -T p , a=a(N n )-a(N p ) Check if N>0 Check if |T|/(N)<1 No   1 sgn T T N      2 1 T T T          sgn , n n D T Q N         2 p N Q T D d          1 sgn b T a        2 p N Q da b D d dN         1 2 b b b      No   1 n p N p N T D d      2 1 T T T          2 sgn , n n D T Q N         2 n N Q T D d          1 n p N p N da b D d dN         2 p N Q da b D d dN         1 2 b b b      Yes YY YN N 113
calculations differ. As illustrated in Fig. 4, we have called these cases YY, YN, and N. In all illustrations, the values of "small force increments" are exaggerated. (YN) which will be discussed later. Indeed, if the area |T| corresponding to T is small enough i.e.

 1 D() (a) (b) -1 Np Nn  1 D() -1 Np Nn (c)  1 D() -1 Np Nn T/ T1/ T2/ A' A T=T1+T2 Q  1 D() (d) -1 Np Nn (e)  1 D() -1 Np Nn -T1/ T2/ A' A T=T1+T2 Q 0 N  0 N  YY YN N   1/ / 1 TN     , (19)
it is possible to define D n () on the new segment N p <<N n such that the force balance equation (1) will be satisfied without any changes to the memory function on the old segment

0<<N p :   / n p N n N D d T      , (20)
as illustrated in Fig. 5 (b). For small increments, linearly approximating the integrand in Eq. ( 20) we write Thus, the condition (I) is already fulfilled, and the use of the rule (III) is not required.

    / 2 n p n n D N D N NT      , (21)

The inequality

  1 nn DN  means that the end of the diagram does not contain a section on which D() equals either +1 or -1 which corresponds to slip-behavior. In other words, in this situation, the tangential force change can be "absorbed" by the system without producing slip, i.e. by a purely elastic deformation of the contacting bodies, and consequently, the slip-free solution Eqs. ( 14), ( 15) is applicable since condition Eq. ( 16) is fulfilled. Hence, the corresponding solution for b is given by Eq. ( 15).

  2 / N Q D d T        , ( 25
)
where Q is the new boundary between the stick and slip zones in the -space, and D n ()=sgn(T)=sgn(T 2 ). In fact, Eq. ( 25) should be considered as an equation for Q.

Once T 1,2 and Q are determined, the solution for the displacement increment corresponding to T 1 easily follows from Eq. ( 11) and reads

  1 sgn b T a      , ( 26
)
while the remaining part b 2 is given by

  2 p N N Q da b D d dN          , ( 27
)
where Q is known (see further comments on Eq. ( 27) below). Finally, it suffices to sum both contributions to obtain the total displacement increment: Case N. We now consider the situation in which the normal force is constant or decreasing, 0 N . In this case, the memory diagram should shrink to the new size N p , and therefore release some tangential force T 1 , as shown in Fig. 5 (d):  

1 n p N p N T D d      , (29)
(T 1 =0 in a particular case when N=0). In order to balance the force equation, we have to account only for the partial force increment 21

T T T      . ( 30
)
As previously, there is no other option than to allow slip to penetrate closer to the contact center (except in a particular case when T 2 coincidentally equals 0). This means that we have to shift point A, starting from position A' (Fig. 5 (e)), while setting

    2 sgn n DT   , (31)
on interval AA', until the force balance equation

  2 n N Q T D d        , ( 32
)
is satisfied. As previously mentioned, various slip scenarios are possible depending on the structure of the memory function on the segment 0<<N n .

Once the values of T 1 and T 2 are determined for this case, the related b-increments

then read   1 n p N p N N da b D d dN        , ( 33
)   2 p N N Q da b D d dN          , ( 34
)
with the total

12 b b b      . ( 35
)
Eqs. ( 29)-( 35) then represent the full solution to the problem in case N.

The solutions that allow slip (i.e. hysteretic), as found in the cases YN and N, have an interesting feature. In the case YN when the normal force increases, the maximal setting |D()|=1 in the new segment is not enough for equilibrating the tangential force, so that we have to engage the older segments. As T 1 and T 2 have the same sign (which also equals the sign of T), the direction of slip always coincides with the direction of the tangential force increment. However, in the case N, corresponding to normal unloading, this is not always true.

Decreasing the normal force releases some force T 1 whose value and sign are defined by the system's history, and therefore completely independent on the increment T. In the situation presented in Fig. 5 (e), the direction of slip is positive since D n ()=+1 at the end of the new diagram, but the released area T 1 / is larger than the area T 2 / generated by slip. This implies that the sum T=T 1 +T 2 is negative, and the slip direction is opposite to the direction of the tangential force change.

The above rules composing the solution procedure only represent a first-order approximation since, in particular, Eq. ( 21) has the precision of   oN  . In addition, the instruction "shift point A until the force balance equation is satisfied" actually implies an implicit, and therefore in practice, approximate character of the procedure. However, despite the fact that the procedure of building up memory diagrams is approximate, the solution for the displacements corresponding to a predetermined memory diagram is exact. This is elaborated on in the next paragraph.

"Reading" memory diagrams

In the previous section, we have formulated the evolution rules for a memory diagram governing the behavior of the force-driven system. The choice of arguments (forces or displacements) depends only on the context in which the contact problem is posed, and does not affect the physics of the contact interaction. In particular, the following reasoning is valid.

Any complex traction distribution within the contact zone corresponding to a certain local tangential displacement field can be represented by a memory diagram D() as described above. Based on D(), we can then use Eq. ( 13) to calculate the tangential force. On the other hand, the solution represented in Eqs. ( 10)-( 11) warrants that the MMD is completely symmetric with respect to the argument choice. Therefore, it is possible to introduce a memory diagram D() in the space of a displacement-related argument parameter  given by =a() where  formally substitutes the argument c in the normal problem solution a(c).

Then, similarly to Eq. ( 13), the tangential displacement is given by the following simple relationship

    0 / a D d b      . ( 36
  0 a a dN T D d da        , ( 37
)   0 N N da b D d dN        . ( 38
)
Obviously, this also leads to the following symmetric representation using the original radial coordinate :

  0 c c dN T D d dc        , ( 39
)   0 c c da b D d dc        . ( 40
)
With this representation, the rest of the MMD equations have to be modified accordingly.

To give an example, in Fig. 6 we plot a traction distribution for a contact between two spheres. In conclusion, the integral formulations Eqs. ( 13), ( 36)-( 40) present a simple way of "reading" the memory diagrams in , , and -spaces. In all cases, the memory functions consist of constant or curvilinear sections, as shown in Fig. 56. The difference only resides in a nonlinear stretching of the horizontal axis.

Numerical implementation and examples

The MMD formulation discussed in the previous paragraphs does not impose any requirements with respect to the numerical implementation procedure. The simplest numerical implementation can consist in introducing a fixed-point grid  i on the -axis and in defining the corresponding function values D( i ) on that grid. However, this method is time consuming and inaccurate. We have implemented a more precise technique that uses the fact that the memory function is constant on certain intervals and thus can be represented only by the pairs of the two boundary values of  defining the intervals. For curvilinear sections, not only the interval boundaries but all intermediate points =N(t i ) are to be memorized (here N(t) is the normal force loading history and t i are discrete observation moments, N=N(t i+1 )-N(t i ), T=T(t i+1 )-T(t i ), etc.). Thus, in general, the function D() can be defined on a non-equidistant and adaptive grid. In those cases where the memory function has to be determined in between

 D() (a) N 1  () (b) 2 3 4 5 6 c 1 3 2 a a
two points of a curvilinear section, we use a linear interpolation. As a result, the complexity of this numerical code exceeds only slightly the complexity of the method itself (see Fig. 5). Fig. 7. Tangential force displacement relationship in which three regimes (cases YY, YN, and N) can be identified (gray, thick black, and thin black lines, respectively). In the inset: loading protocol i.e. time dependencies of normal and tangential forces; the regimes are marked as in the main figure . Four moments of time are selected; memory diagrams at these moments are shown in sets (b)-(e). In all plots, N 0 is a characteristic value of the normal force; a 0 =a(N 0 ).

A simple example illustrating all three cases (or regimes) YY, YN, and N is shown in Fig. 7. for a monotonous time dependence of the tangential force T(t) (see the inset in Fig. 7 (a)), the response b(T) (main Fig. 7(a)) is monotonous as well. However, the normal force N(t) in Fig. 7(a) is not monotonous that makes it possible to see various regimes in the corresponding memory diagrams (see Fig. 7(b-e) corresponding to moments of time t=0.25, 0.4, 0.5, 0.9 in arbitrary units, marked by the dotes lines in Fig. 7(a)). The curve N(t) begins with an increasing segment on which the condition Eq. ( 19) is fulfilled (case YY, thick gray lines in Fig. 7(a)). According to the given explanation, a curved-line segment appears in the memory diagram (Fig. 7(b)). Further, N(t) still increases but the condition Eq. ( 19) is not satisfied anymore thus resulting in appearance of the regime YN (thick black line in Fig. 7 In a more general case of two independent arguments, N and T without functional relationship, the hysteretic behavior differs considerably. Since variations in N are not linked to the T(t)protocol, it occurs that, even for a periodic T(t) loading history, the "loops" are not necessarily closed, on the contrary. Indeed, as illustrated in Fig. 8 (b), the same T(t)-history as used in Fig. 8(a) produces a curve in which all monotonous parts are shifted, bent, etc. It is important to emphasize that the generation of such curves via the direct analysis of the traction and Finally, it can be analytically verified that all classical results, as discussed in [Min-53], for spheres loaded by an oblique force follows in a simple and straightforward manner from the MMD equations. As such we have validated our novel approach in numerous exemplary situations.

Summary: assumptions of the MMD

Chapter IV is devoted to derivation of fundamental solutions for frictional contact mechanics applicable for arbitrary loading protocols. Our objective was to find a solution (i.e. normal (a)

and tangential (b) displacements ) to the frictional contact problem with an arbitrary geometry of contacting profiles and for an arbitrary loading history in terms of the normal (N) and the tangential (T) forces. The solution is given via the method of memory diagrams that uses the following assumptions.

1. Loading is in one plane only (i.e. in 2D) and quasi-static.

2. The Coulomb friction law with friction coefficient  is postulated for contact stresses.

3. |T|< N so that only partial slip is considered; the opposite for quasi-static loading is trivial.

4. Plasticity and adhesion are ignored.

5. All individual contact spots are aligned (i.e. they have the same normal directions) and stay aligned during loading.

6. The normal solution a=a(N) is independent of the tangential loading and is known from previous studies.

7. Dissimilarity effects are neglected.

8. The reduced friction principle is valid for the considered geometry.

The method is valid for a wide range of contact geometries (see section 2). For instance, multiple contact spots can merge and split. The analysis is essentially based on the reduced friction principle which allows us to replace the arbitrary body shapes by a pair of effective axisymmetric profiles having the same normal reaction a(N) which is considered to be known. Once the effective profile is determined, we use the MMD for the effective axisymmetric contact and deduce a computer-assisted analytical solution. The method is formulated both for force-driven and displacement-driven systems.

Chapter V. Modeling for elastic wave propagation in materials with cracks

Force-driven and displacement-driven crack models

Any successful simulation for the elastic wave propagation in a material containing cracks needs a relevant crack model. Such a model should be formulated for a small representative volume of material intersected by a crack. In other words, we introduce a mesoscopic cell in which boundary conditions at crack surfaces are to be defined. Here the term "mesoscopic" means that the cell size is much less than a typical wavelength or other characteristic scale of the macroscopic elastic field, but at the same time the cell can host a large number of microscopic features e.g. asperities present on crack surfaces. The mesoscopic level description includes average stresses and displacements linked with a relationship we are looking for. The normal and tangential displacements denoted here a and , and normal and tangential stresses ( N and , respectively) are considered as lumped but not field variables, in contrast to local (microscopic) stresses and displacements. It is convenient to refer to N and as to forces per unit of nominal contact area or just forces in order to easily distinguish between microscopic field variables and mesoscopic lumped ones. In addition, " forces N and " are traditionally used in contact mechanics since 1950s [Min-53].

The crack model should provide the boundary conditions that represent a link between forces and displacements at crack faces. From the point of numerical simulations, the crack model should exchange force-displacement data with a solid mechanics unit that solves the elasticity equations in the bulk volume. Generally, both crack model and solid mechanics unit can be force-driven or displacement-driven. A possible situation when the crack model is driven by displacements and the solid mechanics drive parameters are forces is illustrated in Fig. 1.

Suppose that the crack model is capable of calculating forces for given displacements in some way. This means that stresses are determined at any relevant point of the structure, and therefore the solid mechanics unit can calculate the proper strains at a new time step. Then, the information about the strain field in the whole material allows one to obtain displacements at the crack faces, and so on.

In Fig. 1 the data exchange is organized in a simple explicit manner. However, as we will see, it is not always the case. In particular, consideration of flat crack surfaces with friction and without adhesion results in an iterative procedure illustrated in Fig. 2. The implicit character of the solution is due to the following reasons. Let us first accept for definiteness that N>0 for compression and, similarly, a<0 for an open crack whose faces do not touch each other at a particular mesoscopic cell. Then, for the normal interaction, two situations are possible: N>0 when the crack is closed (full contact in that cell) and therefore a=0, and when the crack is open i.e. N=0 and a is undefined. Indeed, N=0 means the absence of interaction and corresponds to any crack opening a<0; the actual displacement a has to be determined from external conditions. The third situation with N<0 is not possible since no attraction force between the faces exist in the absence of adhesion. An analogous behavior occurs for the tangential interaction. According to the Coulomb friction law, means stick when the tangential displacement does not change. In the case the crack faces slide so that is undefined in the framework of the crack model and is to be determined from external conditions. Finally, is not possible.

  N N a    , T T b a  Tangential bad adjust , NT , NT CRACK MODEL  N>0 closed  a=0  N=0 open  a<0 undefined  N<0 not possible  TN   stick  b unchanged  TN   sliding  b undefined  TN   not possible Tangential Normal SOLID MECHANICS Input , ab Output , NT , NT new good accept , NT , ab , ab 126
Thus, the force-driven model provides no response in the cases N<0 and . However, a solid mechanics unit that works independently of a crack model can easily provide such data for forces, notably when at a given point an initially closed crack is about to open or when initially stuck faces are about to slide. In these cases, the crack model should "ask" the solid mechanics unit to gradually change a to small negative values or slightly change in order to obtain N=0 or . Since all points in the material are linked via solid mechanics equations, such an operation cannot be done locally and will involve redistribution of stresses and strains in the whole volume. As a result, a complex iterative procedure is required that simultaneously satisfies all slip or contact loss conditions at all mesoscopic cells.

Since the simple force-driven model of the crack with flat faces (Fig. 2) actually entails cumbersome computations, it is of interest to study the case when the faces are not flat. A plausible example would be a crack with rough surfaces. Further, we show that such a modification produces an explicit data exchange procedure in the spirit of Fig. 1, although it makes the crack model more complex.

Normal loading curves for contact of rough surfaces

      2 0 21 a hE N a z dz        . (4)
Eq.( 4) suggests that the normal reaction of a crack section belonging to a considered mesoscopic cell is determined by the gap (aperture) distribution which, in turn, depends on the nature of a crack. Since typical acoustical excitations are can always be considered as a perturbation, the normal reaction in the acoustical strain range mainly depends on the tangent to the curve (z) at z=0. Here three cases can be distinguished: the tangent is vertical, horizontal, has inclination angle between 0 and /2. The vertical tangent means that even for a small compressive displacement a a non-zero contact area will be immediately formed. In roughness brought into contact. We also assume that the Eq. ( 5) is approximately valid for fatigue cracks since the internal stresses release during cracking and associated microscopic displacements result in the loss of conformity at the atomic scale. Therefore low-scale composite roughness that mainly contributes into acoustical response can be considered as uncorrelated. Thus, in this study, we accept Eq. ( 5) as a model equation for the normal reaction curve, although for the proposed theory form of N(a) is not essential.

Certainly, a tangential excitation may considerably displace the contacting asperities in the lateral directions and therefore alter the normal reaction curve N(a). However, this effect can be neglected if we assume that the random gap is a stationary process whose characteristics '

h and   z 
do not depend on the tangential shift. Therefore N(a) is supposed to be a portrait of the system which is not affected by any other interactions and incorporates all geometric information about rough surfaces in contact.

Tangential contact interactions: full sliding and partial slip

The objective of this section is to build up the crack model that works in the spirit of Fig. 1 and does not produce any iterative data exchange procedure. To do this, we use the established normal and tangential load-displacement relations; however, one point is still missing. The matter is that the MMD is only applicable for partial slip, and the case of total sliding remains unaddressed. In order to complete the description, the total tangential displacement should be split into two parts: one, b 0 , corresponding to total sliding and the other one, , to partial slip, .

The idea can be illustrated with the following example. Suppose that tangential loading increases under a fixed normal compression N=N(a) corresponding to a certain normal displacement a. Asperities recede under load in both normal and tangential direction. In addition, the tangential receding , in contrast to the normal one, a, is accompanied with partial slip. According to the MMD or eventually to the Coulomb friction law, can not grow infinitely. Once the maximum value is achieved the asperities can not deform anymore. Then the total sliding process develops where the very last stick point belonging to one face separates with its neighbor at the opposite face. The displacement between those points is denoted b 0 .

129 Certainly, when the normal compression changes, changes too. In addition, the contact at a current mesoscopic cell can be lost at some moment of time so that both stick and slip zones disappear. In a general case when a and b evolve in an arbitrary manner, the concept is illustrated in Fig. 3. This figure explains how to calculate the forces N and T for any given values of displacements a and b. In order to do that, tangential displacements components b 0 and together with the corresponding memory diagram should updated. The update operation for tangential displacement components is denoted using the assignment operator ": " which means that values at a previous step have to be replaced by new ones.

The model of cracks with rough faces assumes three possible contact states: contact loss when no points belonging to the opposite faces are in contact, total sliding when at each contact point slip takes place and TN   according to the Coulomb friction law, and the remaining case of partial slip. Below we consider these cases in more detail.

(i) If 0 a  the contact is lost which obviously means that N=T=0. The repartition in Eq. ( 6) in that case has no sense as well as the memory diagram. However, it useful to formally define these characteristics even at the absence of contact having in mind that at the next time step contact can be reestablished. Then the process should start with a "virgin" memory diagram D()=0 since the contact zone contains no residual tangential stress. The asperities are not strained at this moment i.e.

, hence b 0 =b. Accepting these modifications will guarantee the correct evolution representation once the crack faces will ever get in contact.

(ii) Suppose now that 0 a  and the old (i.e. obtained at the previous step) value of is such that with a new a. Then, the new value of should obviously be corrected, since the maximum tangential displacement corresponding to the elastic deformation of asperities can only be , and the attempt to further increase the tangential action will result in total sliding. In this case we have to set the new with the sign corresponding to the direction of sliding. The remaining part of the tangential displacement corresponds to the total sliding contribution, . In other words, the reference point for measuring the tangential deformation of asperities is shifted. In accordance to the Coulomb friction law, TN   and   1 D   with the same sign corresponding to the sliding direction. Here the magnitude of the normal force N=N(a) is selected at the known normal reaction curve. 130 (iii) Assume now that 0 a  and . In that situation, some points of the contact zone stick, therefore slip, if it occurs, can only be partial. Then the reference value b 0 is not affected which is symbolically expressed as assigning 00 : bb  . Obviously for the remaining part of the tangential displacement. In this regime, the MMD should be used with a given as an argument i.e. . Again, N=N(a). Fig. 3. Three possible contact states in the model of cracks with rough surfaces. For each case, the following information is supplied: conditions under which the case occurs; solutions for components b 0 and b , solutions for forces N and T, memory diagrams.

In the latter case there is, however, the risk that the new will exceed a, even though the old does not. Then the solution should be taken as in the total sliding regime. Such a situation should be additionally checked for.

The algorithm in Fig. 3

0 a  0 a   1 D() -1  1 D() -1 a  D() a   T MMD b  1 -1
131 condition for the displacements in the form instead of the traditional form TN   written for forces. The problem with this approach however, is that this contact formulation is strictly valid only for stationary problems. In case of dynamic problems, such as the elastic wave propagation problems we are interested in, the contact pairs can only be used as long as situations with impact are avoided. This is certainly not the case for cracks, which can exhibit clapping behavior. Moreover, when friction also needs to be included in the contact pair, it becomes even more difficult to obtain a stable (transient) solution.

Numerical implementation of wave-crack interactions

The above feature is a strong argument supporting the usefulness of the crack model formulated in sections 2 and 3 of this chapter. The disadvantage of this crack model, however, is that it is too complex to be incorporated directly into COMSOL. Luckily, this can be circumvented by first implementing the full crack algorithm in MATLAB and then connecting the obtained MATLAB function to COMSOL using the LiveLink for MATLAB

[COM-15b].

Using this approach, the interaction of an elastic wave with a frictional contact interface can be studied using the following steps:

(a) In the Structural Mechanics Module of COMSOL, normal and tangential stresses ( and  ), defined by appropriate internal and external boundary conditions, are used to calculate normal and tangential displacements (u n and u t ) at a particular time step t of the procedure.

The relative normal and shear displacements calculated at the Gauss points (or integration points) on the contact interface are outputted to MATLAB. The link between the mesoscopic stresses  and , relative contact displacements u n and u t , on one hand, and the variables of the contact model is given by

    2 , 2 , , , . nt u a u b N a T a b          (7)
Based on Eq. ( 7), the normal reaction curve N(a) is considered to have a quadratic dependency:

    22 0 N a C a a  , (8)
where C = 6 × 10 10 Pa -1/2 m -1 , similar to the value used in the quadratic pressure-overclosure relation used in the numerical study of contact between two solid blocks of aluminum by

Illustrative example

In order to demonstrate the potential of the above described model, an illustrative example of a shear wave propagating in a 2D aluminum sample with an inclined crack is studied. The study allows to demonstrate the influence of both normal and tangential interaction of the crack faces on their respective displacements, as well as the presence of nonlinear features generated at the crack interface due to these interactions. Finally, Fig. 7 shows the results obtained when using the full crack model with both the normal and tangential forces included. As can be seen, the normal relative displacement and normal contact stress did not change when compared to the previous nonlinear model (Fig. 6).

The tangential behavior, on the other hand, is clearly changed due to the introduction of T in the model. Eq. ( 7) states that the normal force T is depending on both the normal and tangential displacements, with the relation being determined by the algorithm represented in Fig. 3. In case of an open crack (i.e. a < 0 or u n > 0), the tangential force T, and therefore also the tangential stress  , is equal to zero. In case of contact (i.e. a≥0 or u n ≤0), the tangential force differs from zero, with the value depending on the contact state: total sliding or partial slip. In the partial slip case, the tangential force (and tangential stress) is determined by the MMD algorithm, whereas the total sliding case states that the size of the tangential force (tangential stress) can never exceed μ times the size of the normal force (normal stress).

This behavior is clearly visible in Fig. 7 where the tangential contact stress  indeed only differs from zero in case of contact, with the value never exceeding μ times the normal contact stress. Again, the influence of introducing the tangential force T is also obvious in the relative displacements plot, where the relative tangential displacement is changed when compared to the previous nonlinear case (Fig. 6). It can be seen in Fig. 7 that the relative tangential displacement amplitude decreased in case of contact. This is due to the presence of the tangential contact stress, which is linked to friction effects that will reduce the relative distance the crack faces can move in tangential direction. shows a quadratic relation between the normal contact stress and the relative normal displacement, as expected according to Eq. ( 8). The tangential reaction curve, where switching between the three defect states (contact loss, total sliding and partial slip) can be identified, clearly indicates hysteresis behavior with several hysteresis loops due to the presence of friction occurring in the total sliding or partial slip regimes.

Clapping-and friction-induced nonlinear features

The above discussed results illustrate that the proposed crack model is able to simulate the expected behavior of closed cracks with friction. This behavior is mainly characterized by the presence of normal and tangential stresses at those positions where the crack faces are in contact. Since these contact stresses are directly linked to the relative displacement of both crack faces, the movement of both faces will be distorted, resulting in nonlinear features generated at the crack interface, such as for instance the generation of harmonics. present. In order to highlight the nonlinear components, several techniques can be used. One possibility would be to subtract the calculated relative displacement signals obtained using the linear (reference) model from those obtained using the crack model, in order to eliminate all linear components. However, since this method requires having two similar samples, one sample containing a nonlinear crack and another reference sample containing a linear crack, this method is not suitable to be used in real experiments. We therefore prefer to use another method, the Scaling Subtraction Method (SSM) [Sca-08],

based on the distorted scaling of the received signals with increasing excitation amplitude due to nonlinearity. This method only requires the sample with the nonlinear crack which is then excited twice, one time at a low excitation amplitude A low and one time at a high excitation amplitude A high = n A low . By subtracting the scaled calculated relative displacement signals obtained using the low excitation amplitude from the ones obtained using the high excitation amplitude, the linear contribution in the signals will also be eliminated and nonlinearities will be enhanced. Fig. 10 shows the results when using this approach in case of A low = 10 nm and In summary, by studying the time evolution of the nonlinear content in the wave propagation, it is obvious that due to the wave-crack interaction, the crack starts to behave as a nonlinear source. The signals excited by this nonlinear source can for instance be used for defect detection, localization and/or characterization.

A

Conclusions and perspectives

This chapter is concerned with theoretical modeling for waves propagating in a linear elastic material that contains cracks with friction. Contact interactions introduce nonlinearity to the system; in particular, friction makes the nonlinear behavior hysteretic and memory-dependent.

The principal idea of the approach is to describe frictional interaction on the basis of classical physical models such as Coulomb friction law and, correspondingly, avoid phenomenological hypotheses on contact acoustical nonlinearity. A successful description for contact interactions would produce vector load-displacement relationships at internal contacts in the material. The boundary conditions at internal and external boundaries are both necessary for wave propagation and vibration modeling.

Our approach includes a concept of a mesoscopic cell (physically small volume containing a small fragment of a crack). In each cell, contact stresses and displacements are supposed to be uniform. In other words, we consider an auxiliary problem of contact between two bodies loaded by remote forces and seek for a link between lumped parameters instead of elastic fields. This feature makes our approach essentially different from purely numerical methods of contact mechanics in which all structure elements should be properly meshed and no intermediate scale is present.

An attempt to implement the Coulomb friction law as a boundary condition is related to a certain difficulty. The matter is that the Coulomb friction law does not produce an explicit link between contact stresses and displacements. The problem can be solved by adding some phenomenological assumptions or by implementing an iterative procedure that matches stresses and strains in the material in order to satisfy the Coulomb condition for sliding.

However, in this work, another possibility is explored; we introduce roughness to the crack model that makes it possible to account for partial slip regime occurring when some contacting points slide and some do not. In this situation, the previously developed method of memory diagrams is applicable. Its main advantage is that it allows one to link directly contact loads and displacements by expressing them through an internal functional dependency (memory diagram) that stores all memory information in the system. The loaddisplacement relationship produced by the MMD can be combined with analogous solutions available for two other contact regimes, total sliding and contact loss. As a result, the hysteretic vector link between loads and displacements is calculated by a procedure involving the MMD.

The established boundary condition completes the description for elastic waves in materials with cracks. Technically, the code for modeling waves is programmed using an available finite element software (COMSOL) and allows one to calculate and to visualize all elastic fields in the sample.

An important perspective of this study is related to the comparison of modeling results to experimental data. A sample with engineered contacts at known positions can be used for measurements via one of the techniques described in Chapters I-III.

In addition, a known dynamic response of a frictional internal contact makes it possible to calculate instantaneous energy loss due to friction and use it as a heat source for modeling for temperature effects. Modern finite element software packages such as COMSOL Multiphysics have a possibility to add a heat transfer module to the solid mechanics unit and calculate temperature fields generated by frictional losses. This study could be applied as modeling support for the methods of ultrasonic thermography (see example in Section 3.4 of Chapter I). Résumé en anglais : This PhD thesis work contributes to the development of nonlinear elastic methods for non-destructive testing and imaging of contact-type defects in solids. In this work, two modifications of recent nonlinear nondestructive testing methods are suggested: the coda wave interferometry combined with the nonlinear time reversal principle and aircoupled nonlinear ultrasonic imaging. The principal advantage of former technique is in its extremely high sensitivity owing to the fact that weak changes in sample's parameters are accumulated during the formation of the coda wave. The other technique has a complimentary strength and offers a possibility of a remote detection. The developed techniques are tested on samples with artificially fabricated defects at known locations and their performance is accessed. The potential for obtaining robust nonlinear images is demonstrated. The second part of the work is concerned with a theoretical description of contact acoustical nonlinearity and its use for creating of a numerical toolbox capable of simulating wave propagation in complex structures containing internal contacts. A physical model describing the tangential shift of two contacting bodies in the presence of friction has been proposed. Its result is an analytical computer-assisted solution for hysteretic relationships between normal and tangential contact displacements and loads. The contact model and derived load-displacement relationships are used as boundary conditions posed at the internal boundaries (contact surfaces) in a finite element wave propagation model programmed via commercial software.

Mots-clefs : Nonlinear acoustics, Coda wave interferometry, contact nonlinearities

  d'un point de vue expérimental de développer une technique de CND par acoustique non linéaire suffisamment sensible et fiable afin de détecter des défauts « cachés » (non détectable par des techniques d'acoustique linéaire), permettant de faire des images et fonctionnant dans un environnement industriel ; -d'un point de théorique de développer un modèle (ou même un outil) numérique pour la modélisation de la propagation d'onde acoustique dans un milieu contenant des défauts prenant en compte des géométries complexes et reposant sur des modèles physiques réalistes des contacts.

(

  potentiel de Lennard-Jones [Len-24], modèle de Landau [Lan-93], …), et les non linéarités de contact qui apparaissent lorsque le matériau contient des contacts internes comme des fissures, des délaminations, … Pour les non linéarités de contact, les modèles existants sont classés en deux catégories : -Les modèles phénoménologiques comprenant le modèle de « clapping » [Sol-02], le modèle de friction glissante [Bal-02], le modèle de Preisach [Pre-35], [Kra-89], les non linéarités hystérétiques quadratiques [Ale-04], les modèles de Nazarov [Naz-03], Davidenkov [Dav-38] et Granato-Lücke [Gra-56], les modèles de Preisach vectoriel [May-88] et tensoriel [Hel-01], et le modèle de Preisach-Arrhenius [Tor-02], [Gus-05]. -Les modèles physiques comprenant le modèle de Lawn et Marshall [Law-98], le modèle micro-potentiel [Ale-05], le modèle d'adhésion [Sha-94], [Ale-07b], le modèle « Soft-ratchet » [Vak-05], et un modèle reposant sur un ensemble de billes de Hertz-Mindlin [Nih-00].

  s'affranchir des variations qui pourraient provenir d'autres effets que la présence d'un défaut (comme une variation de température, une modification de la géométrie de l'échantillon,…) une version non linéaire de la technique CWI est utilisée [Zha-13]. Dans cette technique, qui repose sur la mesure des variations du décalage temporel ε et de l'amplitude de l'enveloppe Kd de l'onde de coda du fait de son interaction non linéaire avec une onde de pompe basse fréquence d'amplitude croissante, aucun signal de référence réalisé dans un échantillon sans défaut n'est nécessaire. Les variations ε et Kd sont alors des indicateurs de la présence d'un défaut. En effet, en présence des non linéarités de contact des défauts l'onde pompe vient moduler l'onde de coda, alors que dans un échantillon intact aucune interaction ne sera observée. Malheureusement, cette technique CWI ne permet pas de localiser le défaut. Pour palier à ce problème, on concentre l'onde pompe dans une zone localisée de l'échantillon à l'aide du principe de retournement temporel [Fin-89], [Fin-92], [Fin-96] réalisé grâce à un transducteur à cavité chaotique [Mon-04], [Bou-09]. Cela permet, comme le montre la Fig. 1, de réaliser une analyse CWI localisée. La possibilité de contrôler le retard de l'onde sonde par rapport à l'onde pompe est aussi introduite afin d'étudier les effets de dynamique lente.

Fig. 1 .

 1 Fig. 1. Principe de la technique de CWI non linéaire pour la détection de défauts utilisant une onde de pompe focalisée en différnents points de l'échantillon.

  a été développé pour contrôler la focalisation par le transducteur à cavité chaotique. Ce programme permet de tester différentes méthodes de traitement de signal afin d'améliorer le rapport signal sur bruit et la qualité de la focalisation : compression à l'aide de sweep [Mis-05], retournement temporel, retournement temporel sur 1bit [Mon-01], filtre inverse [Zve-04], [Qui-04], et filtre inverse sur 1 bit. Le filtre inverse, dans la version proposée par Quieffin et al. [Qui-04], permettant d'obtenir une meilleur qualité de focalisation et un meilleur contrat, est sélectionné. La bande de fréquence pour l'onde pompe focalisée est choisie entre 50 et 100 kHz, afin d'être découplée de la bande de fréquence choisie pour l'onde sonde, à savoir de 300 kHz à 1.2 MHz. Deux techniques utilisées pour l'analyse des signaux CWI, i.e. la technique dites du doublet (qui utilise le coefficient corrélation par décalage temporel) et la technique d'étirement, sont décrites en détails afin de montrer leurs avantages et inconvénients. La technique d'étirement plus précise est choisie pour traiter l'ensemble des mesures réalisées. Le choix de la fenêtre temporelle utilisée pour le traitement est explicité. Il nécessite entre autre une mesure précise du rapport signal sur bruit sur les signaux de coda qui est réalisé sur chaque mesure.Le montage expérimental mis en place est présenté sur la Fig.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Schéma du montage expérimental pour l'imagerie de défauts par la méthode de CWI non linéaire.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Protocole de mesure utilisé pour le test avec la méthode CWI non linéaire de l'échantillon de verre fissuré thermiquement. Résultat de l'analyse CWI en mode continu : évolutions (b) du paramètre d'étirement ε, et (c) du coefficient de decorrélation Kd. Une mesure comparative sur un échantillon de verre intact est aussi présentée.

Fig. 7 .

 7 Fig. 7. (a) Principe de localisation de défaut combinant la technique CWI et le concept de retournement temporel. Evolution (b) du paramètre d'étirement ε et (c) de Kd quand l'onde pompe n'est pas focalisée (en vert), focalisée sur la fissure (en noir), ou en dehors de la fissure (en bleu).

Fig. 8 .

 8 Fig. 8. Schéma du dispositif expérimental utilisé pour l'imagerie non linéaire de défauts par la méthode d'imagerie non linéaire par couplage air utilisant un post traitement par mise à l'échelle et soustraction.

Fig. 10 .

 10 Fig. 10. Comparaison des images obtenues sur une plaque de composite CRFP contenant une délamination artificielle de 6 mm de rayon (a) par une mesure par ultrasons aériens en transmission et (b) par la méthode NACU utilisant un post traitement SSM.

  Fig. 11 et l'algorithme de la méthode MMD sur la Fig. 12.

Fig. 12 .

 12 Fig. 12. Algorithme complet de la méthode des diagrammes de mémoires.

Fig. 13 .

 13 Fig. 13. Courbes d'évolution du déplacement tangentiel en fonction de la force tangentiel pour deux chargements temporels différents présentés en insert. (a) Une hystérésis pour laquelle N et T sont liés par une dépendance fonctionnelle. (b) Deux hystérésis pour lesquelles N et T sont indépendant. N 0 est une valeur caractéristique de la force normale; a 0 =a(N 0 ).

Fig. 14 .

 14 Fig. 14. Echanges de données entre le modèle de fissure et le module de mécanique du solide dans le cas d'un modèle de fissure piloté par les déplacements.

Fig. 15 .

 15 Fig. 15. Géométrie et maillage de l'exemple choisi.

Fig. 16 .

 16 Fig. 16. Déplacements relatifs u n et u t (au dessus) et contraintes de contact  et  (en dessous) au point central de la fissure pour une excitation par une onde de cisaillement de fréquence 100 kHz et d'amplitude A = 100 nm.

:

  develop a nonlinear acoustical NDT technique  sufficiently sensitive for robust detection of weak or hidden damage  having the potential for imaging  explore the possibility for remote detection and real field applications Theoretical: create a numerical model or, eventually, numerical tool for modeling wave propagation in materials containing defects  taking into account real complex geometries of samples  based on physically plausible contact models Practical:  contribute to the creation of novel self-repairing aeronautical materials Dissertation structure This document is organized as follows. It contains five chapters, including the introductory one (Chapter I) in which the concept and models of contact acoustical nonlinearity are discussed as well as the most recent existing nonlinear NDT methods. The content of the Chapters II-V is original. The first two of them concern two experimental techniques: nonlinear coda wave interferometry and nonlinear air coupled ultrasonic method. In these chapters, the methodology (background, principle, strengths), experimental setup, measurements results are explained in detail. Chapters IV and V are related to a theoretical development. In particular, in Chapter IV a contact model based on roughness and friction is introduced. Chapter V contains the description of the implementation procedure and demonstrates how the contact model was integrated into a standard commercial finite element software (COMSOL). The last part of the manuscript contains summary, conclusions, and perspectives. Each chapter contains sections numbered 1, 2, 3. etc. Some sections include subsections 1.1, 1.2, 1.3, etc. Equations and figures are numbered in consecutive order as Eq. (5), Fig. 7, etc, within each chapter. In cases when it is necessary to refer to a figure or equation from another chapter, they are referred to as Eq. (I.12), Fig. IV.3, etc. Literature references in the document are cited as [Ada-95] for a paper published by Adams et al. in 1995. The complete list of cited references can be found at the end of the manuscript. Frequently used terms, abbreviations, and variables Terms and abbreviations CFRP : carbon fiber reinforced plastic; CWI : coda wave interferometry; GFRP : glass fiber reinforced plastic; LDR : local defect resonance; MMD : method of memory diagrams; NACU : nonlinear air-coupled ultrasonic imaging; NDT : non-destructive testing; PZT : piezoelectric transducer or lead zirconate titanate; SSM : scale subtraction method; TR : time reversal; coda wave : the tails of signals received in samples where the waves repeatedly sample the same region in space due to multiple scattering produced by inhomogeneities in the material; Variables u -the wave field displacement; t c -the center of the time window of length 2T; A tr and t trthe amplitude and the travel time of the wave following trajectory t r ; τ trthe change in the travel time of the wave that propagates along path t r ; f c , Δfthe center frequency and the bandwidth of the coda wave spectrum; δυthe change in the sound velocity υ; l* -the transport mean free path; Kdthe coefficient of decorrelation which indicates the level of decorrelation between two coda waves; εthe stretching coefficient and maximum of dilatation parameter  which corresponding to the relative variation in the sound velocity υ; N, Tthe normal and the tangential forces; σ, τthe normal and the tangential stresses; 29 a, bthe normal and the tangential displacements;  -the radial coordinate; cthe radius of contact zone;

Fig. 1 .

 1 Fig. 1. The Lennard-Jones 6-12 potential approximates the intermolecular interactions of two atoms.A more general form of the nonlinear stress-strain relationship corresponding to the material nonlinearity is easy to obtain [Lan-93] in the following way. The expression for the elastic energy should not depend on the choice of the reference frame, therefore it must contain

  where H(ε) is the Heaviside step function (see [Sol-02], for the clapping or bimodal model, or similar relationships by [Naz-89], and [Ost-91]).

Fig. 2 .

 2 Fig. 2. (a) Hysteretic element characterized by critical values x c and x o . (b) The Preisach space (x o ,x c ) containing three areas: triangles with open and with closed elements, and rectangle where the state of elements is determined by the history of the system. (c) Typical configuration of the history-dependent part of the Preisach space.

  of interrogating and inspecting large areas with a fixed number of permanently attached transducer (see [Jan-90], [Gao-05], [Hay-06]). RAPID stands for Reconstruction Algorithm for Probabilistic Inspection of Damage and it is generally a baseline-dependent technique based on a probability distribution [Gao-05], [She-12],

Fig. 3 .

 3 Fig. 3. (a) Schematic of typical sparse array geometry. (b) Typical result from the SSM baseline-free nonlinear RAPID GWT [Het-16b].

Fig. 4 .

 4 Fig. 4. Experimental setup of the resonant scanning laser vibrometry.

Fig. 5

 5 Fig.5illustrates the result obtained with the RESLV technique in a higher harmonic mode.

Fig. 5 .

 5 Fig. 5. Impact in CFRP (cutout window Airbus 350, left) and its RESLV images (right) [Sol-14].

Frequency

  Fig. 6: RETS principle scheme.

Fig. 6 Fig. 7 :

 67 Fig.6presents a scheme of the RETS technique. The thermal defect response at the LDR frequency is measured and visualized with an infrared camera. To enhance the signal-to-noise ratio of the nonlinear imaging the signals operate in lock-in mode. To this end, the amplitude of ultrasonic excitation at the LDR frequency is modulated at the lock-in frequency. A temperature image sequence of the surface is then recorded with the infrared camera and a discrete Fourier transformation at the lock-in-frequency is applied to compress this image sequence into a pair of amplitude and phase images.

Fig. 8 .

 8 Fig. 8. Scheme of resonant shearosonics.

Fig. 9 .

 9 Fig. 9. Frequency selectivity of the LDR based imaging techniques: RETS (left), RESH (center) and RESLV (right) images of three circular flat-bottom holes in a specimen: the frequencies of maximum output are indicated in the images [Sol-14].

Fig. 1 .

 1 Fig. 1. Schematic representation of coda wave generation due to: (a) multiple diffusion (b) multiple reflections.

Fig. 2 :

 2 Fig. 2: Typical modifications of coda signals: time delays induced by sound velocity changes linked to material parameters perturbations (a) and distortions of the signal shape (b).

  been developed to analyze CWI signals. We will present now the two most used of them: the use of time-shifted correlation coefficient, also called the doublet technique, and the stretching technique.



   where l e is the elastic mean free path (i.e. the mean distance the wave propagates before being scattered by a scatter [Sch-97]), Snieder demonstrated that medium perturbations act as a propagation time change: tr is the change in the travel time of the wave that propagates along path t r . In Eq. (3) any shape modification in individual signals A tr S(t) are ignored which in particularity means that effect of dispersion and attenuation are neglected. Thus the validity of Eq. (3) is related to a quality factor of the material; high quality factors correspond to negligible attenuation and, 52 due to the Kramers-Kronig relations, to negligible dispersion.

  the autocorrelation of the source signal.

  average for the wave paths with arrivals in the time interval (t c -T, t c + T). The correlation coefficient attains its maximum at: average perturbation of path t r . When the perturbations are equal for all the paths, s tr tt  . The maximal value of the time-windowed cross-correlation coefficient is given by:

  ) for the relative velocity perturbation   /   determination is generally called the seismic doublet technique [Pou- 84]. However, in the case / const



  [Sni-06] has shown that the mean perturbation of the travel time vanishes and that the variance of the path length is



  is the dilatation parameter corresponding to the relative variation in the sound velocity υ (see [Lob-03],



  i . As proposed by Larose and Hall [Lar-09], we use spline functions to realize the interpolation. Three examples of comparison between a measured coda wave (green line) and interpolated reference coda waves with various relative velocity changes α (blue dotted line), (a) α = -0.002, (b) α = 0.002, and (c) α = 0.0065, are displayed in Fig. 3. Here  max = 0.002.We search  max in a limited range of values [ 0 ,  1 ]. For a given precision, the speed of the search algorithm can be improved with the use of an iterative technique: We first calculate  max with a step  1 between  0 and  1 . Then we reiterate the calculation in a smaller range of values around [' 0 , ' 1 ] with a smaller step  2, and so on until reaching the desired precision . However, it is necessary to verify that the obtained maximum is an absolute maximum. As shown by, we can verify that the evolution of 56 as a function of  has a sinus cardinal shape, as shown in Fig.4. Moreover, we need to verify, as explained by, that the obtained  max is realistic by comparing it to RMS  given by:

Fig. 3 .

 3 Fig. 3. Principle of the stretching method: calculation of the correlation coefficient , () c S tT R  between the measured coda wave (green line) and the interpolated reference coda waves with various relative velocity changes α (blue dotted line) to determinate the value  max that maximizes , () c S tT R

Fig. 4 :

 4 Fig. 4: Evolution of the correlation coefficient (Eq. (21)) as a function of the relative velocity changes α. Points P1,P2, and P3 correspond to the three values of α for which time domain signals are displayed in Fig. 3.



  obtained for a given time window used for the stretching parameter measurement [t c -T, t c +T]. If the obtained max  

  cc I t f around a time t c and a frequency f c by integrating the square of the mean value of the obtained spectrum in a frequency range Δf around f c [Wea-95]. The precision of this estimation is linked to the RMS value of the fluctuations ofˆ( , )

Fig. 5 .

 5 Fig. 5. Procedure of the signal to noise ratio estimation.

  Reversal (TR) [Fin-89], [Fin-92], [Fin-96] is a well-known technique which has been developed in different fields including medical therapy, diagnostic, and

  defects [Bou-06], [Fin-00]. The size of the focal spot depends on source size and form, and on the frequency of the signal emitted. Concerning NDT applications, TR processes have been applied in several classical ultrasonic inspection methods: C-scan with immersed samples [Cha-95], Rayleigh and Lamb wave propagation in plates and hollow cylinders [Ing-96], [Ing-98], [Pra-98], and structural health monitoring [Wan-03], [Soh-07]. In these studies, it was shown that the TR principle improves the detection of flaws in heterogeneous materials for which the microstructure displays a strong speckle noise that is obstructing the observation of a defect echo in classical ultrasonic inspection. On the other hand, researchers have encountered a serious limitation of the traditional TR technique in the fact that only the strongest scatterer can be imaged. The application of the so-called Décomposition de l'Opérateur de Retournement Temporel method [Pra-98], [Pra-02] and successive TR iterations [Wu-92], [Mon-04] may overcome this difficulty to some extent and may enhance the detection by focusing selectively on weaker scatterers. Using these advanced analysis and signal processing techniques, flaws with sizes even smaller than the wavelength can be detected in highly heterogeneous materials such as titanium alloys [Pra-03], [Bor-03]. Time Reversal Acoustics (TRA) provides the ability to focus ultrasonic waves in time and space, regardless of the position of the initial source and of the heterogeneity of the medium in which the wave propagates [Fin-89], [Fin-92], [Fin-96]. This technique has attracted great scientific and technological interests in different fields including medical therapy, diagnostic, nondestructive testing, and underwater acoustics [Fin-00]. In a standard TRA experiment, waves generated by an acoustic source are firstly measured by an array of reversible piezoelectric transducers located around the source, and then time reversed and reemitted by the same transducers array. To improve the quality of focusing, the transducers should cover a closed surface around the medium in order to obtain the wave front information coming from all directions [Cas-92]. However, this is difficult to implement in practice, and time reversal operation is usually performed on a limited angular area, thus adversely affecting the reversal and focusing quality.

Fig. 6 .

 6 Fig. 6. Principle of one-channel TR focalization in a reverberant solid sample.

A

  Fig. 7. Chaotic cavity transducer principle.

Fig. 9 .

 9 Fig. 9. LabVIEW interface developed for one-channel TR focalization.

Fig. 10 .

 10 Fig. 10. Zoom on the LabVIEW interface developed for one-channel time reversal focusing that shows different kinds of source signals available.

Fig. 11 .

 11 Fig. 11. Zoom on the LabVIEW interface developed for one-channel time reversal focusing that shows different kinds of signal processing techniques available.

  [Tan-01] was shown to improve the focusing quality. Indeed, if linearity and spatial reciprocity assumptions are valid in the medium, the time reversal process corresponds to a spatial and temporal matched filter of the propagation. That is to say, the time reversal process maximizes the received output amplitude signal at a given location and a given time, and corresponds to the signal amplitude received at the focus at a given time, for a given input energy. The classical inverse filter is based on the inversion of the propagation operator relating an array of transducers to a set of control points. This technique allows calculation, both in space and time, of the set of temporal signals to be emitted in order to optimally focus on a chosen control point. This broadband inversion process takes advantage of the singular value decomposition of the propagation operator in the Fourier domain.

  inverse filter technique developed by Quieffin [Qui-04] instead of the one proposed by Zverev [Zve-04].

Fig. 12 .

 12 Fig. 12. Principle of the CWI nonlinear defects detection using a pump wave focused at a given region of the sample with the help of the chaotic cavity transducer.

  1. Coda wave interferometry componentThe CWI part consists of an emitting transducer and a receiving transducer. An AWG 2021 arbitrary waveform generator generates a probe signal. In all experiments a short duration frequency swept signal is used. The probe signal is then amplified (by Amplifier Research type 150A250 or 100A400) and emitted by a Panametrics V103-RM contact transducer or a 10.5 mm diameter KEPO ceramic due to their simplicity and low cost. The transducers are glued on the sample. The coda wave signal is then captured by a piezoelectric transducer that is identical to the one used for emission. The received signal is filtered using a 7th order Chebyshev analog high-pass filter to exclude the low frequency pump components and amplified. The signal is acquired by a LeCroy 64 Xi Oscilloscope, and all the Coda wave interferometry measurements are carried out on a computer.

  Fig. 13. Schematic of the experimental setup for CWI nonlinear defects imaging.

Fig. 14 .

 14 Fig. 14. Temporal sequence of pump and probe signals used in the CWI technique.

  . The probe signal has a frequency range of 0.3 -1.2 MHz and a 5 μs duration. The pump waves are excited in two frequency ranges of 35 -150 kHz and 120 -220 kHz. We used a linear sweep in the indicated frequency range combined with the sweep compression technique for TR.

Fig. 15 .

 15 Fig. 15. Example of coda wave interferometry measurement on a glass plate with a localized defect. (b) Measured signal, and (c) its spectrum.

Fig. 16 .

 16 Fig. 16. Evolution of the energy contrast (ratio between the energy at the focal point and the one at 4 cm away from the focal point) as a function of the delay used in the TR process. The delay corresponds to the length of the time window that is processed (time reversed or inverse filtered). Here, we use a 900 μs sweep in the frequency range 50-250 kHz to excite the transducers. Two signal processing methods were used: Quieffin inverse filter (red line), and time reversal (black line).

Fig. 17 .

 17 Fig. 17. Focal spot measured on a glass plate with a VP-0.5 pinducer for three different frequency ranges. Focusing is obtained with a Quieffin inverse filter.

  Fig.17. A more uniform amplitude of energy is obtained for the commercial transducer but with a maximum amplitude four times lower, and thus with a 16 times lower efficiency in creating variation of the time delay (ε) and the envelope (Kd) of the coda wave induced by nonlinear interactions (that is generally quadratic with the amplitude of the pump wave).

Fig. 18 .

 18 Fig. 18. Comparison of the measured energy as a function of frequency at the focal point in a reverberant glass plate for different sources used in the TR process: 7 mm thickness PZT 27 ceramic (black line), 2 mm thickness PZT 27 ceramics (red line), 3 mm thickness PZT 27 ceramics (blue line), the chaotic cavity transducer shown in Fig. 4(a) (violet line), and a Panametrics V103-RM contact transducer (green line).

Fig. 19 .

 19 Fig. 19. Thermally shocked glass plate and measurements on it.

  Fig. 20, neither ε nor Kd show any significant pump-amplitude dependence for the intact specimen throughout the entire test.

Fig. 21 .

 21 Fig. 21. Evolution of the measured (a) stretching parameters ε and (b) Kd on the highly cracked glass sample vs probe delay showing slow dynamic effects after the end of the pumping (vertical line at 800µs).

Fig. 22 .Fig. 23 .

 2223 Fig. 22. Kd vs T delay for different pump duration: from half period of excitation (6.5µs) to 800µs.

Fig. 24 .

 24 Fig. 24. Evolutions of (a) the stretching parameter ε and (b) the decorrelation coefficient Kd vs T trig with T delay = 1 ms and sweep periods -300 µs.

  (mode 1) and pulsed mode (mode 2) of the evolutions of stretching parameter ε and decorrelation coefficient Kd are shown in Fig. 26. As expected, variations of stretching parameters ε and Kd are stronger in continuous mode: ε changes by 0.025% and 0.1%, and Kd changes by 0.07% and 0.12% for pulsed and continuous mode respectively. After switching off the pump amplitude A pump , ε and Kd demonstrate different behaviors in pulsed and continuous mode.

Fig. 26 :

 26 Fig. 26: (a) Measurement protocol used during the test of the glass sample, containing a crack due to a bullet hit from a pneumatic gun, with nonlinear CWI. Results of CWI analysis in continuous mode (mode 1) and pulsed mode (mode 2) of the evolution of (b) stretching parameter ε and (c) decorrelation coefficient Kd. A comparison is given for an intact glass sample.

Fig. 27 .

 27 Fig. 27. Real time stretching analysis. Evolution of (a, c) stretching parameter ε and (b, d) decorrelation coefficient Kd for the intact (a, b) and damaged (c, d) glass plates.

Fig. 28 :

 28 Fig. 28: (a) Principle of defect localization using CWI technique combined with a TR process.Evolution of (b) the stretching parameter ε and of (c) Kd when the pump wave is not focused (green), focused at crack (black) and outside of it (blue).

  Ultrasonic Imaging (NACU) uses the radiation properties for nonlinear imaging of defects. In this technique, a wide-band (focused if possible to improve the imaging resolution) air-coupled ultrasonic transducer detects the acoustic field over a sample vibrated by a contact transducer or a PZT ceramics. The Fourier analysis of the output signal provides information for computer imaging of the nonlinear frequency components distribution over a specimen. This technique was first proposed and developed by I. Solodov [Sol-07]. To extract the nonlinear frequency components generated at the defect, three different customary post-processing techniques were implemented: harmonic and intermodulation filtering, pulse inversion [Sim-96] and the scaling subtraction method. The first two techniques are well-known while the Scaling Subtraction Method (SSM) is more recent (presented by Scalerandi et al. in 2008 [Sca-08a, Sca-08b]). The idea behind SSM consists in the fact that only a high amplitude excitation can activate nonlinear defects in a sample, whereas a low amplitude excitation generates a linear response. Therefore, the low amplitude response can serve as reference signal as if the medium were fully linear. The difference between the responses calculated after multiplication of the linear one by a proper amplitude factor produces a nonlinear indicator i.e. signature of damage. Strictly speaking, the

Fig. 1 .Fig. 1 .

 11 Fig. 1. Schematic of the experimental setup for nonlinear defects imaging via the NACU technique.

Fig. 3 .

 3 Fig. 3. Scale subtraction NACU image obtained on the CFRP laminate with a delamination of 35x35mm 2 at half thickness. The chosen frequency range is 120-220 kHz. The dotted rectangle shows the delamination position.

Fig. 4 .

 4 Fig. 4. Comparison of the images obtained on the CFRP laminate with a delamination of 35x35mm 2 at half thickness obtained by (a) linear through transmission air-coupled ultrasonic measurements (normalized image of the attenuation) and by (b) Scale Subtraction NACU measurements.

Fig. 5 .

 5 Fig. 5. Scale subtraction NACU image obtained on the CFRP laminate with a delamination of 35x35mm 2 at half thickness. The chosen frequency range is 148-152 kHz. The dotted rectangle shows the delamination position.

Fig. 6 .

 6 Fig. 6. Scale subtraction NACU image obtained on the CFRP laminate with a delamination of 20x20mm 2 at 1/4th thickness. The air-coupled ultrasonic transducer scans the plate from the surface closest to the in-depth delamination position. The chosen frequency range is 70-100 kHz. The dotted rectangle shows the delamination position.

Fig. 7 .

 7 Fig. 7. Linear air coupled ultrasonic image obtained on the CFRP laminate with a delamination of 20x20mm 2 at 1/4th thickness. The air-coupled ultrasonic transducer scans the plate from the surface closest to the in-depth delamination position. The chosen frequency range is 70-100 kHz. The dotted rectangle shows the delamination position.

Fig. 8 .

 8 Fig. 8. Linear air coupled ultrasonic image obtained on the CFRP laminate with a delamination of 20x20mm 2 at 1/4th thickness, when the measurements are made on the side of the plate farthest from the in-depth delamination position. The chosen frequency range is 70-100 kHz. The dotted rectangle shows the delamination position.

Fig. 9 .

 9 Fig. 9. Linear air coupled ultrasonic image obtained on the CFRP laminate with a delamination of 20x20mm 2 at 1/4th thickness, when the measurements are made on the side of the plate farthest from the in-depth delamination position. The chosen frequency range is 100-200 kHz. The dotted rectangle shows the delamination position.

Fig. 10 .

 10 Fig. 10. Linear through transmission air coupled ultrasonic image obtained on a 4 mm thickness GFRP plate containing an artificial delamination with a 6 mm radius.

Fig. 11 .

 11 Fig. 11. Linear air coupled ultrasonic image obtained on a 4 mm thickness GFRP plate containing an artificial delamination with a 6 mm radius for different chosen frequency ranges: (a)[52][53][54][55][56][57][58][59][60]

Fig. 13 .

 13 Fig. 13. Scale Subtraction NACU images on a 4 mm thickness GFRP plate containing an artificial delamination with a 6 mm radius. The chosen frequency ranges are (a) 52-60 kHz, (b) 52-55kHz, and (c) 106-110kHz. The area with the delamination is highlighted.

Fig. 14 .

 14 Fig. 14. Subtraction of two air coupled ultrasonic images made with and without a low frequency pumping in the 52 -60 kHz frequency range on a 4 mm thickness GFRP plate containing an artificial delamination with a 6 mm radius. The chosen frequency ranges are 300-700kHz. The area with the delamination is highlighted.

  history of the contact problem started in 1880s when H. Hertz (see [Lan-93]) published the classical solution for two elastic spheres compressed by a normal force. In absence of adhesion and plasticity, this solution is fully reversible. However, the addition of a tangential force and friction ([Cat-38], [Min-49]) makes the problem hysteretic and memory-dependent. It was noted that even a small tangential force acting on two pre-compressed spheres results in appearance of a slip annulus at the periphery of the contact circle where the surfaces are compressed weakly. The coexistence of the stick (central) and slip (peripheral) zones actually means mixed-type boundary conditions that correspond to zero local tangential displacement in the central region and, in the slip annulus, to the Coulomb friction law written for local tangential () and normal () stresses, = (here  is friction coefficient considered as a constant for two contacting materials). The increase in the tangential force results in the slip propagation towards the contact centre. If now the tangential force starts decreasing, a new slip annulus develops at the contact periphery in which =-. Hence, the same values of the normal (N) and tangential (T) forces can correspond to different distributions of stresses and displacements in the contact zone. This fact explains a complex hysteretic character of the solution. Allowing the normal force to evolve [Min-53] adds a new complexity factor to the problem.The matter is that the slip zone always arises at the contact border where =0 and propagates inward, but, if the normal force increases, the contact border itself propagates outward. The result depends on the value of the derivative dN/dT.Finally, the introduction of a general contact geometry instead of spherical profiles introduces even more complexity into the description. Indeed, for two rough surfaces, the contact zone consists of a multitude of contact spots having random geometry. For changing normal force, those contact spots can merge or split. Further, each of them supports slip and stick zones, and the traction distribution in the stick zones can contain residual stresses left from the previous 103 moments of evolution. With continuously varying normal and tangential loading, this complicated picture continuously evolves.

  elastic friction principle (REFP, see [Jäg-95], [Jäg-97], [Jäg-03], [Cia-98a], [Cia-98b] that, under some restrictions, makes it possible to replace an arbitrary contact geometry by a pair of axisymmetric profiles.

Fig. 1 .

 1 Fig. 1. (a) Forces and displacements in a general contact system (case of rough surfaces is shown). (b) Reduced elastic friction principle for axisymmetric bodies. (c) Equivalent axisymmetric system that has the same normal reaction as an original one (e.g. contact of rough surfaces).

  considered contact types are discussed in [Jäg-95], [Jäg-97], [Jäg-03], [Cia-98a], [Cia-98b].

  and  1,2 are Poisson's ratios of the materials; see, for instance, [Cia98b]). In the latter paper,  1,2 =0 is required in addition. The analysis for dissimilar bodies is much more complicated, since partial slip then may occur even for purely normal loading. However (see [Jäg-95], [Jäg-03]), the error produced by neglecting dissimilarity effects, in many practical cases, is within the variation error of the friction coefficient in Coulomb's law. A more detailed discussion, including quantification of the dissimilarity effects, is presented in [Mun-92] and [Mun-94].

  (a). As required by the Coulomb friction law, the shear stress (traction) within the slip annulus s<<c equals the normal stress times . The corresponding memory diagram is a function D() introduced (see Fig. 3 (c)) in such a way that the slip annulus s<<c in which ()=(c,) corresponds to the segment Q<<N with Q=N(s) on which D()=1. More generally, any memory diagram ending by a segment on which D()=1 indicates the presence of slip in the corresponding annulus of the contact circle. Note that the propagation of slip from the contact boundary =c to a value =s not only creates traction in the slip annulus but also in the stick zone (<s) as well, whereas in the memory diagram function D() on the corresponding

Fig. 3 .

 3 Fig. 3. Basic traction distributions for a contact between two spheres. (a): Constant N and subsequently constant T forces are applied (a) so that slip progresses from the contact boundary c to a certain value s. In the slip zone s<<c the Coulomb friction condition ()=(c,) is fulfilled. (b) The forces N and T change simultaneously so that the condition Eq. (4) is satisfied, so that slip is absent, and |()|<(c,) everywhere in the contact zone. In both cases, the forces have been applied to a traction-free contact (absence of loading history, initial curve). Groups (c) and (b): Memory diagrams corresponding to the basic solutions: (c) corresponds to (a), (d) corresponds to (b).

Fig. 5 .

 5 Fig. 5. Illustrations of the evolution algorithm. (a) Original memory diagram (gray), N increases from N p to N n . (b) Case YY i.e.   1/ / 1, 0 T N N       ; The force increment T can be equilibrated by the proper choice of D() on the new interval N p <<N n . (c) Case YN i.e.   1/ / 1, 0 T N N       ; The increment T is too large to be equilibrated by setting D()=1 only on the new interval N p <<N n ; Therefore, the slip zone with D()=1 penetrates closer to the contact center. (d) Original memory diagram (gray), N decreases from N p to N n releasing force T 1 . (e) The residual force increment T 2 =T-T 1 is equilibrated by a slip propagation. In sets (b), (c), and (e), the thick black lines indicate newly generated parts in the memory diagrams.

  Fig. 5 (b) for the new one, D n ()).Correspondingly, Eq. (21) defines the updated value D n (N n ).

  Case YN. Now we consider the situation when the normal force increases, N>0, but  . Then, in general, the balance of forces in the system (Eq. 1) cannot be achieved by a simple definition of the memory function on the new segment N p <<N n only. slip on the new segment (see rule (III)), this would be only enough to fit T in a degenerated case   1/ / 1 TN     . In a more general case   1/ / 1 TN     , the force increment T can not be equilibrated via a proper selection of the memory function on a new interval only. Defining the memory function according to Eq. (22) will only correspond to a partial force increment T 1 that equals for the moment. In order to take into account the unbalanced part, it is necessary to also include the old interval 0<<N p . The process to properly do this is fixed by the rule (III) which allows us to modify the memory function at any -value starting from the contact boundary value =N n (remember that  has the sense of a generalized radial coordinate) and to use the slip solution expressed in Eqs. (9)-(11). The situation is illustrated in Fig.5 (c) where the case of positive slip with D=+1 is shown for definiteness, while in a general case D should be set to sgn(T 2 )=sgn(T). Firstly, the memory function D() has been set to sgn(T) on the interval N p <<N n , thereby equilibrating the partial force increment T 1 =N. Subsequently, the remaining part T 2 =T-T 1 is compensated for by setting D n ()=1 in the contact zone adjacent to the new segment N p <<N n . In other words, we have to shift point A in Fig.5(c) starting from position A' and setting the memory function to sgn(T) on the interval AA' until the desired value T 2 / is reached for which:

  of the slip zone towards the contact zone can happen via many scenarios. For instance, in the scenario shown in Fig.5(c), the zone <N n ends with a segment where D() is equal to 1, to that setting D() to 1 in this segment does not compensate anything, and we have to address a zone closer to =0. On the other hand, when the memory diagram is such that the old contact area <N n ends with a value of D() different from 1, the update produces an immediate contribution. Furthermore, it is entirely possible that the developing slip zone erases a certain part in the structure of the memory function by setting it to 1. This would happen in Fig.5 (c) if T were large so that point Q would have progressed further to the left thereby completely eliminating the segment with D()=-1 and eventually also other segments. Point Q can even reach =0 resulting in a gross sliding behavior. The strength of the proposed approach is that it automatically complies with the force balance without specifying or classifying all structural changes that can potentially occur to the memory diagram.

)

  At the same time, the memory diagrams D() and D() are related to the traction distribution and local tangential displacement distribution which correspond to the same deformation state of the body. Hence, the two diagrams are structurally identical, and can be obtained by different scale transformations of the original -coordinate. This yields the following relationships:

Fig. 6 .

 6 Fig. 6. Memory diagram (a) and corresponding traction distribution (b) for contact of two identical spheres. The scale transformation is given by   1/ 3 KR   , equation immediately following from the Hertz theory.

  (a)). A typical behavior is shown in Fig. 7(c); a straight-line segment Q<<N on which D()=1 corresponds to slip. Slip propagates inward erasing the previously saved curvilinear segment. At the moment t=0.5 the normal force starts decreasing, therefore the total length of the 0 (future evolution of the memory diagrams is shown with arrows in all sets (b)-(e)), and the system runs in the regime N (Fig. 7(d)). Then, point Q shifts closer to the left end of the memory diagram =0 (Fig. 7(e)); at the moment t=1 Q reaches 0 which means that partial slip transforms into total sliding when T= N.Two more general cases of the MMD application are presented in Fig.8. Fig.8(a) illustrates the particular situation when the two force arguments, N and T, are linked by a functional relationship, N=N(T), so that actually there is only one independent argument but it evolves in a non-trivial manner. The resulting curve b(T) shown in Fig.8(a) is typical for one-parametric hysteresis response: it exhibits closed loops for periodic T(t), partial increase in the argument T on a globally decreasing branch results in the creation of an inner loop, and each completed inner loop has the property of end-point memory, where the curve exits the outer loop with the same tangent as just before entering it, etc.

Fig. 8 .

 8 Fig. 8. Tangential force-displacement curves for some particular loading histories as shown in the insets. (a) One-parametric hysteresis, in which N and T are linked via a functional dependence. (b) Two-parametric hysteresis in which N and T are independent. In all plots, N 0 is a characteristic value of the normal force; a 0 =a(N 0 ).

  of the MMD is an extremely cumbersome task. The formulation in terms of MMD strongly reduces the complexity.

Fig. 1 .

 1 Fig. 1. Data exchange between the crack model and the solid mechanics unit in the case whenthe crack model is driven by displacement.

  Introduction of surface roughness considerably changes the contact behavior since asperities can recede under load. Application of a normal compressive load N>0 results in the appearance of positive a which was not possible in the case of flat surfaces (Fig. 2). Modeling of force-displacement relationship corresponding to normal loading of contacts with rough surfaces is an extensive research topic since 1950s-1960s when pioneering works [Arc-57], [Gre-66] were published. Contact of rough surfaces can be described using three different parameters related to contact area: nominal contact area A n i.e. cross section of the mesoscopic cell by the crack plane, real (atomic) contact area A, and geometric (truncated) contact area A g obtained in a virtual situation when rough profiles can freely penetrate into each other. Following resent experimental estimations [Hyu-07] and theoretical (classical [Bow-39] and recent [Pag-10], [Pag-14], [Per-02]) results, we accept the proportionality between the real contact area A and the applied normal force n AN, equal to the root mean square of the surface slope, and 2   . Further, it is easy to express the geometric area A g obtained when the rough profiles that penetrate into each other at depth a as 127 2 if we assume that the deformed asperities are approximately spherical, and the deformation is elastic (Hertz result, [Lan-93]). Combining Eqs. (1-3) produces the following expression:

  practice, this refers to highly conforming surfaces including flat faces described in the previous sections. Another example is the aperture distribution for rock joints measured by Xray methods (see Fig.11in [Sha-08]). Indeed, the faces of a joint could initially belong to an intact piece of material so that the fault surfaces are globally conforming whereas low-scale roughness is invisible for the method having resolution of order of 1 mm. Secondly, the horizontal tangent refers to an essentially open crack in which points in atomics contact practically do not appear. Finally, an intermediate tangent inclination coefficient k (0<k<∞) same dependency N(a)~a 2 has been measured in [Biw-04] for two aluminum plates in contact; the same result has been obtained by an indirect estimation using nonlinear ultrasound data [Dri-96]. The considered experimental situation suggests that the quadratic dependency is a possible approximation for two globally plane surfaces with uncorrelated

  completes the description of the crack model. The result is the possibility to calculate contact forces N and T per unit area for any normal (a) and tangential (b) displacements. This is the principal difference of the proposed approach in comparison to the simple flat crack model depicted in Fig. 2. The introduction of roughness on internal contacting surfaces and the account for partial slip finally allowed us to organize calculations in the explicit manner as shown in Fig. 1. This advance is related to the fact that the account for roughness or an equivalent axisymmetric contact shape produces the Coulomb sliding

A

  successful model for elastic wave propagation in materials containing cracks requires two components: a crack model and a solid mechanics unit. The crack model should take into account the microscopic normal and tangential contact behavior at the crack interface, whereas the solid mechanics unit is needed for solving the elastic wave equations. In this study, the elasticity equations will be solved using the Structural Mechanics Module [COM-15a] of the commercially available, finite element based, software package COMSOL Multiphysics. This specific module was especially developed for the analysis of mechanical structures that are subject to static or dynamic loads and is therefore extremely suitable for the modeling of elastic wave propagation problems. Exploiting the dedicated physics interfaces and tools of the software package, COMSOL also allows incorporation of particular crack models.When using COMSOL for the modeling of wave propagation problems involving contacts, with or without friction, several approaches can be considered. A first method is based on the use of contact pairs [COM-15a]. A contact pair in COMSOL consists of two sets of boundaries, which are called the source and destination boundaries. In the Structural Mechanics Module of COMSOL, the contact pair environment introduces extra boundary conditions at the internal crack faces, preventing the destination boundaries to penetrate the source boundaries.

  (b) The calculated displacement values at all Gauss points on the contact interface are used as an input in the displacement-driven crack model in MATLAB, in order to determine corresponding stress values at these positions, which on their turn are again outputted to COMSOL to redefine the boundary conditions at the contact interface. (c) Steps (a) and (b) are repeated for the next time step, until the desired calculation time is reached. Using this approach, the final model thus contains two components: the constitutive displacement-driven crack model implemented in MATLAB and the wave propagation model implemented in the Structural Mechanics Module of COMSOL.

Fig. 5 :

 5 Fig. 5: Calculated relative displacements u n and u t (top) and contact stresses  and  (bottom) at the central point on the crack interface in case of a shear wave excitation at 100 kHz with amplitude A = 100 nm. The crack was implemented by means of a linear model where no forces are imposed on the crack faces.

Fig. 6 .

 6 Fig. 6. Calculated relative displacements u n and u t (top) and contact stresses  and  (bottom) at the central point on the crack interface in case of a shear wave excitation at 100 kHz with amplitude A = 100 nm. The crack was implemented by means of a nonlinear model where only the normal force N is imposed on the crack faces.

Fig. 7 .

 7 Fig. 7. Calculated relative displacements u n and u t (top) and contact stresses  and  (bottom) at the central point on the crack interface in case of a shear wave excitation at 100 kHz with amplitude A = 100 nm. The crack was implemented by means of the full crack model with both normal force N and tangential force T imposed on the crack faces.

Fig. 8 .

 8 Fig. 8. Calculated normal and tangential stress-displacement curves at the central point on the crack interface when using the full crack model in case of a shear wave excitation at 100 kHz with amplitude A = 100 nm. The top figure illustrates the quadratic normal behavior in case of contact (i.e. u n ≤ 0), according to Eq. (8). The bottom figure illustrates the tangential hysteresis behavior with several hysteresis loops. The tangential reaction curve also clearly indicates the switching between different defect states (contact loss, total sliding and partial slip). For clarity, the reaction curves are only plotted for values u n smaller than 5 nm in the normal case and values u t larger than -20 nm for the tangential case, since for other displacement values both reaction curves are constantly zero.

Fig. 9 .

 9 Fig. 9. Fourier spectrum of the calculated relative normal and tangential displacements (u n and u t ) obtained for both the linear model and the full crack model at the central point on the crack interface in case of a shear wave excitation at 100 kHz with amplitude A = 100 nm.When using the full crack model, harmonic frequencies are generated.

Fig. 10 .

 10 Fig. 10. Fourier spectrum of the relative normal and tangential displacements signals (u n and u t ) at the central point on the crack interface in case of a shear wave excitation at 100 kHz, after applying the Scaling Subtraction Method (SSM). The displacement signals were thus obtained by subtracting the displacements obtained at an amplitude A low = 10 nm from the signals obtained at amplitude A = A high = 10 A low = 100 nm. Harmonic frequencies are clearly generated.

  high = 10 A low = 100 nm. The generated harmonic frequencies in the resulting SSM signals for both the relative normal and tangential displacements are now clearly visible. The harmonics present in the relative normal displacement signal are mainly caused by clapping effects (i.e. normal interaction between the crack faces), whereas the harmonic frequencies present in the relative tangential displacement signal are mainly generated by friction effects (i.e. tangential interaction between the crack faces). Note that the amplitude of the clappinginduced harmonics is larger than that of the friction-induced harmonics. Moreover, due to the oblique incidence of the shear wave, friction effects produce both even and odd harmonics (in contrast to the normal incidence case where only odd harmonics would be generated [Mez-11]). These results were also obtained in the numerical study of crack-wave interactions by.

Fig. 11 .

 11 Fig. 11. Fourier spectrum of the relative normal and tangential displacement signals (u n and u t ) at 21 positions on the crack interface in case of a shear wave excitation at 100 kHz, after applying the Scaling Subtraction Method (SSM). This displacement signals were thus obtained by subtracting the displacements obtained at an amplitude A low = 10 nm from the signals obtained at amplitude A = A high = 10 A low = 100 nm. Harmonic frequencies in the relative normal displacement signals, generated by clapping effects, are mainly occurring in the left part of the crack interface (i.e. negative positions), whereas harmonics in the relative tangential displacement signals, generated by friction effects, are also occurring near the center of the crack interface.

Fig. 12 :

 12 Fig. 12: Vibration patterns of the scaled, subtracted normal (left figure) and tangential (right figure) displacement components, u n and u t , in the full aluminum sample when filtering around the generated second harmonic frequency f = 200 kHz. The figures clearly illustrate the generation of this second harmonic at the crack interface, with clapping induced nonlinearity mainly occurring in the left part of the crack (left figure) and friction induced effects both in the left part and near the center of the crack (right figure).

Fig. 12

 12 Fig.12shows vibration patterns in the full sample when filtering the scaled, subtracted normal and tangential displacements, u n and u t , around the second harmonic frequency (f = 200 kHz). These vibration patterns allow one to verify that nonlinear features are indeed generated at the crack interface. The left figure shows the result for the scaled, subtracted normal displacement signal u n , where the normal displacement was calculated according to the first Eq. (11), while the right figure shows the result for the scaled, subtracted tangential displacement signal u t , with the tangential displacement being calculated in accordance to the second Eq. (11). The left figure can again be linked to clapping induced nonlinearity, while the right figure is linked to friction induced nonlinearity. In both figures, it is easy to verify that nonlinearities (here second harmonic generation) are indeed being generated at the position of the crack, before they propagate further in the sample. Moreover, these figures also verify the result discussed before (Fig.11) where clapping effects were found to be

Pre-35], [Kra-89]). Suppose

  Preisach model is a mathematical formalism suitable for the description of a wide class of hysteretic dependencies in 1D (see [

	we have a collection of
	bistable elements such that each of them can be open (contribution 1) or closed (contribution
	0). For each of them, there are two critical values of the argument x: x o and x c (x o >x c ). By
	definition, the element (x o , x c ) is open if x>x o and closed if x<x c . (see Fig. 2). A hysteretic
	function y is obtained then as a sum of all contributions of all elements (i.e. contributions of
	the open elements since the closed ones do not contribute) weighted by a 2D weighting
	function (x o ,x c ) called the Preisach density.

  Titre en français : Contrôle non destructif par des méthodes d'acoustique non linéaire pour des applications aéronautiques Résumé en français : Ce travail de thèse est une contribution au développement des méthodes d'acoustique non linéaire pour le contrôle non destructif et l'imagerie de défauts dans les solides. Dans ce travail, des modifications sont proposées pour deux méthodes récentes de contrôle non destructif par acoustique non linéaire : l'interférométrie de coda couplée au retournement temporel, et l'imagerie non linéaire par ultrasons aériens. Le principal avantage de la première méthode est sa sensibilité extrême liée а l'accumulation des effets induits par des changements, même faibles, des propriétés de l'échantillon durant la formation de la coda. La deuxième méthode apporte une approche complémentaire en permettant de réaliser un contrôle sans contact. Les techniques développйes ont été testées sur des échantillons présentant des défauts artificiels à des emplacements connus, et leur performance a été étudiée. La deuxième partie de ce travail porte sur la description théorique des non-linéarités acoustiques de contact et leur utilisation pour le développement d'une boite à outils numériques permettant la simulation d'ondes acoustiques dans des structures complexes contenant des contacts internes. Un model physique décrivant le décalage tangentiel de deux corps en contact en présence de friction est proposé. Il aboutit à une solution analytique pour la relation présentant une hystérésis entre les déplacements de contact normal et tangentiel et les chargements. Ce modèle est utilisé comme condition aux frontières pour les surfaces de contact internes (défauts) dans un modèle de propagation d'ondes acoustiques utilisant un logiciel d'éléments finis commercial.Mots-clefs : Acoustique non linéaire, Interférométrie de Coda, Non linéarités de contact Titre en anglais : Nonlinear acoustic nondestructive testing for aeronautical applications

(a) (b) (c)

Plotting now in Fig. 24 ε and Kd as a function of T trig for the following fixed CWI parameters, duration of the sweep signal is 300 µs and T delay =1 ms, we can see that a large time scale effect plays an important role in nonlinear coda mixing. Indeed, as the delay from each pump pulse remains the same, only the frequency of repetition of the pump pulses changes.

Decreases of ε and Kd values for T trig varying between 0 and 200 ms are 10 and 18 dB, respectively.

In this sample, the effect of nonlinear coda waves mixing is strong and easy detectable, that is why it has been used to improve our comprehension and characterization of fundamental mechanisms of crack behavior under ultrasound excitation. The TR resolution highly depends on the number of transducers used for measuring the received signals. We have demonstrated that even one transducer coupled with a chaotic cavity is enough to have some limited spatial selectivity. In other words, the TR procedure works even in poorest conditions, despite formally it can be guaranteed only when the information necessary for focus retrieval is complete. Thus, a perspective of this work is to use several transducers to improve the focusing resolution and/or to implement more precise imaging techniques based on multi-channel CWI acquisition.

Model specifications

The model geometry consists of a rectangular aluminum domain of 50 mm width and 100 mm height, as illustrated in As illustrated in Fig. 4, the full geometry is meshed using quadratic triangular elements. To reach convergence, COMSOL requires approximately 6 second-order mesh elements per 134 wavelength, which for this particular model corresponds to a maximum element size of approximately 5.2 mm. Smaller mesh elements are generated in the region of the crack, since the MMD algorithm requires a small spatial discretization size in order to obtain accurate solutions. Here, a fixed number of 150 mesh elements at the internal crack boundary is selected, corresponding to an element size of approximately 0.13 mm.

The wave propagation problem is solved using the implicit generalized alpha time-dependent solver, which is the preferred solver to be used for structural mechanics problems in COMSOL. For time-dependent wave propagation problems, COMSOL requires to have at least 20 time steps per wave cycle. However, in order to get accurate solutions, the MMD algorithm requires a much smaller time discretization size. The time step chosen for this particular simulation is t = 50 ns, which corresponds to 200 time steps per wave cycle.

Simulated normal and tangential reaction curves

When solving the model, the excited shear wave will propagate in the aluminum sample until it reaches the crack and starts interacting with it. Part of the energy of the wave will be This free motion is illustrated by the fact that both contact stresses  and  are zero (i.e. stress free situation), and both relative displacements u n and u t are symmetrically oscillating around the initial relative displacement. In Fig. 6, relative displacements and contact stresses are calculated in case the normal force N is included in the model. According to Eqs. ( 7) and ( 8), this normal force differs from zero in case of contact (i.e. 0 a  or 0 n u ).

Consequently, the normal contact stress  will also differ from zero in case of contact. This is clearly the case in Fig. 6, where negative normal stress values are observed at those periods in time where the relative normal displacement u n is also negative. Moreover, the negative stress values also indicate the fact that both crack faces are pushing against each other and are therefore in a compressive state. The influence of introducing the normal force N in the model can also be seen in the relative displacements plot, where the relative normal displacement signal is changed when compared to the curve obtained in the linear (stress-free) model (Fig. 5). Because of the presence of compressive stress in case of contact, the relative normal

General summary and conclusions

This study concerns nonlinear ultrasound nondestructive testing of damage in materials such as aeronautic components and is performed in the framework of the ALAMSA European project devoted to self-repairing materials for aeronautics. The original part of the dissertation contains two major components.

The first (experimental) part describes two methods based on nonlinear ultrasound capable of detecting damage in composite aerostructures. One of them is based on a nonlinear version of the coda wave interferometry. Coda waves in materials with high quality factor accumulate a large amount of information on structures in which they travel. This means that two codas measured in almost identical samples with only tiny differences in properties will differ considerably. A nonlinear version of this principle uses a low frequency pump wave excited in the sample and then filtered out in the measured coda. In that case, if the coda is affected by the pumping wave, this can only occur due to nonlinear frequency mixing. In linear consolidated materials, nonlinear effects are indicators of damage. Therefore, by comparing two coda signals obtained with and without pumping, the presence of damage can be identified. In our study, we combine this method with the time reversal principle that allows one to focus acoustical energy at a given spot by inverting acoustic wave fronts. A usual difficulty of these kind of techniques face is associated with real field applications when structures under study have a complex shape, defects are hardly accessible and so on. In order to develop a suitable technique for aeronautic applications, time reversal has been performed in its poorest conditions possible, i.e. when, instead of a certain number of transducers and multi-channel acquisition, only one transducer is used. To do this, the transducer was coupled with a chaotic cavity that partly compensates the lack of acquisition points by multiple reverberations and therefore a larger amount of information can be obtained. We show that a proper combination of these principles provides robust damage detection and discuss associated effects such as slow dynamic processes induced by vibrations of damaged structures.

Another challenge appearing while testing real aerostructures is the need for remote detection.

To deal with this difficulty, the other experimental method based on air coupled ultrasound was elaborated. The excitation is performed by glued transducers but the detection is fully remote. Using this technique, we produced a number of images of artificially damaged composite samples and compared the results with the location of the actual damage.

147

The theoretical part of the thesis can be seen as modeling support for nonlinear nondestructive testing and imaging. The objective was to model waves in materials containing frictional cracks. The concept we elaborated includes the consideration of roughness of crack faces. We demonstrate how contact mechanical approaches previously developed to deal with contact between rough surfaces help the construct nonlinear hysteretic load-displacement relationships at the crack faces. These relationships are considered as boundary conditions posed at internal boundaries in a sample. They are used by a solid mechanics unit programmed via available finite element software (COMSOL). The developed numerical code is capable of calculating all elastic fields in the sample including their nonlinear components.

The perspectives of this study are related to increasing the resolution of the developed experimental methods. In particular, a multi-channel coda wave interferometry with the subsequent multi-channel time reversal can drastically improve the imaging performance. We also plan to compare nonlinear signatures obtained on engineered samples with known damage locations to modeling results.