Revisiting imperfect interface laws for two-dimensional elastodynamics

We study the interaction of in-plane elastic waves with imperfect interfaces composed of a periodic array of voids or cracks. An effective model is derived from high-order asymptotic analysis based on two-scale homogenization and matched asymptotic technique. In two-dimensional elasticity, we obtain jump conditions set on the in-plane displacements and normal stresses; the jumps involve in addition effective parameters provided by static, elementary problems being the equivalents of the cell problems in classical two-scale homogenization. The derivation of the model is conducted in the transient regime and its stability is guarantied by the positiveness of the effective interfacial energy. Spring models are envisioned as particular cases. It is shown that massless-spring models are recovered in the limit of small void thicknesses and collinear cracks. By contrast, the use of mass-spring model is justified at normal incidence, otherwise unjustified. We provide quantitative validations of our model and comparison with spring models by means of comparison with direct numerical calculations in the harmonic regime.


Introduction
Surfaces of separations between two elastic solids can impact significantly the propagation of waves. In 1967, Jones & Whittier [1] inspected the ability of thin elastic bonds to support guided waves. Considering inertialess layers with increasing stiffness, they exhibited a family of guided waves with limiting cases being the Rayleigh waves and the Stoneley waves. Nine years later, Murty [2] considered the scattering by what he called loosely bonded interfaces; having in mind situations as different as aligned cracks or thin viscous layers, he postulated the simplest contact law able to recover a perfect interface and a fully debonded interface as limiting cases. Assuming that the normal stress σ remains unconditionally continuous across the imperfect interface, his model coincides with the massless-spring model of [1]: in a scalar case, it reads σ = κ [[u]] with [[u]] = (u + −u − ) the jump of the displacement u across the interface and κ a parameter termed bonding parameter or stiffness. Later on, Baik and Thompson have enriched the model by including inertial terms; the resulting massspring model imposes displacements and normal stresses discontinuities [3]. A review of the different models has been analysed by Martin [4] who showed that some of them suffer from non-uniqueness.
The case of thin homogeneous layers is of interest for practical applications; besides the exact solution is available. In a series of papers, Rokhlin and co-workers have developed accurate transmission conditions and inspected the validity of spring models [5][6][7]. This has revealed the limited range of validity of the massless-spring model; in particular, neglecting the inertial terms appears to be unjustified [8][9][10]. The case of thin layer interfaces containing cracks, micropores, voids or faults has also received attention in particular for applications in seismology and in engineering for non-destructive testing. Their academic study started in the 1980s with the works of Achenbach and his co-workers who investigated the scattering properties in many situations including equally spaced collinear cracks [11,12], inclined cracks [13], spherical cavities [14] as well as randomly distributed cracks [15]; a review is presented in [16]. Exact solutions are in general not available and numerical solutions have been sought by means of multiple scattering theory and boundary integral equation method. For zero thickness cracks, the validity of the spring models has been demonstrated [11,12,[17][18][19]. However, a quantitative disagreement for crack like defects, or diffusion bonds, with small but finite thickness has been reported in [19].
The present works aims to derive a model on the basis of rational approximations for defects with non-zero thickness and to assess if the use of spring models is legitimate in some limit. To that aim, we extend the results of Marigo and co-workers [20,21] developed in a static context; the idea is to combine two-scale homogenization theory to treat the periodicity of the inhomogeneities and matched asymptotic technique to deal with the thin layer. In §2, we introduce the full-scale problem of a periodic array of voids and the effective problem whose derivation is detailed in §3. Simplifications of our model for waves at normal incidence and for vanishing thickness of the cracks are envisioned, see §3e. It is shown that for very thin voids, which means with a thickness much smaller than the array spacing, a massless-spring model is obtained which involves two different, tangential and normal, stiffnesses. The model further simplifies for cracks with zero thickness, resulting in a single stiffness whose form is explicit, and conform with the model of Angel & Achenbach [11]. Amusingly, and contrarily to the case of a thin homogeneous layer, keeping in addition an inertial term is unjustified except at normal incidence. Hence the mass-spring model does not apply in general to an array of thin voids (and it does not apply to cracks which are inertialess by construction). Our analysis is set in the time domain which allows us to analyse the energetic properties of the effective problem, see §4. The jump conditions produce a non classical contribution to the flux of the elastic Poynting vector. It is shown that this contribution reads as the time derivative of a positive effective energy which guaranties the stability in the transient regime. Section 5 illustrates the effectiveness of the our model by comparison with full-scale simulations of the two-dimensional elastic problem [22]. (Closed-forms of the scattering coefficients and of the fields are provided.) We take the opportunity in this section to discuss the effectiveness of other formulations (zero-thickness formulation and massless-spring models). We provide in §6 concluding remarks and extensions of the study.

The actual problem and the effective problem
We consider a two-dimensional linear elasticity problem where a periodic array of voids, identical but of arbitrary shape, is embedded in a homogeneous matrix of Lamé coefficients (λ, μ) and mass density ρ, figure 1. With x = (x 1 , x 2 ) the coordinate system, the voids are evenly distributed along x 2 and within x 1 ∈ (0, e) with a spacing h and e is comparable to h. We consider the propagation of elastic waves in the time domain, with t the time. Denoting (u, σ , ε x ) the displacement vector, the stress and the strain tensors, the governing equations in the actual problem read where dot means time derivative and 1 stands for the identity matrix. The above problem is complemented with stress-free conditions σ n = 0 on the boundary of the voids and radiation conditions of the Sommerfeld type. (At this stage, we do not need to specify them.) In the following, with (e 1 , e 2 ) the unit vectors of the basis, we denote u i = u · e i , i = 1, 2 and σ ij = σ · e i ⊗ e j , (i, j) ∈ {1, 2} 2 . We consider a source which imposes a typical, or maximum, frequency ω such that k T h 1 with k T = ω/c T , c T = ω √ ρ/μ. In this subwavelength regime, our study aims to simplify the actual problem. Specifically, the separation of the length scales between the wavelength and the dimensions (e, h) of the array is exploited in an asymptotic analysis detailed in the forthcoming section. As a result, the array of voids is replaced by jump conditions across an effective interface with boundaries at x 1 = 0 and x 1 = e. (It is worth noticing that the region x 1 ∈ (0, e) is not interrogated; only the values of the fields at x 1 = 0 and e are concerned.) These jumps tell us that the in-plane displacements and the normal stress satisfy with (i, j, k) ∈ {1, 2} 3 . In the above transmission conditions, for any scalar field f (x, t), we have defined Eventually for voids with symmetric shape with respect to x 1 or to x 2 , only four effective parameters are involved.

Asymptotic analysis
As previously said, the asymptotic procedure is conducted owing to the separation of the length scales between the typical wavelength imposed by the source and the dimensions (e, h) of the array. We define η = k T h 1 the small parameter measuring the separation of the length scales and k T e = O(η). In this section, we use the dimensionless coordinate x = (x 1 , x 2 ) with x = k T x. Next, for simplicity, we preserve the formulation (2.1) by choosing t = k T t and u(x, t) = k T u(x, t), σ (x, t) = σ (x, t). Doing so, the rescaled formulation of (2.1) reads where dot means now derivative with respect to t. We also define the microscopic coordinate at the scale of the voids where V is the rescaled region of a single void (figure 2). We shall now define different expansions of the solution depending on either we are far or close to the array. Far from the array, the spatial variations of the fields are slow, being attributable to the wave propagation, hence to the largescale 1/k T . Accordingly, we postulate asymptotic expansions for the outer solution of the form In contrast close to the array, the fields have slow variations along x 2 due to phase variations of the waves along the array but they have also fast variations due to the evanescent field excited in the vicinity of the voids. (The evanescent field has spatial variations on the scale of h.) Hence, we postulate two-scale asymptotic expansions for the inner solution of the form The fields (v i , τ i ) are 1-periodic functions w.r.t. y 2 which is formally equivalent to Bloch-Floquet conditions in the limit of infinite wavelengths. The resulting inner problem is set in the elementary cell Y ( figure 2). For convenience, we also define Y m = {y 1 ∈ (−y m 1 , y m 1 ), y 2 ∈ (− 1 2 , 1 2 )}/V, and Y m → Y when y m 1 → +∞. It is worth noting that the boundary condition τ i n = 0 on the boundary of the voids applies to the inner solution (3.3); however there are missing boundary conditions when y 1 → ±∞. Reversely, the radiation conditions for x 1 → ±∞ apply to the outer solution    Doing so, we get at the zero order (3.4) and at the first order and It remains to establish the hierarchy of problems issued from (3.1) and satisfied by the terms of the outer and inner expansions. Given the form of the expansions (3.3), each term satisfies (3.1) namely In contrast for the inner terms, the hierarchy of problems is obtained inserting the expansions (3.3) in (3.1) and using the differential operator ∇ → (e 2 ∂ ∂x 2 + 1 η ∇ y ) (according to the two scales x 2 and y). Collecting the contributions in η i , we get the constitutive behaviours (C) and the equilibrium equations (E). At the zero and first orders, they read and For sake of simplicity, we have defined the following operators: Eventually, the stress-free condition on the boundary of the void applies at each order In the following, we shall solve (3.7)-(3.9) at the first two orders.

(a) Zero-order resolution
Inverting the first relation of (C) in (3.7), we get that ε y (v 0 ) = 0. Hence v 0 is associated to a rigid body motion which is further reduced to a translation due to the periodicity of v 0 with respect to y 2 . Accordingly, we denote In virtue of the matching conditions for v 0 in (3.4), we deduce that where we have used the periodicity of τ 0 with respect to y 2 and the stress-free condition τ 0 n |∂V = 0 from (3.9). Passing to the limit y m 1 → +∞ and making use of the matching condition in (3.4), we eventually get It is worth noting that, from (3.10) and (3.11), ∇ x u 0 is continuous at x 1 = 0 (in the following we denote ∇ x u 0 (0, x 2 , t) the limit). At the dominant order, the array of voids is invisible and we have to go to the first order to capture its effect.

(b) First-order resolution
At this order, we have to solve a problem set on (τ 0 , v 1 ) with respect to the spatial variable y. From (3.7)-(3.8) along with the matching conditions on τ 0 in (3.4), this problem reads In the definition of τ 0 in (3.7), we have used that v 0 = u 0 (0, x 2 , t), from (3.10). Next the limits of ε y (v 1 ) for y 1 → ±∞ follow from (3.4) and using σ 0 = A(ε x 1 (u 0 ) + ε x 2 (u 0 )). The problem (3.12) appears to be linear with respect to the two macroscopic loadings (ε x 1 (u 0 ), ε x 2 (u 0 )). Hence the solution can be expressed as follow where (i, j) ∈ {1, 2} 2 and where the elementary vectors Q ij and associated tensors T ij satisfy the elementary problems where e i ⊗ s e j = 1/2(e i ⊗ e j + e j ⊗ e i ). From (3.14), we see that up to a constant, Q 12 = Q 21 where the constant vectors B ij are effective parameters independent of the macroscopic loading. Now, using the matching conditions (3.5) for the displacement u 1 , it is easy to see from (3.13) that from which we deduce that 3 . Yet we still have to derive the first-order condition for the normal stress. We start by expressing τ 0 as a function of the elementary solutions Q ij . From (3.12)-(3.13), and owing We use the above expression in the second equation of (E), (3.8), along with v 0 = u 0 (0, x 2 , t) from (3.10). We also use from (3.6) that 0 = div x 1 σ 0 + div x 2 σ 0 − ρü 0 . Doing so, we obtain div y τ 1 Now, integrating (3.18) over the cell Y and accounting for the matchings on τ 1 in (3.5) leads to where the vectors C ij are defined by and eventually, for k = 1, 2, Following [20,21], up-to-first order conditions are constructed by recombining the zero and firstorder conditions on (u 0 + ηu 1 ) from (3.10) and (3.11) and on (σ 0 + ησ 1 ) from (3.17) and (3.20). It is worth noticing that these conditions implicitly involve the length h through η = kh. As the array thickness e = O(h) is of the same order of magnitude, it makes sense to account for its contribution too. Introduced in [23], enlarged formulations of the jump conditions has been shown to be stable, see also [24]. They are obtained by Taylor expansions of u 0 (0 + , x 2 , t) and σ 0 (0 + , x 2 , t and with the notations given by in (2.3). Combining the above equations and coming back to the real coordinate and time (x, t) leads to the final conditions in (2.2) which apply to the unique displacement (d) Properties of the effective coefficients and Proof. The symmetries in (3.21) are straightforward since Q ij = Q ji from (3.14), To establish (3.22), we integrate over Y the equilibrium equation div y T ij = 0 in (3.14) after multiplication by y 2 e 1 which is periodic along y 2 . Integrating by part Y m dy y 2 e 1 · div y T ij = 0 and accounting for the periodicity of T ij , the nullity of T ij n on ∂V and the limits lim which provides the values of C ij 1 in (3.22). Next, to show (3.23), we need to introduce D ij = Y dy Aε y (Q ij ). Repeating the procedure after multiplication by y 1 e 1 provides the relations We now integrate div y T ij = 0 after multiplication by Q kl , specifically and we know that A(e i ⊗ s e j )e 1 = λtr(e i ⊗ s e j )e 1 + μ(δ 1i e j + δ 1j e i ).
Eventually, considering q 1112 , q 1122 and q 1222 , it is easy to check that with respect to y 2 . As a result the component Q 22 2 in the same direction is also odd which further implies that its limit B 22 2 = 0. For the same reasons, the elementary problem on Q 12 , associated with a loading e 1 ⊗ s e 2 = ε y (y 2 e 1 ) (and y 2 e 1 is odd with respect to y 2 ), produces an odd Q 12 1 , whence B 12 1 = 0. Suppose now that the defects are symmetric with respect to y 1 (invariance y 1 → −y 1 ). The loading y 2 e 2 associated with the problem on Q 22 is even with respect to y 1 , hence Q 22 2 is also even.
We find the enlarged version of the inertial terms introduced by Baik & Thomson [3] in the jumps of the normal stresses. The jumps of the displacements involve four compliances from which the four stiffnesses can be deduced by inversion. This mass-spring model is as accurate as (2.2) but restricted to normal incidence. The above massless spring model is an approximation of our model up to O(e). It involves two, normal and tangential, stiffnesses which differ as B 11 1 /(λ + 2μ) and B 12 2 /μ do; these later still depend on e as they depend although weakly to the geometry of the thin voids.
-Aligned cracks (e = 0). This case is the limit of the previous one. We have checked that for vanishing e, we have κ 1 = κ 2 = κ with and we recover the model of Angel & Achenbach [11] with a single, explicit, stiffness.  1)) byu and the constitutive behaviour (second equation in (2.1)) by ε x (u), and we sum up. Doing so and after an integration by part, we get

