Scattering of an Ostrovsky wave packet in a delaminated waveguide

We examine the scattering of an Ostrovsky wave packet, generated from a solitary wave, in a two layered waveguide with a delamination in the centre and soft (imperfect) bonding either side of the centre. The lower layer of the waveguide is assumed to be significantly denser than the upper layer, leading to a system of Boussinesq–Klein– Gordon (BKG) equations. Direct numerical modelling is difficult and so a semi-analytical approach consisting of several matched asymptotic multiple-scale expansions is used, which leads to Ostrovsky equations in soft bonded regions and Korteweg–de Vries equations in the delaminated region. The semi-analytical approach and direct numerical simulations are in good agreement with each other and theoretical estimates. The dispersion relations are used to estimate the wave speed and hence classify the length of the delamination, in addition to changes in the amplitude of the wave packet. We also show how to scale the non-dimensional results to material variables and an example for PMMA is presented. These results can provide a tool to control the integrity of layered structures.


Introduction
The propagation of long longitudinal bulk strain solitary waves in elastic waveguides has been observed in a variety of elastic waveguides, such as rods, plates, shells and bars [1][2][3][4][5][6], and more recently the generation of undular bores has been considered [7].The dynamical behaviour of layered structures depends both on the properties of the bulk material in the layers, and on the type of the bonding between the layers.For rectangular bars with a perfect bond between the layers, the displacements can be modelled using Boussinesq-type equations [8,9] and have been observed experimentally [10][11][12][13].When the bonding between the layers is sufficiently soft and the materials of the layers have similar properties, bulk strain solitons are replaced by radiating solitary waves [12,14,15], a solitary wave with a co-propagating oscillatory tail, and theoretical estimates for the tail have been obtained [16].
When the materials of the layers are distinctly different, it is suggested that bulk strain solitons are replaced by Ostrovsky wave packets, following the theory for the corresponding initial-value problem [17].The Ostrovsky equation originally arose in the context of shallow water waves where the rotation of the Earth is taken into account [18] and the evolution of wave packets generated from an initial pulse has been extensively studied [19][20][21].So far the scattering of such wave packets in a region of delamination is yet to be explored, and that will be the focus of this paper.
The Ostrovsky equation requires that any regular localised solution has zero mass, however it can arise in the context of fluids where such a restriction is not enforced.This contradiction is resolved in [22] in the infinite interval and more Fig. 1.Examples of the bi-layer structure with an initial delaminated region for x 0 < x < x 1 , a soft bonded region for x 1 < x < x 2 and either (a) a delaminated region for x 2 < x < x 3 , or (b) a delaminated region for x 2 < x < x 3 followed by another soft bonded region for x 3 < x < x 4 .In both cases we assume that the material in the lower layer is significantly denser than the upper layer.
recently, for the finite interval, the contradiction is resolved for the Boussinesq-Klein-Gordon (BKG) equation, whose leading order weakly-nonlinear solution is the Ostrovsky equation [17,23,24].A similar contradiction is also resolved for a system of coupled Boussinesq equations when the characteristic speeds are close [25].These considerations are for equations on the finite interval and not in the context of scattering problems, where Ostrovsky equations can arise in the case of a layered waveguide when the materials of the layers are sufficiently different.To study the scattering of a wave packet in such a structure while circumventing the zero-mass contradiction, we consider the solution for a class of localised initial conditions that have zero mean.
In this paper we consider the scattering problem for a two layered bar with a delamination in the centre of soft bonding, where the material between the layers is distinctly different.In particular, we assume that the lower layer is significantly more dense than the upper layer, so that the displacements can be assumed to be zero in the lower layer and the strains in the system are described by BKG equations, rather than coupled Boussinesq equations.The paper is organised as follows.In Section 2 we set the problem formulation for the scattering of zero-mass wave packets in a bi-layer where the lower layer is significantly denser than the upper layer.In Section 3 we construct the weakly-nonlinear solution to the problem using an asymptotic multiple-scales expansion, which can be used in a semi-analytical method for the solution of the problem.In Section 4 the direct numerical solution to the problem is compared to the semi-analytical solution constructed here and the results are shown to be in good agreement with each other, as well as theoretical predictions.Two cases are considered: a soft bonded region followed by infinite delamination, and the case of delamination 'sandwiched' between soft bonded regions, which we use to identify methods for classifying the delamination.We also show how to scale the non-dimensional results to a physical setting, with an example for a PMMA bar.We conclude in Section 5.

