Analytical and numerical modelling of non-collinear wave mixing at a contact interface

Journal of Sound and Vibration 468 (2020) 115078

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Analytical and numerical modelling of non-collinear wave mixing at a contact interface L.R. Francis Rose a, Philippe Blanloeuil b, *, Martin Veidt c, Chun H. Wang b a

Defence Science and Technology Group, Melbourne, VIC, 3207, Australia School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW, 2052, Australia c School of Mechanical & Mining Engineering, University of Queensland, Brisbane, QLD, 4072, Australia b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 8 August 2019 Received in revised form 5 November 2019 Accepted 6 November 2019 Available online 7 November 2019 Handling Editor: Dr S Felix

An analytical model is presented for wave mixing of two incident non-collinear plane waves at a plane contact interface characterised by a nonlinear traction law that is representative of rough-surface contact, where the normal response is described by a quadratic nonlinearity whereas the tangential response is linear. The approach relies on a decomposition of the scattered field into contributions exhibiting distinctive symmetries with respect to the interface, combined with a perturbation analysis that treats the nonlinear response as a small correction relative to a linear response. The decomposition is shown to lead to uncoupled governing equations, thereby enabling explicit analytical formulae to be derived for the scattered field amplitudes for a linear interface, which in turn provides the source terms for determining the amplitudes and directions of the nonlinear mixed waves. The characteristics of ultrasound wave mixing at an interface are thus shown to differ in several respects from bulk wave mixing due to material nonlinearity. In particular, interface mixing generates mixed waves at the sum and difference frequencies that propagate in both forward and backward directions (transmitted and reflected), whereas bulk mixing only produces a forward propagating mixed wave at the sum frequency. The directions of propagation for interface mixing are determined from the sum or difference of the components of the incident wave vectors parallel to the interface, coupled with the sum or difference of the frequencies, which constitutes a generalisation of Snell's law. These directions generally differ from the sum of the wave vectors that provides the propagation direction for bulk mixing. A computational model is also presented that confirms these predictions of the analytical model. The implications of both models for the design and interpretation of experimental investigations are briefly discussed. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Nonlinear ultrasonic Scattering Analytical solution Finite element

1. Introduction Non-collinear wave mixing provides an attractive approach for detecting various forms of material degradation or structural damage that give rise to a nonlinear ultrasonic response. This approach was first demonstrated by Croxford et al. [1] to be capable of detecting nonlinearity induced monotonic plastic deformation, as well as low-cycle fatigue damage. They

* Corresponding author. E-mail address: [email protected] (P. Blanloeuil). https://doi.org/10.1016/j.jsv.2019.115078 0022-460X/© 2019 Elsevier Ltd. All rights reserved.

2

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

noted several advantages of non-collinear mixing relative to more conventional approaches based on acousto-elasticity [2e4], or higher harmonic generation [5e7]. These include frequency selectivity (the frequency of the mixed wave can be readily distinguished from harmonics of the incident waves if the driving frequencies are chosen to be unequal), modal selectivity (the mixed wave mode differs from that of the input waves), and directional selectivity (the mixed wave propagates in a different direction from the input waves and their higher harmonics). The experimental work of Croxford et al. [1] was guided by the theoretical results first derived by Jones and Kobett [8] for bulk wave mixing within the framework of the five-constant theory of hyperelasticity, where the two constants characterising linear elastic response in an isotropic material are augmented by three third-order elastic constants that are conventionally chosen to be the Murnaghan constants, or the Landau-Lifshitz constants [7]. This theoretical framework involves no characteristic length or time, and accordingly, the wavespeeds are independent of frequency. Consequently, the formula derived by Jones and Kobett [8] for the interaction angle between the input waves only involves the frequency ratio of the input waves (and not the individual values of frequency), as well as the ratio of wavespeeds (which in turn only depends on the Poisson ratio). This theory for bulk wave mixing has recently been extended to the case of Lamb waves [9e12] interacting with material-nonlinearity over a domain, where the wave modes are now dispersive, and consequently the conditions for generating a mixed-frequency wave can no longer be expressed in terms of a frequency ratio. Wave mixing has also been investigated for other forms of distributed damage [13e16]. More recently, Alston et al. [17] have presented a detailed experimental investigation of wave mixing at a contact interface between two aluminium blocks that are pressed together by a variable clamping force. In the absence of a theoretical model for mixing at a contact interface, the theory for bulk-wave mixing due to material nonlinearity [8] was used to guide this experimental investigation of contact interface nonlinearity. Focussing on a symmetrical configuration relative to the normal to the interface where both inputs were shear waves with equal angle of incidence, Alston et al. [17] sought to establish a characteristic signature (or fingerprint) for interface mixing in a two-parameter space spanned by (i) the interaction angle between the input beams, and (ii) the frequency ratio (between the centre frequencies of the tone bursts used as inputs). They found, in particular, that the mixed wave attained its maximum amplitude at an interaction angle of 75
, which is quite different from the angle of 120
expected from Ref. [8] for bulk-wave mixing, but close to the angle of 78
predicted by a finiteelement computational model first developed by Blanloeuil et al. [18] using a unilateral contact law to model the interface. Alston et al. [17] also investigated the case of bulk-wave mixing in a monolithic aluminium block, verifying the theoretical optimum interaction angle of 120
for that case, and also establishing a baseline amplitude for this mixed wave from the background material nonlinearity which was found to be an order of magnitude smaller than that obtained from interface mixing. The aim of the present work is to provide an analytical model for wave mixing at a contact interface that is characterised by a general traction law between the normal and shear stresses at the interface and the relative normal and tangential displacements across the interface. This interface model is conceptually equivalent to having nonlinear springs across the interface, which has been shown to be an effective model for rough-surface contact in the presence of a compressive stress that prevents loss of contact [19e27], as well as for other interface conditions such as for two substrates bonded by a thin adhesive layer [19,28]. The analytical solution relies on a perturbation analysis, i.e. successive approximations, in which the first approximation (or zeroth-order solution) corresponds to linear springs across the interface. A key innovation is the decomposition of the relative displacement across the interface into contributions exhibiting either Mode I or Mode II symmetry, in the terminology of fracture mechanics. This decomposition was introduced for single-wave incidence in our recent work [29], where it was shown to lead to explicit solutions for the amplitudes of the transmitted and reflected waves, including higher harmonics, whereas previous approaches required numerical solutions to determine these amplitudes [24,30]. Zhang et al. [31] have recently presented analytical and numerical modelling of non-collinear shear wave mixing at an imperfect interface. Their analytical model is designed to exclude a linear-springs response by assuming zero values for the interface compliance coefficients, which precludes a perturbation analysis in the sense presented here in Section 4. Thus, the present work provides the first analytical model directly applicable to the experimental investigations in Ref. [17] for roughsurface contact, which involve a linear-springs response, as do the computational models in Ref. [31]. The presentation is organised as follows. First, the behaviour of rough-surface contact is briefly discussed in Section 2, with the aim of selecting a suitable simplified model for the present work. Next, the interface wave mixing problem is formulated in Section 3, where the new scattered-field representation is also presented. Section 4 outlines the procedure and key results of a perturbation analysis, with the detailed derivations being included in an Appendix for completeness and clarity. Section 5 provides a discussion of these analytical results, noting several important differences with bulk-wave mixing due to material nonlinearity. Section 6 presents a finite-element computational study that provides a validation of the analytical results, as well as exploring a wide range of input parameters. Finally, Section 7 provides some concluding remarks.

2. Traction law for rough-surface contact Fig. 1 (a) illustrates schematically the contact interface between two blocks with nominally flat, but rough surfaces, with the scale of the surface roughness being greatly exaggerated for clarity. The contact mechanics can be visualised as involving distributed springs, of nonlinear stiffness, between the blocks, as illustrated in Fig. 1 (b). As the blocks are brought closer together, contact first occurs when the gap separation reaches a value shown as dmax in Fig. 1 (c), where the gap separation d is

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

3

Fig. 1. Rough-surface contact. (a) Schematic of first contact at asperities, with the scale of surface roughness greatly exaggerated for clarity. (b) Contact mechanics modelled by distributed springs between the substrates with a gap opening d. (c) Relation between the (average) interface stress sTot yy and the gap opening d, also showing the variations in stress and relative displacement due to ultrasonic waves. (d) The nonlinear traction law between syy and Duy for rough surface contact (solid red curve), compared with the linear springs approximation (dashed red line) and the unilateral contact law (dot-dashed green lines). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

defined as the distance between two plane surfaces corresponding to the average surfaces drawn through the rough surfaces of the two blocks, as indicated schematically by dashed lines in Fig. 1 (a,b). It is assumed that as the clamping force is increased, the gap separation decreases continuously, as indicated in Fig. 1 (c), which shows schematically the relation between the gap opening d and the average total normal stress sTot yy at the interface (i.e. the overall transmitted force per unit area of interface, rather than the local stress values at contacting asperities). Prior to ultrasonic probing, the blocks are assumed to be subjected to a clamping force that gives rise to a contact stress of magnitude p0

sTot yy ðy ¼ 0; tÞ ¼  p0 ;

(1a)