Energetic properties of the effective model
E Ω is the classical elastic energy and Φ ∂Ω is the flux of the elastic Poynting vector (σu). Usually, the fieldsu and σ n are continuous or they vanish on ∂Ω. Here, the imperfect interface produces a discontinuity of these fields, hence it makes a contribution Φ Γ to Φ ∂Ω which reads
The positivity of E Γ (t) guarantees the stability of the effective problem, which basically means that u and σ are bounded [23]. In practice, E Γ > 0 prevents from numerical instabilities which happen when an increase in time of E Ω (t) → +∞ is compensated by E Γ → −∞, see an illustration in [25].

A scattering problem, validation and discussion
In this section, we consider the harmonic regime with time dependence e −iωt ; we shall inspect the validity of our effective model by means of comparison with direct numerics based on a multimodal method for two-dimensional elasticity [22]. We restrict ourselves to rectangular voids, whose symmetry makes that, from (3.22) We have in addition from (3.23) that λB 11 1 − C 11 2 = (λ+2μ)B 22 1 − λϕe.

(a) Scattering coefficients and energy fluxes
In two dimensions, the wavefield is defined in terms of the elastic scalar potentials Φ and Ψ such that u = ∇Φ + ∇ × (Ψ e 3 ). The incident wave in the matrix derives from the potentials with (α L , β) = k L (cos θ L , sin θ L ), (α T , β) = k T (cos θ T , sin θ T ), (5.2) and The domain is unbounded and the radiation conditions are accounted for by imposing that the scattered field is outgoing; in the effective problem, the solution simply reads  x 2 (where means real part) and the results holds by summing the two equations for an incident flux For an incident longitudinal wave (I=L) or an incident transverse wave (I=T), we define the normalized energy fluxes in reflection φ r OI and in transmission φ t OI with O=L,T by and the conservation of the fluxes is ensured by