Problem formulation
Let us consider the scattering of a wave packet in a two-layered structure with a soft bond between the layers.Two examples of this structure are shown in Fig. 1, as inspired by the experimental setup in [12], one with a semi-infinite delamination and one with a finite delamination.The scattering of long longitudinal waves in this structure was previously studied in [15] for the case when the characteristic speeds in the layers differed by O (ε), where ε is a small amplitude parameter.In this case, a radiating solitary wave is generated in the bonded sections of the structure and it was shown that the solitary wave separates from its co-propagating tail in the delaminated sections.Here we consider the situation when the difference in the characteristic speeds is more significant.Specifically we model the limiting case where the lower layer has a significantly higher characteristic speed in contrast to the upper layer, so we can assume that the longitudinal displacements in the lower layer are zero.
We assume that a long longitudinal solitary wave with a zero mean pedestal propagates within the split area of the waveguide between x 0 and x 1 , before entering the soft bonded region between x 1 and x 2 , then a delaminated region between x 2 and x 3 , and finally another bonded region between x 3 and x 4 .This is illustrated in Fig. 1(b) and follows the structure used in [15].
Each layer of the waveguide is assumed, without loss of generality, to be the same material throughout, and the lower layer of the waveguide is considered to be of much denser material, such that the displacements are essentially zero.The displacements in the waveguide are modelled by the system of scaled, regularised, non-dimensional Boussinesq or Boussinesq-Klein-Gordon (BKG) equations u (1) tt − u (1)  xx = 2ε ( −6u (1)   x u (1)   xx + u (1)   ttxx ) , tt − u (2)  xx = 2ε ( −6u (2)   x u (2)   xx + u (2)  ttxx − γ u (2)   ) , tt − u (3)  xx = 2ε ( −6u (3)   x u (3)   xx tt − u (4)  xx = 2ε ( −6u (4)   x u (4)   xx + u (4)  ttxx − γ u (4)   ) , with the displacement in each section of the waveguide denoted by u (i) , where i indicates the section of the waveguide, and we have the small amplitude parameter ε.This system of equations is complemented with continuity conditions on the interface between the sections, namely we have continuity of longitudinal displacement and continuity of normal stress We have presented the scaled form of the equations with the coefficients chosen in such a way that the equations derived in the weakly-nonlinear solution are in canonical form.Following the approach used in [15], we write our system in terms of the longitudinal strains rather than the displacements, and so differentiating our system (1)-(4) and setting u (i) x = f (i) , the strains are described by the system f (1) tt − f (1)  xx = 2ε (2) ) , ) , where we rewrite the continuity conditions in terms of strains as and In [15] we showed that seeking a weakly-nonlinear solution to this problem, when the characteristic speeds are close, led to a system of coupled Ostrovsky equations in the bonded sections and uncoupled Korteweg-de Vries (KdV) equations in the delaminated sections of the bar.This is consistent with the time-averaged derivation for the initial-value problem, as studied in [17].It was noted that the coupled Ostrovsky equations require that the difference between the solutions in the layers have zero mean, in a similar way to the single Ostrovsky equation.Recent works have shown that, for the initial-value problem, this can be overcome by explicitly calculating the mean value of the BKG or coupled Boussinesq equations and subtracting this from the function to obtain a system of equations with zero mean.The weakly-nonlinear solution is then constructed in the regime where all functions necessarily have zero mean [17,23].This approach has yet to be adapted for the scattering problem considered here.To circumvent this issue we consider localised initial conditions that already satisfy the zero mean requirement.Such an initial condition was considered in the first example of [15].