and a corresponding gap opening d0 . This constitutes the reference state for the interface during probing. Consider now a normally incident longitudinal plane wave specified by iðkyutÞ uinc ; y ðy; tÞ ¼ Ae

(1b)

using standard notation, with A denoting amplitude, u the angular frequency, and k the wavenumber. The total normal stress at the interface is now given by

sTot yy ðy ¼ 0; tÞ ¼ syy ðy ¼ 0; tÞ  p0 ;

(1c)

where syy denotes the normal stress generated by the ultrasonic field. The perturbation method that will be employed in the present work is restricted to situations where the interface stress is always negative, i.e. where




syy ðy ¼ 0; tÞ
< p0 :

(2)

This restriction implies that the interface always remains in contact, i.e. there is no clapping, unlike the case modelled computationally in Ref. [18]. However, the incident wave gives rise to a relative displacement Duy across the interface defined by

4

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Duy ¼ uy ðy ¼ d0 = 2; tÞ  uy ðy ¼  d0 = 2; tÞ:

(3a)

In keeping with conventional practice for interface modelling [19e27], the reference gap opening d0 is considered to be small enough to be ignored when stating boundary conditions across the interface, i.e. Eq. (3a) is considered to be equivalent to

Duy ¼ uy ðy / 0þ ; tÞ  uy ðy / 0 ; tÞ:

(3b)

As can be seen from Fig. 1 (c), the relative displacement due to the ultrasonic waves must be consistent with the restriction that the gap opening remains within the range ð0; dmax Þ; i.e.




Duy
< Min ½d0 ; dmax  d0 ;

(3c)

so that there is no interpenetration across the interface and no loss of contact. The interface response is characterised by a traction law that is assumed to be monotonic and continuously differentiable, so that it can be approximated in the neighbourhood of the reference state ð p0 ; d0 Þ by a polynomial expansion,

syy ðy ¼ 0; tÞ ¼

∞ X

 mþ1 εm K N ; m Duy

(4a)

m¼0

where ε denotes a small non-dimensional parameter that will be required later for the perturbation analysis, and the superscript N indicates coefficients that pertain to the response normal to the interface. To illustrate the approach, it will be sufficient in the present work to retain only a first-order correction to a linear springs response, i.e. to replace Eq. (4a) by







2

syy ðy ¼ 0; tÞ ¼ K N0 Duy þ εK N1 Duy :

(4b)

From the schematic relation shown in Fig. 1 (c), the linear springs coefficient K N 0 can be expected to vary with the contact stress p0 : Biwa et al. [25] noted that the experimental data of Drinkwater et al. [32] for a contact interface between aluminium blocks, and a contact stress in the range of 16e190 MPa, is consistent with a power-law relation m KN 0 ¼ Cp0 ;

(5a)

with C ¼ 6  1010 Pa1=2 m1 ; m ¼ 0:5; and the reference stress being taken as p0 ¼ 10MPa. This relation will be used below for illustrative purposes, using the preceding values for C and m. Furthermore, the first-order coefficient was also claimed to be consistent with a power law

1 2 2m1 KN ; 1 ¼  mC p0 2

(5b)

so that K N 1 is independent of the contact stress for m ¼ 0:5; which will be assumed to be the case for the numerical results presented later. This traction law in Eq. (4b) is shown schematically in Fig. 1 (d), where it is contrasted with the traction law for unilateral contact that is used in the computational model of Blanloeuil et al. [18]. It can be seen that the latter is discontinuous, and therefore not amenable to an analytical solution based on the perturbation method. It is also worth noting that the two models are complementary in the sense that they apply for different ranges of the non-dimensional load ratio







 z ¼ p0
syy
;

(6)



where
syy
denotes the amplitude of the normal stress due to the ultrasonic field: the analytical model presented below only applies for z > 1; where clapping does not occur, whereas the computational model of [18] covers the range 0 < z < 1; where clapping occurs, but does not give a nonlinear response for z > 1; because the surfaces were assumed to be ideally flat so that the interface provides total transmission for. z > 1: A traction law has also been derived for shear wave incidence, based on a micromechanical model of contacting asperities [24]. Symmetry considerations dictate that K T1 ¼ 0, where the superscript T pertains to the tangential stiffness coefficients, and that nonlinearity in the shear response would provide a correction only at the second order of approximation [24,29]. For simplicity, it will be assumed in the present work that the shear response of the interface is linear, i.e.

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

sxy ðy ¼ 0; tÞ ¼ K T0 Dux ; Dux ¼ ux ðy/0þ ; tÞ  ux ðy/0 ; tÞ:

5

(7)

Consequently, the nonlinear response of the interface in the present work is entirely attributable to the quadratic nonlinearity of the tension/compression springs characterised by the first-order coefficient K N 1 in Eq. (4b). It is emphasized that this simplification is not a necessary restriction for the perturbation analysis, and indeed analytical results for singlewave incidence with cubic nonlinearity for the shear response have been presented in our earlier paper [29]. Micromechanical modelling [24] suggests the following relation between the linear coefficients,

K T0 ¼

2ð1  nÞ N K ; 2n 0

(8)

where n denotes the Poisson ratio. This relation is in agreement with experimental measurements [23], where it has been shown to also be applicable for oblique incidence. Thus, for the detailed numerical results presented below, the interface traction law will be characterised by the following numerical values for the coefficients in the traction law, T N 1 1 2 KN 0 ¼ 0:19 MPa nm ; K 0 ¼ 0:15 MPa nm ; K 1 ¼ 0:9 kPa nm

(9)

if not stated otherwise. 3. Problem formulation for interface wave mixing Consider a 2D plane strain configuration where two plane waves, identified by suffixes A and B, are incident on a plane interface located at y ¼ 0 between two half spaces having the same isotropic and linearly elastic material properties. We shall focus on the case where (i) the interface is characterised by the traction law in Eqs. (4b) and (7a) and , i.e. quadratic nonlinearity for tension/compression only, and (ii) both waves are incident from the same side of the interface. The case where both incident waves are shear-vertical (SV) waves has been the most widely investigated previously, both experimentally and computationally [1,17,18,31], and accordingly the analytical results for that case will be presented first, followed by the case of two incident longitudinal (P) waves, and finally incident P and SV waves. The angle of incidence will be denoted by q when measured from the x-axis, and by j when measured from the y-axis, as indicated in Fig. 2 (a). The angle between the directions of propagation of the incident waves is referred to as the interaction angle f,

  A B f ¼ jA þ jB ¼ p  q þ q :

(10)

For incident SV waves, the direction of polarisation is indicated by an arrow, as shown in Fig. 2 (a). It can be seen from the direction of the arrows in Fig. 2 (a) that there is a phase difference of p between the two incident waves. This is designed to ensure a symmetrical configuration with respect to the y-axis when both waves have the same amplitude, frequency and

Fig. 2. Configuration for non-collinear shear-wave mixing at a plane contact interface, defining the incident wave propagation directions and the interaction angle (a), with the pattern of reflected, transmitted and mode-converted waves for wave A and B incident separately shown in (b) and (c) respectively. Here, k refers to wave vectors and q to angles of incidence relative to the x-axis, with subscript L and T referring to longitudinal and shear waves respectively, and superscript A and B referring to waves associated with the incident waves A and B respectively.

6

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

angle of incidence. This phase difference must be kept in mind below when evaluating the horizontal components of the wave vectors. Each incident SV wave in isolation will, in general, generate reflected, refracted and transmitted waves as indicated in Fig. 2 (b,c). These scattered waves include harmonics of various orders that propagate at the same speed as the fundamental, and therefore at the same angle relative to the interface. The problem considered below is to characterise the additional waves due to interface wave mixing that arise when two incident waves interact at the interface. 3.1. Scattered field representation The fundamental concept in Ref. [29] is that the scattered field generated at an interface can be expressed in terms of the relative displacements Dux and Duy across the interface. Furthermore, a general scattered field can be decomposed into two contributions exhibiting different symmetries with respect to the interface ðy ¼ 0Þ that result in particular combinations of displacement and stress boundary conditions at the interface. First, for a Mode I field, ux is an even function of y, whereas uy is an odd function of y, i.e.

ux ðx; y; tÞ ¼ ux ðx; y; tÞ; uy ðx; y; tÞ ¼ uy ðx; y; tÞ:

(11a)

Consequently, sxy is an odd function ofy, which, combined with the requirement that stresses are continuous across the interface, automatically satisfies the interface condition

sxy ðx; y ¼ 0; tÞ ¼ 0:

(11b)

Furthermore, ux is continuous across y ¼ 0; whereas uy is discontinuous. Thus,

Dux ðx; y ¼ 0; tÞ ¼ 0;

(11c)

and the scattered field can be completely determined by specifying the Mode I relative displacement Duy across the interface. For a time-harmonic field, this relative displacement can be expressed as follows,

Duy ðx; y ¼ 0; tÞ ¼ VeiðxxutÞ ;

(11d)

where V denotes the complex amplitude, which needs to be determined for particular incident wave combinations, u the angular frequency, and x the horizontal wave number. By employing integral transform techniques [33], the displacement field corresponding to the symmetry and boundary conditions in Eq. (11 a-d) can be expressed as follows, for y > 0;