(b) Results
The material and geometrical parameters are ρ = μ = λ = 1, h = 1 and we shall consider different (e, ϕ). We shall inspect the effectiveness of our model in the range k T h ∈ (0, 3) hence well beyond the expected range of validity of homogenization (which assumes η = k T h 1).
In the first exemple for square voids e = ϕ = 0.5, we shall take the opportunity to illustrate the consequence of the enlarged formulation (5.1) in which the jumps and means are expressed for an interface with boundaries at x 1 = 0 and e, (2.3). More usually, the zero-thickness formulation, that is provided directly by the asymptotic analysis, is used for an interface with boundaries at x 1 = 0 − and 0 + . (It involves jumps (f (0 + , x 2 ) − f (0 − , x 2 )) and mean values 1/2(f (0 + , x 1 ) + f (0 − , x 2 )).) This formulation does not account for the Taylor expansions used in §3c. Hence it is obtained by simply removing in (5.1) the contributions due to Taylor expansion, e.g. he∂ Then, we shall move on the case of aligned thin voids, with e h whose limiting case are cracks, i.e. e = 0. This allows us to inspect in addition the effectiveness of the massless-spring models

(i) Scattering by square voids
Here we consider square voids with ϕ = e = 0.5, hence an array with significant thickness. We start by reporting in figure 4 the displacement u 1 (x 1 , x 2 ) for an incident wave of the form (5.2) with k T h = 2, (A L = 3, A T = 1), (θ T = 30 • , θ L = 60 • ) producing k L sin θ L = k T sin θ T . The qualitative agreement between the field computed numerically and that given by (5.3) is excellent. In the actual problem, the near field is composed of many evanescent modes while in the homogenized problem, it is accounted for by the jump conditions. (Hence, in the close vicinity of the array, the variation of the evanescent field is not reproduced by construction, see insets of figure 4.) More quantitatively, we report in figure 5 the normalized fluxes against k T h obtained numerically (plain lines) and from our model (black dashed lines). Our model, in its enlarged formulation, has an excellent accuracy up to k T h ∼ 1 and a reasonably good one up to k T h ∼ 2. By contrast, the validity of the model in the zero-thickness formulation (grey dashed lines) is severely limited to the low frequencies. It is a good news that our model which has good energetic properties has in addition a better effectiveness.    and F cos θ L /d F = G cos θ T /d G = 1/i; we recover a perfect transmission without L-T polarization conversion.
We have thus, in addition to our model (5.1), an approximate massless-spring formulation for e 1, (3.28) involving two stiffnesses, and its limit e = 0 resulting in a single, explicit, stiffness literature for waves at normal incidence, hence θ T = θ L = 0 in (5.10). (This case does not give rise to mode conversion.) We start with this case in figure 8 and move on to oblique incidence in figure 9 (mode conversions occur). The plots show the results for four thicknesses from very thin e = 0.01 to relatively thin e = 0.3. It turns out that the thickness has a strong influence on the variations of the scattering coefficients, which is accurately captured by our model (dashed black lines) up to k T h ∼ 2 whatever the array thinness. Next, the massless-spring model with two stiffnesses (3.28) (dashed grey lines) has roughly the same accuracy for e = 0.01 and it rapidly fails in accuracy for thicker voids. Eventually the limit of a single stiffness (3.29) (providing the curve in dotted grey lines) has a limited range of validity even for the thinner array (this is more noticeable at normal incidence).