Weakly-nonlinear solution
We seek a weakly-nonlinear solution to the equation system ( 7)- (10), in addition to a space-averaging method, as was done for layered waveguides with either perfect bonding [8,9] and soft bonding [15].This is in contrast to the timeaveraging method that was used for the initial-value problem for the cRB equations [17,25,26].We will summarise the derivation below, following the established method outlined in [15].
In each region we assume a weakly-nonlinear solution consisting of a leading order transmitted and reflected wave, with higher-order corrections, while in the last region we have no reflected wave, and in the first region we have an incident wave, not a transmitted wave.We assume that all functions in the weakly-nonlinear expansion and their derivatives are bounded and sufficiently rapidly decaying at infinity.This is substituted into the governing equation in each region and space-averaging is applied to obtain equations for the leading order functions.Finally, the continuity conditions are used to obtain 'initial conditions' for these derived equations.

Region 1: Homogeneous section
The structure presented in Fig. 1 has a delaminated region attached to the bar on the left-hand side, which allows for a solitary wave to be generated in this region and is consistent with experiments [12].Therefore we consider an asymptotic multiple-scales expansion of the type where the characteristic variables are given by ξ = x−t, η = x+t, and the slow space variable X = εx.Here, the functions I and R (1) , represent the leading order incident and reflected waves respectively and P (1) is the higher-order corrections.
Substituting into (7) and averaging with respect to the fast space variable x gives

Region 2: Soft bonded section
Let us now consider the case where the layers have a soft bonding.The solution to the BKG equation ( 8) should take the form of Ostrovsky wave packets [17,18].We assume that there is a weakly-nonlinear solution to (8) of the form where the characteristic variables ξ , η and X are the same as for the first section.The function T (2) represents the transmitted wave and R (2) is the reflected wave, and the higher-order corrections in this section are given by P (2) .
The solution is considered over the period of time until the waves reflected from the boundary x = x 2 , between the second and the third region, reach the boundary x = x 1 , between the first and the second region.Applying the space-averaging we find 2) T (2)  ξ + T (2)   ξ ξ ξ where ( 16) and ( 17) are Ostrovsky equations [18].The Ostrovsky equation was initially derived to describe long surface and internal waves in a rotating ocean [18,27], and the long-time behaviour has been studied on the infinite domain [19].
The appearance of Ostrovsky equations agrees with the time-averaged derivation for the initial-value problem [17].We obtain the higher-order corrections as where φ (2) , ψ (2) are arbitrary functions.

Region 3: Delamination
In the third region of the bar, where we model a delamination (no bonding between the layers), we look for a weakly-nonlinear solution to (9) of the form where the characteristic variables have the same meaning as before.Substituting the weakly-nonlinear expansion into (9) and averaging, we obtain two KdV equations of the form

R
(3) and the higher-order corrections are given by where again φ (3) and ψ (3) are arbitrary functions.
Higher-order corrections are given by where φ (4) , ψ (4) are arbitrary functions.Radiation conditions state that if the incident wave is right-propagating and there is no boundary for reflected waves, then there should be no left-propagating waves in the transmitted wave field [8,9,15].
We have now defined the weakly-nonlinear solution in each section of the structure shown in Fig. 1.If we were to extend the existing structure, for example adding another delaminated region, the weakly-nonlinear solution in the new region can be formed using the results in Sections 3.1-3.4.