ux ðx; y; tÞ ¼ 

  x  2 2 iay iby 2 e VeiðxxutÞ ; x  k þ 2 ab e T 2ak2T

  1  2 2 uy ðx; y; tÞ ¼  2 2x  k2T eiay  2x eiby VeiðxxutÞ ; 2kT qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ k2L  x2 ; b ¼ k2T  x2 ; kL ¼ u cL ; kT ¼ u cT ;

(12)

where cL ; cT denote the longitudinal and transverse wavespeed respectively. The displacement field for y < 0 follows from Eq. (11a), and the corresponding stresses are readily derived from these displacements [29]. Secondly, for a Mode II field, ux is an odd function of y, whereas uy is an even function, i.e.

ux ðx; y; tÞ ¼ ux ðx; y; tÞ; uy ðx; y; tÞ ¼ uy ðx; y; tÞ;

(13a)

so that the associated normal stress at the interface vanishes,

syy ðx; y ¼ 0; tÞ ¼ 0: Furthermore, uy is continuous across y ¼ 0; whereas ux is discontinuous. Thus,

(13b)

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Duy ðx; y ¼ 0; tÞ ¼ 0;

7

(13c)

and the scattered field can be completely determined by specifying the Mode II relative displacement Dux across the interface. For a time-harmonic field, this relative displacement can be expressed as follows,

Dux ¼ UeiðxxutÞ ;

(13d)

where again the complex amplitude Uhas to be determined for particular incident wave combinations. The associated displacement field for y > 0 is given by

ux ðx; y; tÞ ¼

  1  2 iay  2 2x e  2x  k2T eiby UeiðxxutÞ ; 2 2kT

   x  2 2 iay iby 2 e UeiðxxutÞ : ab e þ 2 x  k uy ðx; y; tÞ ¼ T 2bk2T

(14)

The displacements for y < 0 follow from Eq. (13a), and the corresponding stresses are readily derived [29]. It is emphasized that the representations in Eqs. (12) and (14) apply for all constituents of the scattered field, including in particular the mixed-frequency waves (if present), provided that the appropriate values are used for the frequency u and horizontal wavenumber x: A fundamental principle for single-wave scattering at a plane interface is that, for all individual scattered waves, the horizontal wavenumber (i.e. the component of the wave vector parallel to the interface) must be equal to the horizontal wavenumber of the incident wave. This principle can be recognised as Snell's law [33,34] for the case of a welded interface; it also applies for a contact interface governed by a general traction law [29]. In particular, for incident or scattered P waves, x ¼ ±kL cosqL ; a ¼ kL sinqL ; whereas for SV waves x ¼ ±kT cosqT ; b ¼ kT sinqT ; where q denotes the angle between the wave vector and the x-axis, as indicated in Fig. 2. When employing the representations in Eqs. (12) and (14) for the individual waves constituting the scattered field due to a contact interface, one must determine three quantities, viz. (i) the angular frequency u, which includes the mixed frequencies for the case of two or more incident waves, as well as harmonics of the input waves; (ii) the horizontal wavenumber x; which determines the angle of propagation according to the wave mode (P or SV), as noted above; and (iii) the complex amplitudes U, V, of the Mode I and Mode II contributions. The final results for various combinations of input waves are summarised next, prior to discussing their practical implications; the details of the derivations are relegated to an Appendix for clarity of presentation. 4. Analytical results for wave mixing at a nonlinear contact interface Consider the interaction of two incident plane waves at a contact interface, as indicated in Fig. 2 (a) for the case of SV waves. The total acoustic field u can be written as the sum of an incident field uInc , which is prescribed, and a scattered field uS , containing both Mode I and Mode II contributions, which is to be determined,

u ¼ uInc þ uS :

(15a)

Once the amplitudes U, V, of the Mode I and II contributions have been determined, the corresponding amplitudes ASL and AST

of the longitudinal (P) and transverse (SV) scattered waves can be obtained from Eqs. (12) and (14) as follows, for y > 0;




 kL  2 kL x


2 2x  kT V þ 2 U
; ¼

2ak2 kT
T  


x
2x2  k2T

S AT ¼
V þ U
:
kT
2bkT ASL

(15b)

For an interface characterised by a general traction law as in Eq. (4a), it can be assumed that the relative displacements across the interface can also be expressed as polynomials in ε;

Dux ¼ Duy ¼

∞ X n¼0 ∞ X n¼0

ðnÞ

εn Dux ; (16a) ðnÞ

εn Duy :

8

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

A perturbation analysis then consists in a successive approximation scheme obtained by grouping terms in the traction law according to successive powers of ε [35], as detailed in the Appendix. The first approximation (or zeroth-order solution), obtained by retaining only the terms involving ε0 , is given by









Inc sSxy Duxð0Þ  K T0 Duð0Þ x ¼ sxy ;

(16b)

sSyy Duyð0Þ  K N0 Duyð0Þ ¼ sInc yy :

This corresponds to the case of linear springs, for which linear superposition applies, i.e. the relative displacements are simply the sum of the relative displacements produced by each incident wave separately. The second approximation (or firstorder correction) is governed by the following equations,









sSxy Duxð1Þ  K T0 Duð1Þ x ¼ 0;



2

sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 Duyð0Þ :

(16c)

Before proceeding to state the results for specific combinations of incident waves, it is pertinent to note two important features of these equations. First, the stress boundary conditions in Eqs. (11b) and (13b) and for Mode I and II fields imply that the two equations in Eq. (16b) are uncoupled, as are the two equations in Eq. (16c). It is this feature that enables explicit analytical solutions to be derived. Secondly, it can be seen from Eq. (16c) that the zeroth-order solution provides the source terms for the first-order correction. As noted above, the zeroth-order solution corresponds to linear springs, and therefore involves no wave mixing or higher harmonics, both of which arise from the nonlinear terms in the traction law. However, the amplitude and other characteristics of the nonlinear response are very much dependent on the linear-springs response, because of its presence in the source terms of Eq. (16c). In particular, the linear springs introduce a characteristic frequency whose implications will be further discussed below. 4.1. Two incident shear waves Consider first the case of two incident SV waves, identified as A and B, as indicated in Fig. 2(a). The incident field for wave A is given by

uAx ¼ sinqA AA eiðxA xbA yuA tÞ ; uAy ¼ cosqA AA eiðxA xbA yuA tÞ ;

(17a)

xA ¼ kA cosqA ; bA ¼ kA sinqA ; kA ¼ uA =cT : For this incident wave in isolation, the zeroth-order solution for the complex amplitudes of the relative displacements across the interface is given by Ref. [29], ð0Þ UA

ð0Þ

  2 2kA bA 2xA  k2A . AA ; ¼ FðxA Þ þ 2ik2A bA K T0 m

VA ¼

4kA xA aA bA FðxA Þ þ 2ik2A aA K N 0

. AA ;

(17b)

m

2  2 2 FðxA Þ ¼ 4xA aA bA þ 2xA  k2T ; where FðxÞ denotes the Rayleigh function (i.e. FðxR Þ ¼ 0 for xR equal to the wavenumber for a Rayleigh wave of frequency u), and m denotes the shear modulus. Some care is required when evaluating aA , and hence also FðxA Þ, for angles of incidence in the range 0 < qA < qc ; where qc ¼ cos1 ðcT =cL Þ denotes the critical angle of incidence. In that range, the scattered longitudinal waves (both reflected and transmitted) propagate as inhomogeneous plane waves, with the direction of propagation aligned with the interface, and with the amplitude decaying exponentially with distance normal to the interface. Consequently,

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 k2L  xA ; qc < qA < p 2; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ i x2A  k2L ; 0 < qA < qc ; kL ¼ uA cL :

9

aA ¼

(17c)

The incident field for wave B is given as follows, keeping in mind the direction of polarisation as shown in Fig. 2 (a),

uBx ¼ sinqB AB eiðxB xbB yuB tÞ ; uBy ¼ cosqB AB eiðxB xbB yuB tÞ ;

(18a)

xB ¼ kB cosqB ; bB ¼ kB sinqB ; kB ¼ uB =cT ; and the zeroth-order complex amplitudes for this incident wave in isolation are now given by ð0Þ UB

ð0Þ

  2 2kB bB 2xB  k2B . AB ; ¼ FðxB Þ þ 2ik2B bB K T0 m

VB ¼ 

4kB xB aB bB . AB ; FðxB Þ þ 2ik2B aB K N 0 m

(18b)

2  FðxB Þ ¼ 4x2B aB bB þ 2x2B  k2T ; ð0Þ

where aB is evaluated using an obvious modification of Eq. (17c). It can be seen that the tangential displacement U B

is

ð0Þ UA :

In particular, for the symmetrical configuration where both incident waves have the same opposite to frequency, amplitude and angle of incidence, the combined tangential displacement is zero, i.e. the scattered field is purely Mode I. The source terms for the first-order correction in Eq. (16c) can now be obtained from Eqs. (17,18), keeping in mind that Eq. (16c) involves physical quantities that are necessarily real-valued,

h

