On radiating solitary waves in bi-layers with delamination and coupled Ostrovsky equations

We study the scattering of a long longitudinal radiating bulk strain solitary wave in the delaminated area of a two-layered elastic structure with soft (`imperfect') bonding between the layers within the scope of the coupled Boussinesq equations. The direct numerical modelling of this and similar problems is challenging and has natural limitations. We develop a semi-analytical approach, based on the use of several matched asymptotic multiple-scale expansions and averaging with respect to the fast space variable, leading to the coupled Ostrovsky equations in bonded regions and uncoupled Korteweg-de Vries equations in the delaminated region. We show that the semi-analytical approach agrees well with direct numerical simulations and use it to study the nonlinear dynamics and scattering of the radiating solitary wave in a wide range of bi-layers with delamination. The results indicate that radiating solitary waves could help us to control the integrity of layered structures with imperfect interfaces.


I. Introduction
][5][6] Initially developed as a purely analytical technique, 7,8 a) Accepted for publication in Chaos (2017): http://dx.doi.org/10.1063/1.4973854b) Corresponding author: K. R. Khusnutdinova.Tel: +44 (0)1509 228202.][10] An efficient IST-based numerical approach to solving the KdV equation was also developed. 11cently, the method has found a new application in our studies of the scattering of long longitudinal bulk strain solitons in a symmetric perfectly bonded layered bar with delamination. 12,135][16][17][18] The theoretical predictions 12 agreed well with experimental studies, 19 and were also confirmed by numerical simulations. 130][21] The exceptional stability of bulk strain solitons 20,22,23 makes them an attractive candidate for the introscopy of layered structures, in addition to the existing methods. 24,25e dynamical behaviour of layered structures depends not only on the properties of the bulk material, but also on the type of the bonding between the layers.In particular, if the materials of the layers have similar properties and the bonding between the layers is sufficiently soft ('imperfect bonding'), then the bulk strain soliton is replaced with a radiating solitary wave, a solitary wave with a co-propagating oscillatory tail. 26The radiating strain solitary wave has recently been observed in laboratory experiments. 24More generally, experimental studies of the excitation of the resonant radiation by localised waves have been a prominent theme in nonlinear optics and a number of other physical settings, see, for example, the reviews 27,28 and the references therein.
Indeed it was shown, within the framework of a complex lattice model, that long nonlinear longitudinal bulk strain waves in a bi-layer with a sufficiently soft bonding can be modelled with a system of coupled regularised Boussinesq (cRB) equations 26 (in non-dimensional and scaled form): Here, f and g denote the longitudinal strains in the layers, while the coefficients c, α, β, δ, γ are defined by the physical and geometrical parameters of the problem 26 (see Section II for details).
In the symmetric case (c = α = β = 1) system (1) admits the reduction g = f , where f satisfies the equation The Boussinesq equation ( 2) has particular solitary wave solutions: where v is the speed of the wave.However, in the system of cRB equations (1), when the characteristic speeds of the linear waves in the layers are close (i.e.][32] Figure 1 illustrates the pure solitary wave solutions of system (1) with δ = γ = 0 (for a fixed value of v), and subsequent evolution of this initial condition into a radiating solitary wave.4][35][36][37][38][39] The radiating solitary waves emerge due to a resonance between a solitary wave and some linear wave, which can be deduced from the analysis of the relevant linear dispersion relation.When considering the linear dispersion relation for the system (1), it is assumed that the coefficients in (1) are perturbed compared to the symmetric case, but remain positive. 26The dispersion relation has the form where k is the wavenumber and p is the phase speed.A typical plot of (3) is shown in Figure 2. A significant difference with the linear dispersion curve of the reduction (2) is the appearance of the second branch, remaining above the first branch for all k, and going to infinity as k → 0, while approaching zero as k → ∞.The pure solitary waves of the single Boussinesq equation (2) arise as a bifurcation from wavenumber k = 0 of the linear wave spectrum, when there is no possible resonance between the speed v of the solitary wave and the phase speed p of some linear wave.Radiating solitary waves arise in the case when there is a possible resonance for some finite non-zero value of k.For example, a possible resonance is shown in Figure 2 for v = p = 1.3.The aim of this paper is to study the scattering of this type of solitary wave in delaminated areas of imperfectly bonded layered structures (see Figure 3).We develop a semi-analytical approach which leads to coupled Ostrovsky equations in bonded regions of a bi-layer, and to uncoupled Korteweg-de Vries equations in the delaminated area.The Ostrovsky equation was originally derived to describe long surface and internal waves in a rotating ocean, 40,41 but recently it transpired that the equation, as well as the coupled Ostrovsky equations, can also describe nonlinear strain waves in layered elastic waveguides with soft interfaces. 29We also develop a direct numerical approach and verify that the semi-analytical method produces the correct results in the cases where we can use both methods.However the direct numerical simulations are expensive, therefore we then use our semi-analytical method to study the scattering of radiating solitary waves in a wide range of complex imperfectly bonded bi-layers with delamination, giving an elaborate description of the possible dynamical effects.The paper is organised as follows.In Section II we state the problem formulation for the generation and the scattering of a radiating solitary wave in an imperfectly bonded bi-layer with delamination.In Section III we develop a weakly nonlinear solution of this scattering problem and discuss the related semi-analytical approach.In Section IV, in a case study, we compare the results obtained using the semi-analytical approach with the results of direct numerical simulations, and show that the results are in good agreement.We then continue to use the semi-analytical approach to study the scattering of radiating solitary waves for a wide range of configurations of the complex structure.We summarise our findings and discuss the results in Section V.
3: Bi-layer with two homogeneous layers for x < x a , a bonded two-layered section for x a < x < x b and a delaminated section for x > x b .
We consider the generation and the scattering of a long radiating solitary wave in a two-layered imperfectly bonded bi-layer with delamination, shown in Figure 3.The model is inspired by the experimental setup. 24Two identical homogeneous layers (the section on the left) are 'glued' to a two-layered structure with soft bonding between its layers (in the middle), followed with a delaminated section (on the right).The materials in the bilayer are assumed to have close properties, leading to the generation of a radiating solitary wave in the bonded section. 26We study the scattering of this wave by the subsequent delaminated region.

+ 2u
(2) as well as some natural initial and boundary conditions, which will be imposed later.
To find the weakly nonlinear solution of the complicated scattering problem we consider the equations ( 11) - (13).
We use several asymptotic multiple-scale expansions, and develop a space-averaging method instead of the timeaveraging method used for the homogeneous initial-value problem. 29All functions present in our expansions and their derivatives are assumed to be bounded and sufficiently rapidly decaying at infinity (these assumptions agree with our numerical simulations).
In the regions where the behaviour is governed by uncoupled regularised Boussinesq equations, the previous work shows that to leading order the weakly nonlinear solution is described by the appropriate KdV equations. 12,13herefore, we will omit the majority of the derivation in these regions and instead focus on the coupled regularised Boussinesq equations in (12).Finally we will use the continuity conditions ( 7) - (10) to obtain 'initial conditions' for the derived leading order equations.
A. First Region: two homogeneous layers In the first region, the equation is identical in both homogeneous layers and therefore we assume the same incident wave in both, and consider asymptotic multiple-scale expansions of the type where the characteristic variables are given by ξ = x − t, η = x + t, and the slow variable X = x.Here, the functions I and R (1) , G (1) represent the leading order incident and reflected waves respectively and P (1) , Q (1) are the higher order corrections.Substituting the asymptotic expansion into the first equation in (11), the system is satisfied at leading order, while at O ( ) we have and a similar equation can be obtained for the second layer.We average equation (14) with respect to the fast space variable x using lim . . .dx, in the reference frame moving with the linear speed of rightand left-propagating waves (at constant ξ or η).Assuming that all functions and their derivatives remain bounded (in order to avoid secular terms in asymptotic expansions) and decay sufficiently rapidly at infinity we have, for example at constant ξ, The same result can be obtained for P (1) at constant η.Therefore, we can average (14) at constant ξ to obtain Similarly, averaging (14) with respect to x at constant η gives R X − 6R (1) R (1) η + R (1) ηηη η = 0.
In each case, we can integrate with respect to the relevant characteristic variable and, recalling that there are no waves at infinity, we obtain

R
(1) Substituting ( 18) and ( 19) into ( 14) and integrating with respect to the characteristic variables, we obtain an expression for P (1) of the form where φ (1) 1 are arbitrary functions.Similarly, we obtain the equations for the waves in the second layer, in addition to (18).The first radiation condition requires that the solution to the problem should not change the incident wave. 12,13or the case of an incident solitary wave this implies that φ 1 = 0 and φ

B. Second region: bi-layer with soft bonding
We assume that the layers have close properties, so that c − 1 = O ( ).In this case the cRB equations admit solutions in the form of radiating solitary waves. 26,29Thus, we assume that and therefore we can make the rearrangement This allows us to use one set of characteristic variables for f (2) and g (2) .Let us assume that there is a weakly nonlinear solution to (12) of the form The characteristic variables ξ, η and X are the same as before, T (2) and S (2) represent the transmitted waves in the second section of the bar, where T is for the top layer and S is for the bottom layer.Similarly, R (2) and G (2) are the reflected waves, and the higher order corrections in this section are given by P (2) and Q (2) , for the top and bottom layers respectively.The solution is considered over the period of time until the waves reflected from the boundary x = x b , between the second and the third region, reach the boundary x = x a , between the first and the second region.Moreover, the boundary x = x b is assumed to be sufficiently far away from the boundary x = x a , allowing us to use the averaging lim ing the asymptotic expansion into the equation for f (2) in ( 12) the equation is satisfied at leading order, while at O ( ) we have For the equation governing g (2) we have X + We average equations ( 24) and ( 25) with respect to the fast space variable x as defined earlier.In each case, we average at constant ξ or η and note that P (2) and Q (2) are both zero when averaged at either constant ξ or constant η.Averaging ( 24) and ( 25) at constant ξ gives T

S
(2) and therefore ( 26) and ( 27) form a system of coupled Ostrovsky equations. 29We note that the Ostrovsky equation was initially derived to describe long surface and internal waves in a rotating ocean. 40,41Coupled Ostrovsky equations appear naturally in the description of nonlinear waves in layered waveguides, both solid and fluid. 29,42imilarly, averaging ( 24) and ( 25) at constant η gives R (2)

G
(2) respectively.Therefore, to leading order, the transmitted and reflected waves are described by two systems of coupled Ostrovsky equations.This result is consistent with the time-averaged derivation. 29ubstituting ( 26) and ( 28) into ( 24) and integrating with respect to the characteristic variables, we obtain where φ (2) 1 are arbitrary functions.Similarly, substituting ( 27) and ( 29) into ( 25) and integrating with respect to the characteristic variables, we obtain where φ 2 are arbitrary functions.

C. Third region: delamination
We now consider the third region, where the same bilayered waveguide does not have a bonding layer, modelling delamination.The motion in this region is governed by two uncoupled regularised Boussinesq equations, but with differing coefficients in each layer.We look for a weakly nonlinear solution to (13) of the form where we now use two sets of characteristic variables ξ = x − t, η = x + t, and ν = x − ct, ζ = x + ct, while X = x.Substituting this into system (13) the equation is satisfied at leading order, while at O ( ) we have

−2Q
(3) We define the averaging in this region as Averaging at constant ξ and ν and integrating with respect to the appropriate characteristic variable, we obtain two KdV equations of the form

S
(3) Substituting the results for (34) into (32) and integrating with respect to the characteristic variables, we obtain where φ 1 are arbitrary functions.The second radiation condition states that if the incident wave is rightpropagating, then there should be no left-propagating waves in the transmitted wave field. 12,13Thus, ψ (3) 1 = 0. Similarly, substituting (35) into (33) and integrating with respect to the characteristic variables, we obtain where φ are arbitrary functions.By the same argument as above ψ (3) 2 = 0.

D. Matching at boundaries: continuity conditions
In order to find 'initial conditions' for the equations derived by the averaging, we collect the expressions for the weakly nonlinear solutions and substitute them into the continuity conditions ( 7) - (10).We first consider the continuity conditions for displacements for the time interval when the waves have not yet reached the third region.The displacement at negative infinity is assumed to be constant.Differentiating the continuity conditions (7) with respect to time at x = x a , and recalling that x , we obtain the following conditions in terms of the strain rates: Substituting the weakly nonlinear solutions obtained in Section III into (38) and noting that the reflected waves R (2) and G (2) in the second section are not yet present, we obtain at leading order We can integrate to obtain an expression at x = x a by noting that integration with respect to x can be reduced to integration with respect to the characteristic variable, as x appears linearly in the expressions for the characteristic variables.By the assumption that the strain waves are localised in the first two regions, when evaluated at either x = −∞ or x = x b the expression will be zero.Therefore, from (40) we obtain Similarly, from (39) we obtain Next, we consider the continuity conditions for displacements for the time interval when the localised strain waves are present in all three regions, but the waves reflected from the boundary x = x b between the second and the third region have not yet reached the boundary x = x a between the first and the second region.The displacement at positive infinity is assumed to be equal to zero (the waves propagate into the unperturbed media).Differentiating the continuity conditions (8) with respect to time at x = x b , and recalling that x , we obtain the following conditions in terms of the strain rates: Then, similarly to the previous considerations, we obtain We now make use of the continuity conditions for normal stress and, substituting the relevant weakly nonlinear solution into (9) (noting that u ) we obtain to leading order Similarly, substituting the relevant weakly nonlinear solutions into (10) we have, to leading order, We can now find 'initial conditions' for the systems describing transmitted and reflected waves in each section of the bar, expressed in terms of the given incident wave.
For the top layer we have T I| x=xa , R (2) where Similarly for the bottom layer we have T = 1.These coefficients agree with previous work for a perfectly bonded waveguide 12,13 and are intuitive, as we would expect a wave to be, to leading order, wholly transmitted when travelling along a layer of the same material.If the value of c varies between sections of the bar i.e. the material in a single layer varies across the bar, then the coefficients should be calculated using the respective values of c.
We note that the functions which remained undefined in the higher order corrections can be found by considering higher order terms in the equations of motion and the continuity conditions, similarly to the solution of the initial-value problems. 29,43However, this is beyond the scope of our present paper.

IV. Numerical modelling
In a case study we compare the results of the semianalytical numerical modelling, based on the leading order terms of the weakly nonlinear solution of the previous section, with the results of direct numerical simulations for the problem ( 4) - (10).We also compare numerical results with theoretical predictions for the amplitude of the lead soliton in the delaminated region, made using the IST.We solve the original Boussinesq equations using the finite-difference method described in Appendix A, and the weakly nonlinear solution derived in Section III using the pseudospectral method described in Appendix B. For the finite-difference method, in all cases, we use step sizes of ∆x = ∆t = 0.01 and, for the pseudospectral method, we use ∆ξ = 0.3 (the same step size is used for all characteristic variables) and ∆X = 0.001.In all cases, we assume α = 1.05, β = 1.05, c = 1 + /2 and = 0.05.We note that, similarly to the single Ostrovsky equation, the coupled Ostrovsky equations ( 26) -( 27) and ( 28) -( 29) imply zero mass constraints: Therefore we first use initial conditions for the incident strain solitary wave which include a pedestal term, 43 to guarantee zero mass, and then show that for this class of problems, one can also work with initial conditions in the form of pure strain solitary waves, without the pedestal terms.Indeed, in the latter case the zero mass constraints are still approximately satisfied, by the nature of the solutions of the problem, which we established in the direct numerical simulations using the finite-difference method.Thus, we use the following initial conditions for the displacements in (4) (corresponding localised initial conditions for the strains in (11) are given by the derivatives of these functions) where we have σ can be zero or one, 2L is the length of the domain used in the weakly nonlinear modelling, x 0 is an arbitrary phase shift, κ = ∆t, and the corresponding strain has zero mass (for σ = 1).The constant of integration is chosen so that the waves propagate into an unperturbed medium.The amplitude of the pedestal for the corresponding strains can be reduced by increasing the value of S. In all cases discussed here, S = 10 and we set x 0 = 0.For w(x, 0) and w(x, κ) we use the same expressions.If the initial condition was not given by an explicit analytic function, then we could deduce the second initial condition for the scheme by taking the forward difference approximation (simulations have shown that either case is sufficiently accurate).
For the weakly nonlinear solutions we take the exact solitary wave solution of the equation ( 18) governing the incident wave, with the same pedestal term (differentiated with respect to x) as used in (51), where Ã = − v1 2 , Λ = 2 √ v1 and where v 1 is related to v by the approximate relation v = 1+ v 1 +O 2 .For the initial conditions in other sections of the bi-layer, we make use of the relations in Section III D to obtain the initial conditions in terms of (52).

A. Solitons in the delaminated section
To obtain quantitative predictions for parameters of solitons in the delaminated section we use the IST applied to the KdV equations ( 34) and ( 35) derived in Section III C. Firstly we recall that the solution of an initial-value problem for the KdV equation on the inifinite line, for a sufficiently rapidly decaying initial condition U 0 (χ), is related to the spectral problem for the Schrödinger equation where λ is the spectral parameter. 2In particular, parameters of solitons are defined by the discrete spectrum of the equation (54).In our previous studies of the scattering of an incident solitary wave in the delaminated area of the perfectly bonded waveguide, the discrete spectrum could be found analytically. 12,13However, in the present study, we are dealing with the scattering of radiating solitary waves, generated in the two-layered bar with soft bonding, and scattered in the delaminated region of the bar.Therefore, we have to find the spectrum of the Schrödinger equation numerically.
To achieve this we implement a shooting method. 44We consider the potential U 0 (χ) for the Schrödinger equation, which is the initial condition in the KdV problem (53).It is well known that the discrete spectrum is bounded by the minimum of the initial condition (negative value) and zero. 45Since the potential U 0 (χ) is localised, the eigenfunctions have the asymptotic behaviour where λ = −r 2 .We rewrite the Schrödinger equation (54) in the form Ψ χ = Φ, Φ χ = [U 0 (χ) − λ] Ψ, and solve this system from the boundary χ = a to χ = b.We normalise the solution by setting Ψ(a) = 1 and Φ(a) = r.
The system is then solved using the standard Runge-Kutta 4 th order method.The ratio of the values of these two functions is tested at the right boundary against the relation Φ(b)/Ψ(b) = −r to determine if r is an eigenvalue.We start with the least possible value for an eigenvalue (the minimum of the initial condition) and iterate to zero in sufficiently small steps in order to determine the eigenvalues to the desired accuracy.In our present study we consider the cases when in each layer there is only one soliton in the delaminated region.We use the method described here to determine the parameters of this soliton (in each layer), and compare with the solitons emerging in our modelling.In other settings, multiple solitons can be generated by a single incident soliton. 12,13

B. Delamination of semi-infinite length
We first consider the bi-layer shown in Figure 3 and use the initial conditions (51) and (52) with zero mass, i.e. with σ = 1.The comparison between the two numerical approaches in this case can be seen in Figure 4.A radiating solitary wave is formed in the bonded section of the bar, shown at t = 600.When this radiating solitary wave enters the delaminated section of the bar, the soliton separates from the tail and becomes a classical soliton with dispersive radiation following behind.The agreement between the weakly nonlinear solution and the direct finite-difference technique is good for the solitons and reasonable for the tail, with a small phase shift introduced.The agreement is improved by reducing , and this has been tested for a number of values.
If we take the same initial conditions with non-zero mass, i.e. σ = 0, we obtain a similar result to the previously discussed case, as can be seen in Figure 5.The radiating solitary waves generated in the two layers are close to each other, and therefore the zero mass constraints for the difference of the two solutions are approximately satisfied (see Appendix B).We now apply the IST framework to the waves entering the delaminated section of the bar, as the behaviour of the transmitted waves in the two layers in this section is governed by two separate KdV equations.We numerically solve the scattering problem for the related Schrödinger equation, as discussed in Section IV A, to obtain the eigenvalues.Since there is only one discrete eigenvalue for each layer of the waveguide, the long time asymptotic behaviour of the solution of the appropriate KdV equation consists of one soliton and dispersive radiation, which in the canonical form ( 53) is given by where r is defined by the eigenvalue λ = −r 2 , and χ 0 is the phase shift.We use the theoretical predictions to justify the numerical schemes in Appendices A and B. In each layer, the height of the soliton found using these schemes has been compared with the theoretical prediction using the IST, to confirm that the numerical schemes resolve the behaviour of the system correctly.The theoretical (IST) predictions and the numerical results for the height of the soliton are compared in Table I.
In the case with zero mass initial condition, the prediction of the heights using the IST underestimates the numerical solution, as the solitons have not yet fully sepa-  Comparison of amplitudes for solitons in the delaminated area for both layers with the predicted value using the IST, for zero mass (σ = 1) and non-zero mass (σ = 0) initial conditions.rated from the negative pedestal.In the case with initial condition having non-zero mass, the agreement between the theoretical predictions and the numerical results is excellent.

C. Delamination of finite length
Let us now pose a question.Is it possible to determine if there is a delamination in some part of the bar, between two bonded regions?A graphical representation of this structure is shown in Figure 6, and all considerations of Section III are extended to this situation.We know that transmitted waves will propagate in the delaminated area For the finite-difference method, the full computational domain is [−1000, 1000].In the pseudospectral method, N = 16384.
x < x a 0 < x < x b x a < x < 0 x > x b O z y x FIG.6: Bi-layer with two homogeneous layers for x < x a , a bonded two-layered section for x a < x < 0, a delaminated section for 0 < x < x b and another bonded two-layered section for x > x b .
with speeds close to the characteristic speeds of the linear waves, and therefore the time it will take for the wave to travel through a delaminated region, of length l, can be estimated as T i ≈ l/c i , where i represents the layers in the bar.Indeed, when the radiating solitary waves enter the delaminated region, as seen in Figure 7, the solitons propagate with speeds close to their respective characteristic speeds.When these solitons enter the second bonded region they again generate radiating solitary waves.If the separation between the two solitons is sufficiently large when they enter the second bonded region, we see a distinctive double-humped wave of significantly reduced amplitude -a clear sign of delamination.
However if the delamination area is shorter, then the solitons in the delaminated section will not be fully separated.In this case, the radiating waves in the second bonded region overlap and generate a new single-humped radiating solitary wave.Irrespective of the separation, this process of creating a new radiating solitary wave is accompanied by some additional radiation, and therefore the amplitude of the new radiating solitary wave is reduced in both layers when compared to the radiating solitary wave propagating in a fully bonded waveguide, with no delamination.Furthermore, as the radiating solitary wave is not supported by the KdV equation, in the delaminated region the radiation separates from the soliton and the periodicity observed in the tail disappears as the tail transforms into a wave packet.This feature gives another indication that delamination is present.
In order to investigate this behaviour more fully, we consider several cases with different delamination lengths, as measured in wavelengths of the solitary wave.The wavelength is measured using the common measure Full Width at Half Magnitude (FWHM).In this case, the FWHM of the incident soliton measures approximately 5 units.We present results for delamination of length 10, 20, 40 and 60 FWHM, and the case where there is no delamination.Note that Figure 7 is for a delamination length of approximately 60 FWHM.
025, σ = 0 and = 0.05: direct numerical simulations (solid line) and weakly nonlinear solution (dashed line).Two homogeneous layers, of the same material as the upper layer, are on the left, and the waves propagate to the right.For the finite-difference method, the full computational domain is [−600, 1000].In the pseudospectral method, N = 8192.
We note that the results presented in Figure 9 are obtained using the semi-analytical method.The finitedifference method solves for two sections of the bar at a time and therefore we must wait until the wave and its tail are fully contained in the region before moving the calculation domain.However, for a shorter delamination i.e. 20 FWHM or less, the leading wave front will reach the boundary of the calculation domain before the tail has fully entered this region.Therefore, the wave will either reflect and interfere with our solution, or we will lose part of the tail when we move the calculation domain.This is a natural limitation for the use of the finite-difference method in its present form.This could be remedied by solving for all sections of the bar simultaneously, however this will be much more expensive.We see from Figure 9 that there are some key differences between the model without delamination and the model with delamination.We only show the waves in the 'top' layer as the waves in the 'bottom' layer are similar.Firstly, as the length of delamination increases, the amplitude of the radiating solitary wave created in the second bonded region is reduced.This can be explained by the fact that the waves in the delaminated section of the bar travel at different speeds in each layer and will be incident on the second bonded region at different times, so the energy exchange between layers results in the generation of a radiating solitary wave of reduced amplitude.A graph of the amplitudes against the delamination length, in FWHM, is presented in Figure 10.We can clearly see that after an initial growth period (up to 20 FWHM), the dependence is close to linear.The presence of the doublehumped solution, as seen in the image for 60 FWHM, clearly identifies a delamination.Further numerical experiments have shown that this double-humped structure is emerging at around 45-50 FWHM.Furthermore, the small hump behind the lead soliton in the 40 FWHM image is the start of a double-humped solution, but the second 'hump' has a similar amplitude to the radiation and therefore is not distinct.The speed of the waves in the delaminated region is different to the bonded region, and therefore when the new radiating solitary wave is formed in the second bonded region, it will have a phase shift.Measuring from the minima of the waves, we calculate a phase shift of 0.2, 0.8, 2.7 and 3.6 for 10, 20, 40 and 60 FWHM respectively, growing with the delamination length as expected.
The radiating solitary wave is not a solution of the KdV equation and therefore, in the delaminated region, the radiating tail forms a wave packet, breaking the regularity of the tail region.This feature is again more pronounced for a larger delamination area, however it can already be clearly identified for a short delamination, such as in the case of 10 FWHM as seen in Figure 9.Further experiments have identified this behaviour for 5 FWHM, however the amplitude is similar to the rest of the tail and therefore this is difficult to identify visually.We summarise these observations as follows.For very short delamination areas i.e. 5 FWHM, these differences are negligible and suggests that the soliton does not take care of delaminations shorter than a threshold value.For delamination areas that are greater than 5 FWHM, we can use our observations above for the amplitude reduction and phase shift to identify the presence of delamination.Modifying the FWHM value of the incident soliton can help identify shorter delaminations.

D. Further experiments
From a physical standpoint, we would like to test a given structure in as many ways as possible to obtain all possible information.Let us assume that we have a structure such as that in Figure 6 but with the two homogeneous layers removed.Given this structure, there are four natural tests that we can conduct: with two homogeneous layers of either the same material as the top or bottom layers, and attaching the homogeneous layers to either the left-hand side or right-hand side of the structure, with the waves propagating to the right or to the left, respectively.Examples using the same material as the top layer are shown in Figure 7 and Figure 8 for the homogeneous layers being attached on the left-hand side and right-hand side respectively.We discuss all results here but omit the other cases for brevity.Firstly we note that the double-humped structure is present in the second bonded region in all cases, confirming that each test identifies the presence of a sufficiently long delamination area, from our observations from Section IV C.There is a small phase shift between the results for the different materials in the homogeneous layers, arising from the higher characteristic speed of the material of the bottom layer.We also note that, for the homogeneous layers being present on the right, the first bonded region is longer and this leads to a longer radiating tail.This tail becomes a wave packet in the delaminated region and we observe that the larger amplitude wave packet is closer to the double-humped structure for the case where the tail is longer, i.e. when the homogeneous section is on the right-hand side.Indeed, the length of the bonded FIG.9: A comparison between the case without delamination (solid lines) and with delamination (dashed lines) of differing lengths, measured in FWHM of the incident soliton.The model is the same as that used in Figure 7 with the same parameters, and all images are for t = 1200.region after the delamination is shorter in this case and therefore we would expect the wave packet to be closer to the leading wave.This gives us an indication of where the delamination is present in the bi-layer, i.e. if the radiation wave packet is closer to the leading wave when sending the waves from the right, then the delamination is closer to the left-hand side of the structure, and vice versa.
Another useful feature is that the generated wave is of a larger amplitude in the case when the homogeneous layers are of the same material as the bottom layer (with a larger characteristic speed), and therefore the FWHM measure is smaller.In addition, the amplitude difference to the case with no delamination is even clearer.

V. Conclusions
In this paper we studied the scattering of a long radiating bulk strain solitary wave in a delaminated bi-layer with a soft bonding between the layers.The modelling was performed within the framework of the system of coupled regularised Boussinesq equations (1), which were derived to describe long nonlinear longitudinal waves in a twolayered waveguide with a soft bonding layer ('imperfect interface') from a layered lattice model with all essential degrees of freedom of a layered elastic waveguide 26 .For a single layer, the model leads to a 'doubly dispersive equation' (DDE) 20,21 , earlier derived for the long longitudinal waves in a bar of rectangular cross-section using the nonlinear elasticity approach 12 .In dimensional variables, the DDE for a bar of rectangular cross-section σ = 2a × 2b has the form where The direct numerical modelling of this type of problem is difficult and expensive because one needs to solve several boundary value problems linked to each other via matching conditions at the boundaries.Therefore, we developed an alternative semi-analytical approach based upon the use of several matched asymptotic multiplescale expansions and averaging with respect to the fast space variable.The developed approach is an extension of our earlier work, 13 where we considered a simple bilayer with perfect bonding.Unlike our earlier work, the bi-layer with the soft ('imperfect') bonding does not support the usual solitary waves.They are replaced by radiating solitary waves, with a one-sided oscillatory tail, as discussed in the Introduction.We modelled the generation and subsequent scattering of these radiating solitary waves in a number of complex waveguides with and without a delamination area, as well as predicting the parameters of the lead solitons generated in the delaminated area using the IST framework for the relevant KdV equations.
The developed direct finite-difference scheme and the scheme for the weakly nonlinear solution show good agreement in all regions of the bi-layer, with a small difference in the amplitude and minor phase shift between the results.This could be remedied by the inclusion of higher-order corrections in the weakly nonlinear scheme, similarly to initial-value problems. 29,43We also note that the direct finite-difference scheme is expensive in comparison to the weakly nonlinear scheme.
Our study has revealed key features of the behaviour of radiating solitary waves in such delaminated bi-layers, for different lengths of the delaminated area compared to the wavelength (FWHM) of the incident soliton.If the delaminated area is sufficiently long (≥ 25 FWHM), then there is a significant reduction in the amplitude of the transmitted radiating solitary wave (≥ 10 %).In fact, the incident radiating solitary wave undergoes a complicated process of shedding a tail and propagating with slightly different speeds along the two layers in the delaminated region, followed by generation of a new radiating wave in the second bonded region.For shorter delamination regions (< 25 FWHM), the key dynamical effect manifesting the presence of a delaminated region in the structure is the appearance of a wavepacket in the regular tail of the radiating solitary wave.The waves are not sensitive to very short delamination regions, comparable to the wavelength of the incident soliton.In practice, using an admissible incident soliton with smaller wavelength (and higher amplitude), would increase the sensitivity to shorter delamination regions.If the delaminated region is longer (≥ 45 FWHM), the separation of solitons, propagating in two layers in the delaminated region, leads to the emergence of a double-humped radiating wave in the second bonded region.We did not show the cases with delamination areas greater than 60 FWHM.The dynamical behaviour in these cases is simpler, leading to the emergence of two distinct radiating solitary waves in each layer of the second bonded region -a very clear sign of delamination.However, such cases are likely to be uncommon in real-world applications because of the dissipation processes which have not been accounted for in our modelling.Typical values of elastic moduli for the PMMA (polymethylmethacrylate) and PS (polystyrene) and experimental data for solitons in layered PMMA/PS bars of 10 × 10 mm cross-section have been reported previously. 24,46The typical amplitude of the strain for the observed compression solitary waves is O 10 −4 , and the soliton velocity is about 5−7% greater than the linear longitudinal wave speed. 23It has been reported that, in PMMA and PS, solitons can propagate to distances tens of times greater than their width without significant decay. 22,23he generation of a radiating bulk strain solitary wave and subsequent disappearance of the 'ripples' in the delaminated area of a two-layered PMMA bar with the PCP (polychloroprene-rubber-based) adhesive has been observed in experiments. 24Our numerical modelling motivates further laboratory experimentation with a wide range of materials used in practical applications.It also paves the way for similar studies in other physical settings supporting radiating solitary waves and radiating dispersive shock waves, for example, in nonlinear optics. 27,28he solution interval is discretised by N equidistant points with spacing ∆ξ = 2π/N , where N is a power of 2. This allows us to use the Discrete Fourier Transform We now introduce the Runge-Kutta 4th-order method.Assume that the solution at X is given by ûj = u(k, jκ) and ŵj = u(k, jκ), where κ = ∆X.Then the solution at X = (j + 1)∆X is given by ûk,(j+1)κ = ûk,jκ + 1 6 (a 1 + 2b for i = 1, 2. The functions F i are found as a rearrangement of (B2) to contain all non-time derivatives.Explicitly we have For the cases in Section IV C, where the waves re-enter a coupled region after a delamination, we need to introduce a linear damping region ('sponge layer') and add this at each end of the domain to prevent radiated waves reentering the region of interest and interfering with the main wave structure. 48The sponge layer is defined as for some constants ν, K.We choose K so that KL = 12 and ν is chosen so that damping occurs quickly.The sponge layer is incorporated as (u X + uu x + u xxx + r(x)u) x = δ (u − w) , (w X + w x + ww x + w xxx + r(x)w) x = γ (w − u) , (B5) and is treated in the same way as nonlinear terms.To remove aliasing effects, we use the truncation 2/3-rule by Orszag in Boyd. 50This effect is due to the pollution of the numerically calculated Fourier transform by higher frequencies due to the series being truncated.Finally, for the particular non-zero mass initial conditions used in this paper, in the Fourier space the zero mode (for k = 0) is approximated as a small constant O 10 −2 , defined by the initial condition.The maximum of the zero mode of u − w is O 10 −3 for the cases considered in the paper, and therefore this approximation introduces only a small error, approximately satisfying the equations (B2) (multiplied by k) for zero modes.In the real space, the average 1 2L L −L (u − w)dξ is O 10 −5 , showing that the zero mass constraint is approximately satisfied.A more accurate approach is required for general initial conditions with non-zero mass. 43

FIG. 7 :
FIG. 7:The waves f (top row) and g (bottom row) in the various sections of the bi-layer, for α = β = 1.05, c = 1.025, δ = γ = 1, v = 1.025, σ = 0 and = 0.05: direct numerical simulations (solid line) and weakly nonlinear solution (dashed line).Two homogeneous layers, of the same material as the upper layer, are on the left, and the waves propagate to the right.For the finite-difference method, the full computational domain is [−600, 1000].In the pseudospectral method, N = 8192.

FIG. 8 :
FIG.8:The waves f (top row) and g (bottom row) in the various sections of the bi-layer, for α = β = 1.05, c = 1.025, δ = γ = 1, v = 1.025, σ = 0 and = 0.05: direct numerical simulations (solid line) and weakly nonlinear solution (dashed line).Two homogeneous layers, of the same material as the upper layer, are on the right, and the waves propagate to the left.For the finite-difference method, the full computational domain is [−400, 1200].In the pseudospectral method, N = 8192.

FIG. 10 :
FIG.10:The percentage decrease in amplitude for a given delamination measured in FWHM.

TABLE I :
ρ is the density, E is the Young's modulus, ν is the Poisson's ratio, while l, m, n are the Murnaghan's moduli.Nondimensionalisation, regularisation of the dispersive terms and scaling bring the equation to the form (2).
1 + 2c 1 + d 1 ) , + 2b 2 + 2c 2 + d 2 ) , (B3)where a i , b i , c i , d i are functions of ξ at a given moment in time, X, and are defined as