Matching at the boundaries
The continuity conditions ( 11) and ( 12) can be used to determine 'initial conditions' for the equations derived in Sections 3.1-3.4by substituting the weakly-nonlinear solution into the continuity conditions and expressing the value of the functions in the subsequent region, at the interface between the two regions, in terms of the value of the appropriate function in the former region.This follows the approach given in [15] so we summarise here.
For the first region, we consider the conditions for displacements in the time interval when the waves have not yet reached the third region and we assume that the displacement at negative infinity is zero, which is consistent with taking zero mean initial conditions.Differentiating (11) with respect to time at x = x 1 , we obtain a relation between the strain rates of the form As we assume that the waves have not reached the interface between the second and third sections of the bar, the reflected wave R (2) can be omitted from the calculation at this boundary when considering right-propagating waves.Using the weakly-nonlinear solutions derived in Sections 3.1 and 3.2, we substitute the expansions into ( 24) and obtain at leading order where we use the fact our waves are localised.From the second continuity condition we obtain, to leading order, Therefore, the 'initial conditions' for the reflected wave in the first section of the bar and the transmitted wave in the second section of the bar, expressed in terms of the given incident wave, are Here we have the reflection coefficient C (1) R = 0 and transmission coefficient C (1)   T = 1, which is consistent with previous works for a perfectly bonded waveguide [8,9] and for a layered bar with an imperfect bond [15].Indeed these coefficients describe full transmission of the incident wave at leading order, where the material in the bar is the same in each section for a given layer, as expected.Following the same approach for the interfaces between the second and third sections of the bar, and the third and fourth sections of the bar, assuming that the waves are travelling in these sections and have yet to reflect back to the left boundary, gives the 'initial conditions' as T T (3)   | x=x 3 . ( As the material is the same in all the sections of the bar, we find that C (2,3) R = 0 and C (2,3) T = 1, as we found for the first and second sections.This method could, in theory, be extended to higher-order functions but this is beyond the scope of our paper.

Numerical results
We now compare the direct numerical simulations of the equation system ( 1)-( 4), complemented with the continuity conditions ( 5) and ( 6), to the derived weakly-nonlinear solution in Sections 3.1-3.4.We will use a semi-analytical method, as the KdV and Ostrovsky equations were derived analytically but must be solved numerically.To compute the direct numerical simulations of the system (1)-( 4) we make use of the method presented in [28], with ∆x = 0.01 and ∆t = 0.01.
For the semi-analytical method we will solve the KdV and Ostrovsky equations using a pseudospectral method, based upon the method for the KdV equation in [15].For both the KdV and Ostrovsky equations, we take ∆ξ = 0.1 (or ∆η = 0.1) and N = 131, 072.
The initial condition is assumed to be the solitary wave solution of the KdV equation in the first region of the waveguide, which is not bonded.As the solution in the bonded section of the waveguide must have zero mean because of the Ostrovsky equation, we follow the method used in [15,23] to add a pedestal term that gives zero mean.Namely we take where A = −v/2, Λ = 2/ √ v, and we set Here x 0 is the initial phase shift, x 1 is chosen so that the pedestal terms are close to the initial condition and S is a scale factor on the width of the pedestal.In all cases we choose S = 10 and x 1 = 5.For the initial condition of ( 1)-( 4) we take the integral of the strain solitary wave initial condition for the Boussinesq equation in ( 7), so we have where we have [15] ) .
Note that here u(x, 0) refers to the initial condition for the entire waveguide, where the initial condition for u (i) is taken as u(x, 0) in the appropriate domain.We will consider two cases: a two layered bar with infinite delamination on the right and a two layered bar with delamination in the centre, which are graphically represented by Fig. 1.

Amplitude of solitons in delaminated region
It is well known that localised initial conditions for the KdV equation will either produce a series of solitons and some dispersive radiation (unless the continuous spectrum is zero) or they will only generate dispersive radiation, as shown by the theory of the Inverse Scattering Transform (IST) [29].If the initial condition is a soliton, then the parameters of the generated solitons can be found analytically [8,30].In our case we expect the wave entering the delaminated region of the bar will be either an Ostrovsky wave packet or the earlier evolution of this wave packet [18].Therefore, to obtain quantitative predictions for the parameters of the generated solitons, we make use of the IST applied to the KdV equation (19) derived in Section 3.3, in a similar way to what was done in [15] for a radiating solitary wave.We summarise the method below.
Let us recall that the solution of an initial-value problem for the KdV equation, on the infinite line, for a sufficiently rapidly decaying initial condition U 0 (χ), is related to the spectral problem for the Schrödinger equation given as where λ is the spectral parameter [29].The parameters of the solitons are defined by the discrete spectrum of (35) and therefore we aim to describe this spectrum numerically, using a shooting method [31].As we have negative nonlinear coefficient, the discrete spectrum will be bounded by the minimum of the initial condition and zero [32], thus giving a range over which to find the discrete eigenvalues.The eigenfunctions of the Schrödinger equation have the asymptotic behaviour where λ = −r 2 .Therefore we rewrite (35) as a boundary value problem in the form where we solve this system from the boundary χ = a to χ = b.The solution is normalised by setting Ψ (a) = 1, Φ(a) = r, and the system is solved with a Runge-Kutta 4th order method.The ratio of Φ and Ψ is tested at the right boundary against the relation Φ(b)/Ψ (b) = −r to determine if r is an eigenvalue.