i

h

i

Duxð0Þ ¼ < U Að0Þ eiðxA xuA tÞ þ U Bð0Þ eiðxB xuB tÞ ;

(19a)

Duyð0Þ ¼ < V Að0Þ eiðxA xuA tÞ þ V Bð0Þ eiðxxB uB tÞ ; where < denotes the real part. It will prove convenient to introduce the notation ð0Þ

A

ð0Þ

A

ð0Þ

B

ð0Þ

B

U A ¼ U A0 eic0 ; V A ¼ V A0 eiz0 ;

(19b)

U B ¼ U B0 eic0 ; V B ¼ V B0 eiz0 :

Equation (19a) can then be written in terms of real-valued quantities as follows,









Duxð0Þ ¼ U A0 cos xA x  uA t þ cA0 þ U B0 cos xB x  uB t þ cB0 ; Duyð0Þ ¼ V

A 0

  A cos xA x  uA t þ z0 þ V

B 0

  B cos xB x  uB t þ z0 :

Thus, Eq. (16c) for the first-order correction takes the form,

(19c)

10

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078









sSxy Duxð1Þ  K T0 Duxð1Þ ¼ 0;

(20a) h

sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 V  2 V A0 2  2 V B0

A 0

2 .  2 . 2 þ V B0 2þ

  A cos 2 xA x  uA t þ z0 þ (20b)

  cos 2 xB x  uB t þ zB0 þ

V A0 V

B 0

2 o n cos ðxA þ xB Þx  ðuA þ uB Þt þ zA0 þ zB0 þ

V A0 V

B 0

n oi A B cos ðxA  xB Þx  ðuA  uB Þt þ z0  z0 :

It is clear from Eq. (20a) that there is no Mode II contribution, which is a consequence of K T1 ¼ 0 (as discussed in Section 2), whereas the source term for the Mode I contribution in Eq. (20b) consists of the sum of six terms: (i) two zero-frequency terms corresponding to static displacements; (ii) 2 s-harmonic terms at twice the incident frequency of wave A and wave B respectively; and (iii) two terms at the sum and difference frequencies. It is these last two terms that are characteristic of wave mixing due to interface nonlinearity. Accordingly, we shall focus here on characterising the scattered field due to these mixed-frequency source terms, but it is emphasized that the other terms will also contribute to the observed physical response. 4.1.1. Mixed wave at sum frequency for SV/SV incidence By retaining only the sum-frequency source term in Eq. (20b), the governing equation for the Mode II contribution can be rewritten as follows, exploiting the fact that the equation is linear,





h

i

sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 V A0 V B0 < eiððxA þxB ÞxðuA þuB Þtþz0 þz0 Þ : A

B

(21)

This equation has the same form as for the zeroth-order problem (cf. Appendix), which leads to the following solution for ð1Þ

the complex amplitude V ABþ of the Mode I contribution,

 þ 2 2iaABþ kAB KN m T 1

ð1Þ

V ABþ ¼

FðxABþ Þ þ 2iaABþ



 þ 2 kAB KN T 0

V A0 V B0 ;

m

(22)

.

þ

xABþ ¼ xA þ xB ; kAB ¼ ðuA þ uB Þ cT ; T where aABþ has to be evaluated in accordance with Eq. (17c) depending on the direction of the mixed wave. It is recalled that a Mode I relative displacement generates both P and SV waves, in accordance with Eq. (12). The directions of propagation of the ð1Þ

mixed P and SV waves arising from V ABþ are given by þ



.

þ



þ

.

þ



.

þ



þ

.

; kAB qAB ¼ cos1 xABþ kAB ¼ ðuA þ uB Þ cL ; L L L

(23)

; kAB qAB ¼ cos1 xABþ kAB ¼ ðuA þ uB Þ cT ; T T T ð1Þ

and their amplitudes are obtained by substituting V ABþ from Eq. (22) into Eq. (15b). These amplitudes and propagation angles complete the characterisation of interface wave mixing for the configuration in Fig. 2 (a). gimes can be identified, as follows: (i) The wave pattern for the sum-frequency mixed waves is shown in Fig. 3 (a). Three re þ

for 0 < jxABþ j < kAB L ; both P and SV waves are generated at the sum frequency uA þ uB , with propagation angles given by Eq. þ

þ

þ

þ

(23) in the ranges 0 < qAB < p=2; 0 < qAB < qc ; (ii) for kAB < jxABþ j < kAB L T L T ; a propagating SV wave is generated, with propagation þ

angle in the range qc < qAB < p=2; as well as an inhomogeneous longitudinal plane wave propagating along the interface, and T (iii) for

þ jxABþ j > kAB T ;

no propagating waves are generated. This classification is the analogue for interface wave mixing of the

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

11

Fig. 3. Wave pattern for mixed waves from SV/SV incidence. Direction of propagation of the mixed longitudinal and shear waves at (a) the sum of the incident frequencies, and (b) the difference of the incident frequencies. Here, k refers to wave vectors and q to angles of incidence relative to the x-axis, with subscript L and T referring to longitudinal and shear waves respectively, and superscript AB± referring to mixed waves at the sum or difference frequency respectively. The sum and the difference of the x-component of the incident shear wave vectors xA and xB defines the x-component of the wave number for the mixed waves xA ± xB , and hence their direction.

Jones and Kobett [8] formula for the interaction angle for bulk wave mixing. Before discussing the differences between these two cases in detail (see Section 5), it will be convenient to record the results for various other cases of wave mixing. 4.1.2. Mixed wave at difference frequency for SV/SV incidence Consider next the source term involving the difference of the incident frequencies, where it is assumed, without loss of generality, that uA > uB . Following the same approach as outlined above, we can derive the following solution for complex amplitude of the Mode I relative displacement at the difference frequency

ð1Þ

V AB ¼

  2 2iaAB kAB KN m T 1

A B   2 V 0 V 0 ; AB N FðxAB Þ þ 2iaAB kT K0 m

(24)

.



xAB ¼ xA  xB ; kAB ¼ ðuA  uB Þ cT : T The amplitudes of the resulting P and SV waves can again be obtained by substituting Eq. (24) into Eq. (15b), and the directions of propagation are given by 



.







.





.







.

; kAB qAB ¼ cos1 xAB kAB ¼ ð uA  uB Þ c L L L L

(25)

; kAB qAB ¼ cos1 xAB kAB ¼ ð uA  uB Þ c T : T T T

gimes can again be identified, as discussed above for the sum frequency. The wave pattern is shown in Fig. 3 (b). Three re 4.2. Mixed wave at sum and difference frequencies for P/P incidence Consider next the case of two incident longitudinal waves A and B as shown in Fig. 4 and described by Eq. (26 a,b).

uAx ¼ cosqA AA eiðxA xaA yuA tÞ ; uAy ¼ sinqA AA eiðxA xaA yuA tÞ ;

xA ¼ kA cosqA ; bA ¼ kA sinqA ; kA ¼ uA =cL ;

(26a)

12

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Fig. 4. Wave pattern for mixed waves from P/P incidence. Direction of propagation of the mixed longitudinal and shear waves at (a) the sum of the incident frequencies, and (b) the difference of the incident frequencies. Here, k refers to wave vectors and q to angles of incidence relative to the x-axis, with subscript L and T referring to longitudinal and shear waves respectively, and superscript AB± referring to mixed waves at the sum or difference frequency respectively. The sum and the difference of the x-component of the incident shear wave vectors xA and xB defines the x-component of the wave number for the mixed waves xA ± xB , and hence their direction.

uBx ¼ cosqB AB eiðxB xaB yuB tÞ ; uBy ¼ sinqB AB eiðxB xaB yuB tÞ ;

(26b)

xB ¼ kB cosqB ; bB ¼ kB sinqB ; kB ¼ uB =cL : The solution of the zeroth-order problem for the relative displacements is now found to be









Duxð0Þ ¼ U A0 cos xA x  uA t þ cA0 þ U B0 cos xB x  uB t þ cB0 ; Duyð0Þ ¼ V

A 0

  A cos xA x  uA t þ z0 þ V ð0Þ

A

B 0

  B cos xB x  uB t þ z0 : ð0Þ

(27)

A

where the relevant amplitudes U A ¼ U A0 eic0 , V A ¼ V A0 eiz0 etc. for each of the incident waves in isolation are given in the Appendix, Eqs. (A.5,6). Following the same approach as for SV/SV incidence (see Appendix for details), the complex amplitude solutions for the Mode I opening is found to be

ð1Þ

V AB± ¼

 ± 2 2iaAB± kAB KN m T 1 FðxAB± Þ þ 2iaAB±



 ± 2 kAB KN 0 T

V A0 V B0 ;

m

(28)

.

±

xAB± ¼ xA ±xB ; kAB ¼ ðuA ±uB Þ cT : T The amplitudes of the mixed-frequency P and SV waves are obtained by substituting Eq. (28) into Eq. (15b), and the directions of propagation are given by ±



.

±



±

.

±



.

±



±

.

; kAB qAB ¼ cos1 xAB± kAB ¼ ðuA ±uB Þ cL ; L L L

(29)