Conclusion
In this paper, we have derived effective jump conditions for an imperfect interface composed of an array of defects, e.g. voids or cracks. The conditions are expressed across an interface having the same thickness as the actual ones (and whose interior is not resolved); we have shown that this enlarged formulation enjoys good energetic property, as the interface has a positive energetic contribution. It is worth noticing that this property prevents from possible severe instabilities of any explicit time-discretization numerical scheme. Besides, by comparison with direct numerical calculations, we have shown that the enlarged formulation is more efficient than the zero-thickness formulation expressed across a zero-thickness interface. In the reported examples, we have observed that it remains efficient for any incidence and far beyond the long wavelength limit, usable almost up to the appearance of the higher modes when k T > β + 2nπ/h, n integer. Spring models have been recovered: as a particular case, the mass-spring model is valid at normal incidence, and as limiting cases spring-models approximate our model for vanishing void thickness up to O(e/h).
Some extensions of the present study are straightforward, others less. Technically, the extension to three-dimensional elasticity is cumbersome but it does not present additional difficulties. The same holds if the voids are replaced by elastic inclusions with material parameters comparable to those of the matrix. The case of highly contrasted inclusions able to produce internal, low-frequency, resonances is more challenging. Already considered in the scalar, antiplane, case by Pham et al. [26] and Touboul et al. [27], it is difficult to anticipate how the effect of these resonances translates in two-or three-dimensional elasticity. From a practical point of view, the study of guided waves by such imperfect interfaces and their link with Rayleigh, Stoneley and Love waves in realistic configurations is of interest in a geophysical context. In this spirit, the effects of the viscous fluids within the imperfect interface or the influence of non-linear contact laws, are interesting paths to explore.
Data accessibility. The numerical simulations of the direct problem were obtained using the modal method published in Maurel & Pham [22].
Authors' contributions. The authors worked together on the formulation of the problem and its resolution. A.M. and K.P. conducted the numerical simulations. The content of the paper was discussed and approved by all authors.
Competing interests. We declare we have no competing interests. Funding. K.P. thanks the support of the Agence de l'Innovation de Défense (AID) from the Direction Générale