Infinite delamination
We first consider the bi-layer waveguide with infinite delamination shown in Fig. 1(a) and compare the direct numerical simulations to the semi-analytical results.The comparison between the methods is shown in Fig. 2 and we can see that there is good agreement in all regions.The incident soliton is evolving into an Ostrovsky wave packet in the bonded region, agreeing with the expected results from the initial-value problem in [17].This wave packet fissions into solitons in the delaminated region and there is a leading soliton which is clearly separated from the radiation behind, while the other peaks are yet to separate, escape the well and move on a zero background, so cannot be clearly identified as solitons yet.The agreement could be improved by including higher-order terms or reducing the value of ε.
The theoretical predictions using the IST, as outlined in Section 4.1, suggest there are three solitons, with the lead soliton's amplitude predicted as −1.016, in direct agreement with the semi-analytical results on which the prediction is based.The second and third solitons are predicted as having amplitudes −0.392 and −0.042 respectively.We can potentially see a second soliton in Fig. 2, which has a larger amplitude than the prediction, however it has yet to escape the well and will lose some amplitude as it does so.The calculation would have to be run for a very large time to attain the other predicted values.

Linear dispersion relation
To further understand the behaviour of the wave packet, we consider the linear dispersion relation for the KdV and Ostrovsky equations.It was shown in [19] that the Ostrovsky wave packet moves at a speed close to the maximum of the group speed, while the solitary wave solutions of the KdV equation exist in the spectral gap.We can use these results to help us determine the position of the wave packet in the waveguide, and thus determine if a delamination occurs.We derive the dispersion relations here for the moving reference frame.
The linear dispersion relations are obtained by seeking solutions of the linearised form of ( 13) and ( 16) (or equivalently (19) and ( 22)), by removing the nonlinear terms.We seek solutions in the form I ≈ I 0 e i(kξ −ωX) and T 2 ≈ T 0 e i(kξ −ωX) , where ω is the wave frequency and k is the wavenumber.Substituting into the linearised form of ( 13) and ( 16) gives the dispersion relation The same relationship can be found for ( 19) and ( 22) respectively.The phase speed is found as p = ω/k, so Table 1 The predicted speed of the Ostrovsky wave packet in the bonded regions, v b , and the numerically observed speed of the waves in each region, denoted v 1 to v 3 .The speed of solitons in the delaminated region is found in the spectral gap, namely for p > 0 and close to p = 0.The group speed is found as The peak of this group speed for T 2 is c g = −2.449at k 0 = 0.639.If the bar is not delaminated the wave packet will move with this speed, while if the bar is delaminated the solitary wave solutions in the delaminated region will move with speed p > 0. This means that, for each unit in X within the delaminated region, we have a minimum phase shift of x s = 2.449 (assuming soliton of amplitude close to 0).
Translating back to (x, t) we find the speed of the waves in the bonded region is predicted as Table 1 denotes the wave speeds in each region of the waveguide modelled in Section 4.2 and the predicted speed from the dispersion relation.We can see that the solitary waves move at a speed greater than 1, while in the bonded region we have a good agreement with the predicted speed from the dispersion relation.