; kAB qAB ¼ cos1 xAB± kAB ¼ ðuA ±uB Þ cT : T T T

Again, three regimes can be identified for the scattering of the mixed wave at the sum or difference frequency, as shown in Fig. 4 (a,b). It is important to note that although the expressions above have the same form as the ones obtained in Section 4.1 for the incident shear waves, the values of xA and xB are not the same in the case of two incident P waves, and therefore the directions of propagation are not the same. Similarly, the 0th order normal and tangential gap opening amplitudes V

A 0

and

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

13

V B0 are different, leading to different scattered wave amplitudes although the 1st order solution is expressed in the same form. 4.3. Mixed wave at sum and difference frequencies for P/SV incidence Finally, consider the case where wave A is a shear wave given by Eq. (17a) and wave B is a P wave given by Eq. (26b). Proceeding as above, the amplitude of the Mode I relative displacement characterising the mixed waves is found to be

ð1Þ

V AB± ¼

 ± 2 2iaAB± kAB KN m T 1

A B  ± 2 V 0 V 0 ; N FðxAB± Þ þ 2iaAB± kAB K m T 0

(30)

.

±

xAB± ¼ xA ±xB ; kAB ¼ ðuA ±uB Þ cT : T ð0Þ

A

where the relevant amplitudes U A ¼ U A0 eic0 , etc. for the incident SV and P wave in isolation are given in the Appendix. The amplitudes of the mixed P and SV waves are obtained by substituting Eq. (30) into Eq. (15b), and the directions of propagation are again given by Eq. (29), but with xAB± as given in Eq. (30). 5. Characteristics of interface wave mixing compared with bulk mixing The preceding analysis indicates a fundamentally different mechanism for wave mixing due to interface nonlinearity relative to bulk nonlinearity. To discuss the difference in more detail, it is first recalled that bulk mixing relies on an internal resonance that requires a spatial matching of the spacing between the intersecting wave crests of the incident waves and the wavelength of a mixed wave at the sum frequency uA þ uB. The direction of this mixed wave is thus given by the wave vector [8]. þ

kAB ¼ kA þ kB :

(31a)

This resonance condition leads to the following formula for the interaction angle between the incident waves [8],

cos f ¼ k2 

o  n 1  k2 1 þ ðuB =uA Þ2 2uB =uA

, ; k ¼ cT

cL :

(31b)

This formula only involves the frequency ratio uB =uA , and not the individual frequencies, because the bulk linear response does not entail a characteristic length or characteristic time, and hence does not entail a characteristic frequency. Finally, there is a cumulative build-up of amplitude within the interaction domain, which leads to the observed amplitude of the mixed wave being proportional to the size of that interaction domain. None of these characteristics apply to interface mixing. First, the mixing does not involve an internal resonance. Secondly, the mixing does not give rise to a single wave mode at a single frequency and propagating in a single direction. Instead, both P and SV mixed waves are generated, with the waves propagating in both the forward direction (transmitted waves) and the backward direction (reflected waves), and at both the sum and the difference frequencies, as documented in Section 4. Thirdly, the direction of propagation of these mixed waves is now determined by the sum or difference of the x-components of the incident wave vectors (i.e. the components parallel to the interface), instead of the sum of the wave vectors as in Eq. (31a). However, the difference between these two directions may not be large in practice, as will be seen in Section 6. Finally, there is no longer a simple formula for the interaction angle giving rise to the maximum amplitude for any particular mixed wave. Before considering in more detail the dependence of amplitude on interaction angle, it is worth noting that for both bulk and interface mixing, the analytical approach employs a perturbation analysis, i.e. the nonlinear response is assumed to be a small correction whose amplitude depends on a source term generated by the linear response. However, the response of an interface bridged by linear springs involves one or more characteristic frequencies (unlike the bulk linear response, as noted above). For the case of interest in the present work, where the nonlinearity is restricted to the tension/compression response at the interface (Section 2), the relevant characteristic frequency is most readily identified by considering a normally incident P wave, for which the reflection and transmission coefficients (R and T) derived in Refs. [19,21,25] can be re-written as follows,

1 1 R ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; T ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 2 1 þ ðu0 =uÞ 1 þ ðu=u0 Þ2 2K N u0 ¼ 0 ; rcL

(32)

14

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Fig. 5. (a) Relations between the linear stiffness coefficient K N 0 , the characteristic frequency f0 and the contact stress p0 for rough-surface contact. (b) Reflection and transmission coefficients (R and T) at an interface bridged by linear springs for a normal incident P wave, as a function of normalised frequency f = f0 .

where r denotes the density of the bulk material. The characteristic frequency f0 ð ¼ u0 =2pÞ depends on the stiffness coefficient K N 0 ; which in turn depends on the contact stress p0 for the case of rough-surface contact (Section 2), as shown in Fig. 5 (a). The linear stiffness coefficient K N 0 is often determined experimentally [23,25] by matching the measured reflection and/or transmission coefficients to the relations in Eq. (32), which are shown in Fig. 5 (b). Thus a systematic investigation of interface mixing needs to document the relevant characteristic frequency f0 (or equivalently the stiffness coefficient K N 0 , or the contact stress p0 ), as well as the individual input frequencies. To illustrate how the present analysis can guide the design and interpretation of experiments, consider SV/SV incident waves, as shown in Fig. 2 (a), for the special case where the two waves have the same amplitude, frequency and angle of incidence qA ¼ qB ¼ qinc : This symmetrical configuration has been investigated experimentally for both. Bulk and interface mixing [1,17]. In both cases, the only propagating mixed wave is a P wave at double frequency (i.e. sum of incident frequencies), but interface mixing generates both forward and backward propagating mixed waves, of equal amplitude, propagating in the direction of the negative and positive y-axis, as noted in Section 4.1.1, whereas bulk mixing only generates a forward propagating mixed wave. Furthermore, for interface mixing, the amplitude of this mixed wave depends on both the interaction angle fð ¼ p 2qinc Þ and the normalised frequency f =f0 ; whereas for bulk mixing there is no characteristic frequency, and consequently no dependence on frequency if the incident stress amplitude is used as the basis for comparison (noting that the stress amplitude is proportional to the frequency for a fixed displacement amplitude). The variation of the mixed P wave amplitude with interaction angle for interface mixing is shown in Fig. 6 for the particular choice of incident amplitude AA ¼ AB ¼ 10 nm, frequency fA ¼ fB ¼ 5 MHz; contact stress p0 ¼ 0:2 MPa; and bulk properties representative of aluminium alloy: Young's modulus E ¼ 70 GPa; Poisson ratio n ¼ 0:33; density r ¼ 2700 kg=m3 : For this load and material, the characteristic frequency is f0 ¼ 3:2 MHz. The scattered wave attains its maximum amplitude at an interaction angle of 70
, as shown in Fig. 6. This variation differs radically from that expected for bulk mixing, which would show a resonance peak for f ¼ 120
according to Eq. (31b), (cf. Fig. 4 in Ref. [17]), but the overall shape agrees with recent experimental observations for interface mixing (cf. Fig. 9 in Ref. [17]).

Fig. 6. Mixed P wave amplitude characteristics for symmetrical SV/SV incidence as a function of the interaction angle, with zero amplitude at the critical angle fc ¼ p  2qc ¼ 60
, followed by a maximum shortly after obtained for fopt ¼ 70
.

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

15

The distinctive features in Fig. 6 are (i) the amplitude falls to zero at fc ¼ 60+ ; corresponding to the critical incident angle

qinc c

¼ 60+ for the chosen value of Poisson ratio, and (ii) the amplitude reaches a maximum Amax for an optimal interaction angle fopt ¼ 70
just beyond the critical angle. An important feature revealed by the results presented in Section 4.1.1 is that both Amax and fopt should depend on the normalised frequency. This dependence is shown in Fig. 7, where the optimal interaction angle and the corresponding maximal mixed P wave amplitude are evaluated for the symmetric configuration considered above. In experimental work involving rough-surface contact, the parameter that is most easily controlled is the contact stress p0 ; and it was reported in Ref. [17] that varying p0 changes both Amax and fopt . However, a change in p0 also changes the linear stiffness parameter K N 0 and the characteristic frequency f0 ; and hence the normalised frequency f = f0 : Thus, the present analysis suggests that future experimental work should keep track of the change in f0 and present the results in the form shown in Fig. 7. It is worth noting that the perturbation analysis is restricted to weakly nonlinear traction law and it is thus expected that the amplitude of the mixed waves is accurately predicted for relatively low level of nonlinearity. However, the other characteristics of the interface wave mixing presented in the previous sections also apply for larger nonlinearity, including the case of clapping. In particular, it remains true that (i) the scattered field can be expressed in terms of the relative displacements across the interface, which can be decomposed into Mode I and Mode II contributions; (ii) the wave mixing gives rise to source terms involving both sum and difference frequencies; and (iii) the mixed frequency waves propagate both forwards and backwards (corresponding to transmitted and reflected waves), with the angles of propagation being determined from the sum or difference of the components of the incident wave vectors parallel to the interface, which is a generalisation of Snell's law to interface wave mixing. The computational results presented in the next section will be seen to confirm the predictions of the analytical model presented above.