Appendix A. Positivity of the interfacial energy
We aim to show that the interfacial energy E P in (4.3) is definite positive. We divide the proof in two steps. In the first step, we show that It follows that E P is positive and satisfies the inequality in (4.5). The two steps require the use of an energy minimization principle, or variational principle, expressed in terms of displacements (step 1) or in terms of stresses (step 2).

(a) First step
To show (A 1) we shall prove that ∀(a, b) ∈ R 2 , (λ+2μ)B 11 1 a 2 + 2(λ+2μ)B 12 1 ab + μB 12 2 and At this stage, (a, b) ∈ R 2 are not specified and to show (A 3) we consider the displacement field Q ab = aQ 11 + bQ 12 and its associated stress tensor T ab . From (3.14), (Q ab , T ab ) solve the problem div y T ab = 0, T ab = A ε y (Q ab ) + ae 1 ⊗ s e 1 + be 1 ⊗ s e 2 and (Q ab , T ab ) y 2 -periodic, T ab n = 0 on ∂V, lim The strong formulation (A 5) can be derived from a weak formulation. This latter is obtained using that Y dy Q · div y T ab = 0 for any Q being an admissible test function, y 2 -periodic and whose gradient tends to 0 at when y 1 → ±∞. Integrating by parts leads to Y dy ε y (Q) · Aε y (Q ab ) + a Y dy (λ+2μ) where δ(Q k ) = ( lim Such weak formulation is itself derived from a variational principle i.e. an energy minimization principle. In that perspective we introduce the space D of admissible displacement fields defined by D = Q : Q y 2 -periodic, continuous, lim The potential energy E ab : D → R then reads The variational principle is given by and it can be shown that the first-order stationary condition of the potential energy E ab corresponds to the weak formulation (A 6). The potential energy of the solution is obtained taking Q = Q ab in (A 7) along with (A 6) and we find Eventually, applying the variational principle (A 8), it is straightforward to obtain (A 3). It is also straightforward to see that, with the choice of (a, b) in (A 4), hence a = ε 11 + λ λ+2μ ε 22 and b = 2ε 12 , (A 3) provides (A 1).

(b) Second step
We now move on the second step, where we aim to prove (A 2). For that we consider the specific displacement field Q * = a 1 Q 11 + a 2 Q 22 with (a 1 , a 2 ) given by and a 2 = λ+2μ 2(λ+μ) .