Finite delamination
The next logical case to consider is when the delamination is no longer semi-infinite and we 'sandwich' the delamination region with two soft bonded regions, as illustrated in Fig. 1(b).Our aim is to determine if a delamination region is present and, if so, can we infer its length and position using the observed wave.
We modify the infinite delamination case in Fig. 2 so that we have an additional soft bonded region located at 1600 < x < 2400.The output from this calculation for the final soft bonded region is shown in Fig. 3.We can see that the leading solitary wave is evolving into an Ostrovsky wave packet, as was observed in the first bonded region, and there is radiation behind, which can potentially mask any further wave packets generated by other solitons.Thus, we will focus on the leading wave packet.We also note that the agreement is good around the leading wave packet and weaker away from the wave packet, due to the lack of higher-order correction terms.
To understand how the length of the delaminated region affects a bonded structure, we must investigate the effect of varying the delamination length on the waves in the final bonded region.Firstly we consider the structure outlined in Section 4.4 with the same parameters.The delamination length is varied from D = 0 (assume the bar has no delamination) in intervals of 10 up to D = 360.The maximum time of the calculation is taken to be t = 2000.
We consider two cases, ε = 0.005 and ε = 0.001.The corresponding waves incident on the delaminated region for these two cases are shown in Fig. 4. As the Ostrovsky equation has 'time' appearing as X = εx, we expect that a lower value of ε will mean less evolution, so ε = 0.001 is a wave that is closer to a soliton, while for ε = 0.005 it is a wave packet, giving us two types of behaviour to study.
Firstly we consider the change in amplitude as the delamination length varies.The amplitude of the wave packet varies as it travels, within a wave packet which has a peak amplitude at its centre [19].Therefore, to determine the peak amplitude, we find the maximum amplitude of the wave packet in the second bonded region, between t = 1600 and  For the case when ε = 0.001 we can see an almost linear relationship between delamination length and amplitude change.This is consistent with results for the solitary wave and radiating solitary wave [15,28].For ε = 0.005 we initially find that the amplitude is increasing, then as the delamination length increases the amplitude decreases, however in a nonlinear fashion.This is expected as the wave packet contains one or more peaks that will evolve into solitons in the delaminated region, however they are contained within a wave packet and so their amplitude initially can be artificially amplified by the radiation.For larger delaminations, the generated solitons start to separate from the wave packet and become distinct, giving the expected amplitude drop as they evolve into Ostrovsky wave packets in the second bonded region.
For the phase shift, we measure the position of the front of the wave packet.We see that, after an initial period of adjustment from D = 0 to D = 10, the case of ε = 0.001 evolves mostly linearly, with a small phase shift overall.In contrast, for ε = 0.005 there is a clear increase in phase shift as delamination length increases.The relationship is nonlinear, although after an initial settling region it behaves linearly for a period, beginning to taper for large delamination length.Referring to Section 4.3, for ε = 0.005 we expect a minimum phase shift of x s = 2.449 for D = 1/ε = 200, which is close to the gradient in the centre of the line, although the behaviour is nonlinear so we expect some deviation in the results.For the case when ε = 0.001 we expect a minimum phase shift of x s = 2.449 for D = 1000.The phase shift line in Fig. 5(b) for ε = 0.001 is close to this gradient as well.

Material parameters
Experiments have been performed for PMMA (polymethylmethacrylate) and PS (polystyrene), where long longitudinal waves have been observed [10][11][12][13].Therefore, we naturally consider whether we can take material parameters from a physical setting and model the generation of these wave packets.
Firstly, we note that the long longitudinal displacements in a homogeneous waveguide can be described by a 'doubly dispersive equation' (DDE) [1,2].This DDE can be derived in a bar of rectangular cross-section using nonlinear elasticity As the length of the delamination region increases, the amplitude of the Ostrovsky wave packet decreases, in a mostly linear fashion for small values of small wave parameter ε and nonlinearly for larger values of ε, with an initial growth period before decaying.Similarly, we found an increasing phase shift with relation to the fully bonded case, which behaves linearly for small values of ε and nonlinearly for larger values of ε.The linear dispersion relation was used to predict the speed of the waves in each section and this prediction is in good agreement with the numerically observed results, although the behaviour is more complicated.Experiments for PMMA and polystyrene layered waveguides have been reported previously [12,13] and showed that solitons can be detected at larger distances than linear waves [34].Thus we simulated a PMMA bar of cross-section 10 mm × 10 mm and of length 600 mm.We found an Ostrovsky wave developed and, in a delaminated region, solitons begin to form.Our numerical modelling motivates laboratory experimentation with a wide range of materials to determine if this measure can be used to control delamination.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3.The solution at t = 2400 in the final section of Fig. 1(b), with region length 1600 < x < 2400.The parameters are the same as Fig. 2, comparing the direct numerical simulation (blue, solid line) and semi-analytical solution (red, dashed line).

Fig. 5 .
Fig. 5. Plots of the change in (a) amplitude, measured as a percentage change in amplitude to the fully bonded case, and (b) phase shift, in reference to the fully bonded case, for ε = 0.001 (blue, solid line) and ε = 0.005 (red, dashed line).

Fig. 6 .
Fig.6.Example of Ostrovsky wave packet generation and scattering in a bi-layer structure, using material parameters for PMMA.We assume that ε = 0.1 and we have converted all measurements back to dimensional variables.