6. Finite element modelling The objective of this section is to present finite element (FE) simulations that provide an independent validation of the interface wave-mixing characteristics discussed in Section 5 based on the analytical model. Unlike the analytical model, which relies on a perturbation analysis, these FE simulations are not restricted to weak nonlinearity. The difference between the two will therefore be seen to become more apparent with increasing input amplitude. A representative FE model (2D plane strain) is shown in Fig. 8. It consists of two blocks that are in contact across a plane interface located at y ¼ 0: The interface traction law in Eqs. (4b) and (7) is implemented by using the subroutine VUINTER in the commercial FE package ABAQUS [36], with the values of the stiffness coefficients being determined prior to each simulation using the relations given in Section 2. The penalty contact algorithm is used for the solution with default parameters. The incident waves A and B are generated by imposing normal (P wave) or tangential (SV wave) displacements on a 120 mm face located 100 mm away from the origin and oriented to emit the waves with the desired angle of incidence. The dimension of the input face and the use of spatial apodisation ensure near plane wave conditions, as shown in Fig. 8(a). The geometry of the model is modified for each considered case, but the above geometrical characteristics are kept unchanged. The incident waves are 10 nm amplitude N-cycles tone burst, where the number of cycles N and the centre frequency f (chosen in the range of [0.4e2] MHz) will be specified explicitly later on for each considered case. The rising and fading of the tone burst are defined by a half Hann window covering ¼ of the pulse duration as shown in Fig. 9(a) for the particular case of a 12-cycles 1 MHz excitation. This tone burst excitation ensures a constant amplitude over several cycles, hence providing a

Fig. 7. Wave mixing characteristics for symmetrical SV/SV incidence. Variation of optimum interaction angle and corresponding maximum amplitude of the mixed P wave at the sum frequency versus normalised frequency f =f0 varied through changes of the load p0 .

16

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Fig. 8. Example of FE model configuration for the case of two incident SV waves A and B, showing the position of the interface, dimensions, output segments and boundary conditions. (a) and (b) show respectively the displacement field before and after scattering at the interface.

Fig. 9. (a) Example of a 12-cycles 1 MHz tone burst with time window, (b) corresponding spectrum.

better approximation of the constant amplitude harmonic excitation assumed in the analytical model than a simple Hann windowed tone burst. The output nodal displacements are recorded over 40 mm long segments located 40 mm and 60 mm away from the origin to capture respectively the scattered mixed SV and P waves, as shown in Fig. . These output segments are oriented perpendicularly to the expected directions of propagation from Section 4. To establish that the scattered mixed waves are propagating along the expected directions of propagation, a synthetic beam steering is performed at reception using successive groups of output nodes, with 25 mm aperture and 1 mm step. For each node within the group, the tangential or radial displacement is computed from the nodal displacements depending on whether a mixed SV or P wave is considered. A Hann window is used to select the time-domain data, based upon the expected time of flight for the considered mixed wave. This signal is then filtered to keep only the mixed wave frequency, using a 100 kHz bandwidth. Then, a time delay is applied to the signal (performed in the phase domain) depending on the position of the node relative to the centre of group and its location with respect to the origin. The delayed signals are then summed node-wise to complete the beam steering operation. Finally, the envelope of the remaining signal is obtained from a Hilbert transform, and the amplitude of the scattered mixed wave is evaluated from the peak value of the envelope. These processing steps and the beam steering principle are summarised in Fig. 10. As a consequence, the mixed wave amplitude is computed for multiple directions defined between the origin and the centre of each group of output nodes, which provides the scattering pattern for the mixed wave. The blocks are assumed to be isotropic, with material properties for aluminium alloy, as given above in Section 5. The linear and quadratic bulk material viscosity were set to zero in order to remove unwanted damping in view of the comparison with the analytical solutions, and to ensure energy conservation. Adequate spatial discretization is essential in the FE method. In order to have an accurate solution for the mixed P and SV waves at the sum or difference frequency, the corresponding wavelengths were discretised into at least 10 elements [37]. The maximal element dimension was set at 0.08 mm, which ensures at least 15 elements for the shortest considered wavelength (mixed SV wave at the sum frequency of 2.6 MHz), and therefore provides a sufficient discretization for all other considered frequencies. The FE mesh consists of 4-node bilinear plane strain quadrilateral elements CPE4R with reduced integration and hourglass control, and 3-node linear plane strain triangle elements CPE3 where required. The time explicit resolution involves a Courant-Friedrichs-Lewy stability condition regarding the definition of the time step dt  a=cL , where a corresponds to the smallest element dimension and cL to the

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

17

Fig. 10. (a) Post-processing steps to obtain the scattered wave amplitude for each group of nodes used for synthetic beam steering at reception, (b) beam steering principle.

longitudinal wave velocity in the medium. However, the presence of the large interface stiffness imposes a more stringent stability condition for the time step. A convergence study in time has been carried out and the actual time step was taken as dt ¼ 2:5 ns to ensure convergence for all the considered cases. The simulations are run over 70 ms whenever one of the incident waves is a shear wave, and 40 ms otherwise. 6.1. Existence of mixed waves and directions of propagation Six different configurations including the three types of incident mode pairs SV/SV, P/P and SV/P are considered here to demonstrate the existence of the mixed P and SV waves at the sum or difference frequency, and to verify that the analytical solution provides the correct direction of propagation and amplitude. An FE model was constructed for each of these 6 cases using the parameters indicated in Table 1. The choice of the incident angles and frequencies were selected to obtain different scattering regimes based on the value of xAB± as discussed in Section 4.1.1. The compression load at the interface is set at p0 ¼ 0:5 MPa, and both incident waves have a 10 nm amplitude. In the cases where two different incident wave modes are considered (i.e. SV/P or P/SV), the generation of the P wave is delayed in order to ensure a concurrent arrival at the interface with the SV wave. The number of cycles in the excitations is selected in accordance with the frequency and angle of incidence in order to ensure a sufficient overlapping of the incident waves at the interface. The nodal displacements acquired over the output lines are used to perform the beam steering at reception and plot the scattering patterns for the mixed waves, which are normalised by the incident wave amplitude and shown in Fig. 11(aef) respectively for each case of Table 1. The corresponding analytical solutions (direction and amplitude of the mixed waves) are reported for each scattering pattern in Fig. 11. These results demonstrate that for all considered cases, the expected scattered mixed waves are generated, and the direction of propagation and the amplitude of the mixed waves are correctly predicted by the analytical solutions in Section 4. In particular, for SV/SV-1, where two incident shear waves interact to generate a mixed P wave at the sum frequency, the numerical results confirm that the directions of propagation are defined solely by the x-component of the wave vector along the interface, as noted in Section 4, and are not defined by the sum of the incident wave-vectors as would be the case for bulk nonlinearity [8], and as assumed in previous work on interface mixing [17,18]. However, these directions differ by only a few degrees, and the difference decreases when the configuration tends to be symmetric (same angle of incidence and same

Table 1 Definition of the incident waves A and B, including the mode pair, angles of incidence, frequencies and numbers of cycles. For each case, the directions of propagation of the mixed waves at the sum or difference frequencies are given when existing, with X indicating that the corresponding mixed wave is not allowed.

SV/SV-1 Fig. 11 (a) SV/SV-2 Fig. 11 (b) P/P-1 Fig. 11 (c) P/P-2 Fig. 11 (d) SV/P-1 Fig. 11 (e) P/SV-2 Fig. 11 (f)

±

±

qIA (
)

qIB (
)

fA (MHz)

fB (MHz)

NA

NB

fAB± (MHz)


qAB L ( )


qAB T ( )

50

40

1

0.6

12

7

80

35

2

0.4

30

9

50

60

1

0.6

10

6

70

40

2

0.6

20

6

50

40

1

0.6

20

6 (delay)

70

40

1.5

0.6

20 (delay)

12

0.4 1.6 1.6 2.4 0.4 1.6 1.4 2.6 0.4 1.6 1.1 2.1

X 76.9 33.1 89.1 X 77.6 35.2 85.1 X 59.3 X 101

X 83.4 65.1 89.5 X 83.8 65.7 87.5 X 75.1 37.1 95.5

18

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

Fig. 11. Scattering patterns of the mixed P and SV waves at the sum and/or difference frequency for each case considered in Table 1 obtained from FE simulations. The solid green disks correspond to the analytical results. The dashed line in (a) represents the direction of propagation of the mixed P wave for bulk material nonlinearity. (Color online). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

frequencies). The second important observation is that a mixed SV wave is also generated at the sum frequency in the case of a nonlinear contact interface, as opposed to the case of bulk material nonlinearity. However, the amplitude of this mixed SV wave is much lower than that of the mixed P wave, which remains the most attractive option for practical applications. For the second case SV/SV-2 shown in Fig. 11 (b), the configuration of the incident waves is such that both SV and P mixed waves at the difference frequency are also generated in addition to the scattered wave at the sum frequency. The scattering of the mixed waves at the difference frequency for two incident SV waves is not predicted for classical material nonlinearity [8]. The amplitude of the scattered SV and P wave at the difference frequency are slightly larger than the mixed P wave amplitude

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

19

at the sum frequency for this particular case, but much lower than the mixed P wave amplitude obtained previously in Fig. 11 (a). The amplitude of the mixed SV wave at the sum frequency is negligible here. The third case P/P-1 shown in Fig. 11 (c) corresponds to two incident P waves that are defined for the generation of the scattered waves at the sum frequency only, as indicated in Table 1. Again, both mixed SV and P waves are generated, which does not occur for material nonlinearity where only a mixed SV a wave at the difference frequency is predicted [8]. This interaction results in scattering patterns very similar to the ones obtained previously in Fig. 11 (a). This suggests that single face access wave mixing can be realised with two P waves in the case of a contact interface, which can be very interesting in practice as larger mixed wave amplitude are obtained in Fig. 11 (a) and (c). However, because the mixed wave mode is the same as the incident waves, the mode selectivity feature (i.e. the generation of a different mode) is lost. Therefore, if one wants to exploit the mixed P wave at the sum frequency, it is advisable to use different angles of incidence and different frequencies to retain the benefit of the direction and frequency selectivity. Next, the two incident P waves are chosen to generate the mixed waves at the sum and difference frequency, as shown in Fig. 11 (d). The amplitude of the mixed waves at the difference frequency are slightly larger than the amplitude of the mixed P wave at the sum frequency. The generation of the mixed S wave at the difference frequency ensures the mode selectivity and may also be an interesting candidate for wave mixing with only single face access. The results for the interaction between SV and P waves are presented in Fig. 11 (e) and (f). It is again possible to generate both mixed P and SV waves at the sum or difference frequency depending on the configuration. Fig. 11 (f) shows an example where the scattering of the mixed waves at the sum frequency occurs for scattering angles above 90
whereas the scattering of the mixed SV wave at the difference frequency occurs for angle below 90
. It is noted that mixed SV wave at the difference frequency has a larger amplitude than any of the other mixed waves obtained in Fig. 11(aee). However, the direction of propagation of the mixed wave and the required angles of incidence may preclude the interest of this particular configuration for practical use, and more adequate parameters are needed to exploit the mixed SV wave at the difference frequency. These numerical results validate the analytical solutions regarding the existence of mixed waves, and their directions of propagation and amplitudes. These results also confirm the differences with bulk mixing noted in Section 5. Wave mixing at a contact interface offers more flexibility regarding the choice of the incident wave modes, frequencies and angles of incidence than is the case for bulk wave mixing. In practice, where the single face access is often a requirement for NDE applications, the use of two incident S waves or two incident P waves for the generation of the mixed P at the sum frequency and a SV wave at the difference respectively, seem to be the most interesting configurations considering the directions of propagation and the required angles of incidence. Returning to the particular case SV/SV-1, consider now the response for different compression loads (and hence different T linear stiffness K N 0 and K 0 computed from Eqs. (5a) and (6)) and different incident-wave amplitudes, but keeping the amplitudes of both waves equal. Fig. 12 (a) shows an excellent agreement between the analytical and numerical results for the evolution of the mixed-wave amplitudes with the applied contact pressure, for an incident amplitude of 10 nm, whereas Fig. 12 (b) shows that for a fixed contact stress p0 ¼ 0:5MPa, the mixed-wave amplitudes are correctly predicted by the analytical solution up to incident wave amplitude of 12.5 nm for the considered case. The discrepancy observed for the higher incident amplitudes can be attributed to a limitation of. The perturbation approach as an asymptotic approximation which is appropriate in the limit of weak nonlinearity or small perturbation (i.e. small amplitude). However, for stronger nonlinearity, it is still expected that the analytical solutions retain the correct trends for the dependence on input parameters such as wave mode and interaction angle.

Fig. 12. Evolution of the mixed P and SV wave amplitude at the sum frequency normalised by the incident wave amplitude, as a function of the contact pressure (a) for an incident amplitude of 10 nm, and as a function of the incident amplitude (b) for p0 ¼ 0:5 MPa. The solid lines correspond to the analytical predictions whereas the discrete markers correspond to the FE results. (Color online). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

20

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

7. Conclusion An analytical model for wave mixing at a contact interface has been presented, based on a perturbation analysis for a nonlinear traction law representative of rough-surface contact. This has been shown to lead to several new fundamental insights into the characteristics of interface mixing, which differ significantly from those of bulk wave mixing due to material nonlinearity. In both cases, the nonlinear response is assumed to be a small correction to the linear response. However, the response for an interface bridged by linear springs involves a characteristic frequency, related to the linear spring coefficient for the interface traction law (which can in turn depend on the contact stress), whereas the bulk linear response involves no characteristic frequency. This means that experimental characterisation of interface mixing due to rough-surface contact should first seek to identify that characteristic frequency, which is most readily achieved by investigating the reflection and/or transmission coefficients for normal wave incidence. A systematic characterisation of interface mixing should be presented in terms of the normalised frequency f =f0 : This approach has been used here to show the variation in the optimum interaction angle and the normalised amplitude of the mixed wave for the case of symmetrical SV/SV incidence. However, the characteristics of interface mixing predicted by the analytical model have been confirmed by a computational model which also served to identify promising combinations of incident waves for applications in NDE. It is hoped that the present work will provide a useful basis to guide the design and interpretation of future experimental investigations. Declaration of competing interest None. Acknowledgments This work was supported by an Australian Research Council (ARC) Discovery Grant (DP180102658). Appendix A. Detailed derivations of analytical results The objective in this Appendix is to indicate the main steps in deriving the analytical results presented in Section 4. The interface traction law is assumed to be given by ∞ X

sxy ¼ syy ¼

εm K Tm Dumþ1 ; x

m¼0 ∞ X

(A.1) mþ1 εm K N ; m Duy

m¼0

where ε denotes a small non-dimensional parameter, reflecting an assumption that the nonlinear response can be regarded as a small perturbation relative to the leading order linear term, and that the solution can therefore be obtained by the method of successive approximations [35]. Thus, the relative displacements are also expressed as power series,

Dux ¼ Duy ¼

∞ X n¼0 ∞ X

ðnÞ

εn Dux ; (A.2) n

ε

DuyðnÞ ;

n¼0 ðnÞ

ðnÞ

where Dux and Duy denote the n-th order normal and tangential relative displacements. By substituting Eq. (A.2) into Eq. (A.1) and collecting terms according to their power of ε , we obtain the following governing equations at the zeroth and first order, respectively,









Inc sSxy Duxð0Þ  K T0 Duð0Þ x ¼ sxy ;

(A.3)

sSyy Duyð0Þ  K N0 Duyð0Þ ¼ sInc yy ; where superscripts I and S identify quantities pertaining to the incident and scattered field, and









sSxy Duxð1Þ  K T0 Duð1Þ x ¼ 0;



2

sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 Duyð0Þ :

(A.4)

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

21

It can be seen that the solution of the zeroth-order equation (A.3), which correspond to linear springs across the interface, provide the source term for the first-order equation (A.4). Accordingly, we shall first record the results for the linear-springs response for various combinations of incident waves. The principle of superposition applies for linear springs, and therefore the response to two incident waves is simply the sum of the responses to each wave separately. Therefore, it will be sufficient to record the results for incident SV and P waves separately. It was shown in Ref. [29] that analytical solutions can be derived for the linear-springs responses characterised by Eq. (A.3) by expressing the scattered field in terms of Mode I and Mode II contributions, because this decomposition leads to a decoupling to the two equations in (A.3). In particular, for an incident SV wave A shown in Fig. 2 (a,b) and specified by Eq. (17a), the amplitudes of the Mode I and II relative displacements are given by Eq. (17b), whereas for SV wave B, with the polarisation vector as shown in Fig. 2 (a,c), and specified by Eq. (18a), the relative displacements are given by Eq. (18b). These relative displacements can then be substituted into Eqs. (12) and (14) to derive the scattered fields, from which the stress components can also be obtained. Next, for an incident P wave shown in Fig. 4 (a) specified by Eq. (26a), the Mode I and II relative displacements are found to be [29]. ð0Þ

UA ¼

ð0Þ VA

kAL

4aA bA k2T xA  . AA ; FðxA Þ þ 2ik2T bA K T0 m

  2 2aA k2T 2xA  k2T . AA ; ¼  kAL FðxA Þ þ 2ik2T aA K N 0 m

(A.5)

with FðxÞ defined in Eq. (17b), whereas for the P wave specified by Eq. (26b), the relative displacements are found to be ð0Þ

UB ¼

ð0Þ VB

kBL



4aB bB k2T xB

. AB ; FðxB Þ þ 2ik2T bB K T0 m

  2 2aB k2T 2xB  k2T . AB : ¼  kBL FðxB Þ þ 2ik2T aB K N 0 m

(A.6)

ð0Þ

ð0Þ

Recalling that xB is negative, this leads to a tangential gap opening U B that has an opposite sign in comparison with U A ð0Þ VB

ð0Þ VA

obtained in Eq. (A.5), whereas the normal gap opening and have the same sign as one may expect. These results for single-wave incidence can now be combined to provide the relative displacements for various combinations of incident waves to obtain the source term in Eq. (A.4). At this point, some caution is required when using the complex exponential representation for the incident and scattered waves. It must be recalled that Eq. (A.4) involves physical ð0Þ

quantities, and therefore Duy must be interpreted as the real part of the solution obtained from the complex representation ð0Þ 2

prior to evaluating the source term ðDuy Þ : Accordingly, the correct source term is that shown in Eq. (20), with the relevant amplitudes and phase constants being determined for the relevant incident wave A or B in isolation, as summarised above. Thus, for the case of SV/SV incidence, for example, the relevant equations for determining the first-order correction corresponding to the sum frequency source term are found to be









sSxy Duxð1Þ  K T0 Duxð1Þ ¼ 0 h

sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 V A0 V

B 0

i  cos ðxA þ xB Þx  ðuA þ uB Þt þ zA0 þ zB0 :

(A.7)

However, now that the source term has been correctly identified, it is advantageous to exploit the fact that Eq. (A.7) is linear and to resort again to a complex exponential formulation, i.e. to seek the solution of Eq. (A.7) as the real part of the solution of









sSxy Duxð1Þ  K T0 Duxð1Þ ¼ 0; sSyy Duyð1Þ  K N0 Duyð1Þ ¼ K N1 V A0 V B0 eiððxA þxB ÞxðuA þuB Þtþz0 þz0 Þ : A

B

(A.8)

This equation now has the same form as the zeroth-order equation for which an explicit solution can be derived [29], which leads to the following complex amplitudes for the first-order correction

22

L.R.F. Rose et al. / Journal of Sound and Vibration 468 (2020) 115078

ð1Þ

U ABþ ¼ 0; ð1Þ V ABþ

¼

. 2iaABþ k2T K N 1 m

. V A0 V B0 :

(A.9)

FðxABþ Þ þ 2iaABþ k2T K N 0 m

The same approach can be used for the other combinations of incident waves. The approach can also be extended to the case of nonlinearity for the shear response of the interface, as was done in Ref. [29] for single-wave incidence. References [1] A.J. Croxford, P.D. Wilcox, B.W. Drinkwater, P.B. Nagy, The use of non-collinear mixing for nonlinear ultrasonic detection of plasticity and fatigue, J. Acoust. Soc. Am. 126 (2009) 117e122. [2] Y.H. Pao, W. Sachse, H. Fukuoka, Acoustoelasticity and ultrasonic measurement of residual stress, in: W.P. Mason, R.N. Thurston (Eds.), Physical Acoustics, vol. XVII, Academic Press, New York, 1984. [3] N. Gandhi, J.E. Michaels, S.J. Lee, Acoustoelastic Lamb wave propagation in biaxially stressed plates, J. Acoust. Soc. Am. 132 (2012) 1284e1293. [4] J.H. Cantrell, K. Salama, Acoustoelastic characterisation of materials, Int. Mater. Rev. 36 (1991) 125e145. [5] P.B. Nagy, Fatigue damage assessment by nonlinear ultrasonic materials characterization, Ultrasonics 36 (1998) 375e381. [6] J.H. Cantrell, W.T. Yost, Nonlinear ultrasonic characterization of fatigue microstructures, Int. J. Fatigue 23 (2001) 487e490. [7] J.H. Cantrell, Fundamentals and applications of nonlinear ultrasonic nondestructive evaluation, in: T. Kundu (Ed.), Ultrasonic Nondestructive Evaluation, CRC Press, Boca Raton, 2003. Chap. 6. [8] G.L. Jones, D.R. Kobett, Interaction of elastic waves in an isotropic solid, J. Acoust. Soc. Am. 35 (1963) 5e10. [9] M. Hasanian, C.J. Lissenden, Second order harmonic guided wave mutual interactions in plate: vector analysis, numerical simulation, and experimental results, J. Appl. Phys. 122 (2017), 084901. [10] M. Hasanian, C.J. Lissenden, Second order harmonic guided wave mutual interactions in plate: arbitrary angles, internal resonance, and finite interaction region, J. Appl. Phys. 124 (2018) 164904. [11] Y. Ishii, K. Kiraoka, T. Adachi, Finite-element analysis of non-collinear mixing of two lowest-order antisymmetric Rayleigh-Lamb waves, J. Acoust. Soc. Am. 144 (2018) 53e68. [12] H. Cho, M. Hasanian, S. Shan, C.J. Lissenden, Nonlinear guided wave technique for localized damage detection in plates with surface-bonded sensors to receive Lamb waves generated by shear-horizontal wave mixing, NDT E Int. 102 (2019) 35e46. [13] A. Demcenko, R. Akkerman, P.B. Nagy, R. Loendersloot, Non-collinear wave mixing for non-linear ultrasonic detection of physical ageing in PVC, NDT E Int. 49 (2012) 34e39. [14] A. Dem cenko, V. Koissin, V.A. Korneev, Noncollinear wave mixing for measurement of dynamic processes in polymers: physical ageing in thermoplastics and epoxy cure, Ultrasonics 54 (2014) 684e693. [15] M.E. McGovern, W.G. Buttlar, H. Reis, Characterisation of oxidative ageing in asphalt concrete using a non-collinear ultrasonic wave mixing approach, Insight Non-Destructive Test. Cond. Monit. 56 (2014) 367e374. [16] A. Dem cenko, L. Mainini, V.A. Korneev, A study of the noncollinear ultrasonic-wave-mixing technique under imperfect resonance conditions, Ultrasonics 57 (2015) 179e189. [17] J. Alston, A. Croxford, J. Potter, P. Blanloeuil, Nonlinear non-collinear ultrasonic detection and characterisation of kissing bonds, NDT E Int. 99 (2018) 105e116. [18] P. Blanloeuil, A. Meziane, C. Bacon, 2D finite element modeling of the non-collinear mixing method for detection and characterization of closed cracks, NDT E Int. 76 (2015) 43e51. [19] J.-M. Baik, R.B. Thompson, Ultrasonic scattering from imperfect interfaces: a quasi-static model, J. Nondestruct. Eval. 4 (1984) 177e196. [20] S.I. Rokhlin, Y.J. Wang, Analysis of boundary conditions for elastic wave interaction with an interface between two solids, J. Acoust. Soc. Am. 89 (1991) 503e515. [21] P. Nagy, Ultrasonic classification of imperfect interfaces, J. Nondestruct. Eval. 11 (1992) 127e139. [22] P.P. Delsanto, M. Scalerandi, A spring model for the simulation of the propagation of ultrasonic pulses through imperfect contact interfaces, J. Acoust. Soc. Am. 104 (1998) 2584e2591. [23] A. Baltazar, S.I. Rokhlin, C. Pecorari, On the relationship between ultrasonic and micromechanical properties of contacting rough surfaces, J. Mech. Phys. Solids 50 (2002) 1397e1416. [24] C. Pecorari, Nonlinear interaction of plane ultrasonic waves with an interface between rough surfaces in contact, J. Acoust. Soc. Am. 113 (2003) 3065e3072. [25] S. Biwa, S. Nakajima, N. Ohno, On the acoustic nonlinearity of solid-solid contact with pressure-dependent interface stiffness, J. Appl. Mech. 71 (2004) 508e515. [26] S. Biwa, S. Hiraiwa, E. Matsumoto, Stiffness evaluation of contacting surfaces by bulk and interface waves, Ultrasonics 47 (2007) 123e129. [27] C. Pecorari, Spring boundary model for a partially closed crack, Int. J. Eng. Sci. 46 (2008) 182e188. [28] Z. An, X. Wang, M. Deng, J. Mao, M. Li, A nonlinear spring model for an interface between two solids, Wave Motion 50 (2013) 295e309. [29] P. Blanloeuil, L.R.F. Rose, M. Veidt, C.H. Wang, Analytical and numerical modelling of wave scattering by a linear and nonlinear contact interface, J. Sound Vib. 456 (2019) 431e453. [30] T. Nam, T. Lee, C. Kim, K.-Y. Jhang, N. Kim, Harmonic generation of an obliquely incident ultrasonic wave in solid–solid contact interfaces, Ultrasonics 52 (2012) 778e783. [31] Z. Zhang, P.B. Nagy, W. Hassan, Analytical and numerical modeling of non-collinear shear wave mixing at an imperfect interface, Ultrasonics 65 (2016) 165e176. [32] B.W. Drinkwater, R.S. Dwyer-Joyce, P. Cawley, A study of the interaction between ultrasound and a partially contacting solid-solid interface, Proc. R. Soc. Lond. A 452 (1996) 2613e2628. [33] J. Achenbach, Wave Propagation in Elastic Solids, Elsevier, 1973. [34] J.L. Rose, Ultrasonic Guided Waves in Solid Media, Cambridge University Press, 1999. [35] A.H. Nayfeh, Perturbation Methods, John Wiley & Sons, 2008. [36] M. Smith, Abaqus User Subroutines Reference Guide, Simulia, Providence RI, 2016. [37] L.L. Thompson, A review of finite-element methods for time-harmonic acoustics, J. Acoust. Soc. Am. 119 (2006) 1315e1330.