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

TL;DR

This paper develops a semi-analytical method to study the scattering of radiating solitary waves in layered elastic structures with delamination, providing insights into wave dynamics and potential structural integrity control.

Contribution

A semi-analytical approach using matched asymptotic expansions to analyze radiating solitary waves in bi-layered structures with delamination, validated against numerical simulations.

Findings

Semi-analytical method agrees well with numerical simulations.

Radiating solitary waves influence the integrity of layered structures.

Wave dynamics can be controlled through understanding of scattering processes.

Abstract

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.

Figures

Tables

Table 1. Table 1 : Comparison of amplitudes for solitons in the delaminated area for both layers with the predicted value using the IST, for zero mass ( σ = 1 𝜎 1 \sigma=1 ) and non-zero mass ( σ = 0 𝜎 0 \sigma=0 ) initial conditions.

Regime	Layer	Numerical	Theoretical
$\sigma =1$	1	-0.2545	-0.2473
$\sigma =1$	2	-0.2301	-0.2192
$\sigma =0$	1	-0.2979	-0.2979
$\sigma =0$	2	-0.2680	-0.2680

Equations

$$f_{tt}-f_{xx}$$

$$g_{tt}-c^2{g}_{xx}$$

$$f_{tt}-f_{xx}=\frac{1}{2}(f^2{)}_{xx}+f_{ttxx}.$$

$$[k^2(1-p^2)-k^4{p}^2+\delta ][k^2(c^2-p^2)-\beta k^4{p}^2+\gamma ]=\gamma \delta ,$$

$$u_{tt}^{(1)}-u_{xx}^{(1)}$$

$$w_{tt}^{(1)}-w_{xx}^{(1)}$$

$$u_{tt}^{(2)}-u_{xx}^{(2)}$$

$$-\delta \left(u^{(2)}-w^{(2)}\right)],$$

$$w_{tt}^{(2)}-c^2{w}_{xx}^{(2)}$$

$$+\gamma \left(u^{(2)}-w^{(2)}\right)]$$

$$u_{tt}^{(3)}-u_{xx}^{(3)}$$

$$w_{tt}^{(3)}-c^2{w}_{xx}^{(3)}$$

$$u^{(1)}{|}_{x=x_a}$$

$$u^{(2)}{|}_{x=x_b}$$

$${u_x^{(1)}+ϵ\left[-6{\left(u_x^{(1)}\right)}^2+2u_{ttx}^{(1)}\right]|}_{x=x_a}=$$

$${u_x^{(2)}+ϵ\left[-6{\left(u_x^{(2)}\right)}^2+2u_{ttx}^{(2)}\right]|}_{x=x_a},$$

$${w_x^{(1)}+ϵ\left[-6{\left(w_x^{(1)}\right)}^2+2w_{ttx}^{(1)}\right]|}_{x=x_a}=$$

$${c^2{w}_x^{(2)}+ϵ\left[-6\alpha {\left(w_x^{(2)}\right)}^2+2\beta w_{ttx}^{(2)}\right]|}_{x=x_a};$$

$${u_x^{(2)}+ϵ\left[-6{\left(u_x^{(2)}\right)}^2+2u_{ttx}^{(2)}\right]|}_{x=x_b}=$$

$${u_x^{(3)}+ϵ\left[-6{\left(u_x^{(3)}\right)}^2+2u_{ttx}^{(3)}\right]|}_{x=x_b},$$

$${c^2{w}_x^{(2)}+ϵ\left[-6\alpha {\left(w_x^{(2)}\right)}^2+2\beta w_{ttx}^{(2)}\right]|}_{x=x_b}=$$

$${c^2{w}_x^{(3)}+ϵ\left[-6\alpha {\left(w_x^{(3)}\right)}^2+2\beta w_{ttx}^{(3)}\right]|}_{x=x_b};$$

$$f_{tt}^{(1)}-f_{xx}^{(1)}$$

$$g_{tt}^{(1)}-g_{xx}^{(1)}$$

$$f_{tt}^{(2)}-f_{xx}^{(2)}$$

$$-ϵ\delta \left(f^{(2)}-g^{(2)}\right),$$

$$g_{tt}^{(2)}-c^2{g}_{xx}^{(2)}$$

$$+ϵ\gamma \left(f^{(2)}-g^{(2)}\right)$$

$$f_{tt}^{(3)}-f_{xx}^{(3)}$$

$$g_{tt}^{(3)}-c^2{g}_{xx}^{(3)}$$

$$f^{(1)}=I\left(\xi ,X\right)+R^{(1)}\left(\eta ,X\right)+ϵP^{(1)}\left(\xi ,\eta ,X\right)+\mathcal{O}\left({ϵ}^2\right),$$

$$g^{(1)}=I\left(\xi ,X\right)+G^{(1)}\left(\eta ,X\right)+ϵQ^{(1)}\left(\xi ,\eta ,X\right)+\mathcal{O}\left({ϵ}^2\right),$$

$$-2P_{\xi \eta }^{(1)}={\left(I_X-6I{I}_{\xi }+I_{\xi \xi \xi }\right)}_{\xi }+(R_X^{(1)}-6R^{(1)}{R}_{\eta }^{(1)}$$

$$+{R_{\eta \eta \eta }^{(1)})}_{\eta }-6\left(2I_{\xi }{R}_{\eta }^{(1)}+R_{\eta \eta }^{(1)}I+I_{\xi \xi }{R}^{(1)}\right),$$

$$\underset{\chi \to -\infty }{\lim }\frac{1}{x_a-\chi }\int {\text{displaylimits}}_{\chi }^{x_a}{P}_{\xi \eta }^{(1)}\mathrm{d}x=\underset{\chi \to -\infty }{\lim }\frac{1}{2(x_a-\chi )}\int {\text{displaylimits}}_{2\chi -\xi }^{2x_a-\xi }{P}_{\xi \eta }^{(1)}\mathrm{d}\eta$$

$$=\underset{\chi \to -\infty }{\lim }\frac{1}{2(x_a-\chi )}{\left[P_{\xi }^{(1)}\right]}_{2\chi -\xi }^{2x_a-\xi }=0.$$

$${\left(I_X-6I{I}_{\xi }+I_{\xi \xi \xi }\right)}_{\xi }=0.$$

$${\left(R_X^{(1)}-6R^{(1)}{R}_{\eta }^{(1)}+R_{\eta \eta \eta }^{(1)}\right)}_{\eta }=0.$$

$$I_X-6I{I}_{\xi }+I_{\xi \xi \xi }$$

$$R_X^{(1)}-6R^{(1)}{R}_{\eta }^{(1)}+R_{\eta \eta \eta }^{(1)}$$

Full text

radiating solitary wave; multiple-scale expansions; coupled Ostrovsky equations

**Long longitudinal bulk strain solitary waves observed in elastic waveguides, such as rods, bars, plates and shells, can be modelled with Boussinesq-type equations. Radiating solitary waves, that is solitary waves radiating a co-propagating one-sided oscillatory tail, emerge in layered elastic waveguides with soft (‘imperfect’) bonding between the layers. In this paper we study the scattering of a radiating solitary wave in delaminated areas of imperfectly bonded bi-layers within the scope of the coupled Boussinesq equations. We develop direct and semi-analytical numerical approaches and demonstrate that radiating solitary waves undergo changes which could be used to control the quality of the interfaces. **

I Introduction

The discovery of solitons as extremely stable localised coherent structures ZK is intrinsically linked with the discovery of the Inverse Scattering Transform (IST) for the Korteweg-de Vries (KdV) equation - the method for the solution of a large class of initial-value problems on the infinite line. GGKM The latter has shown that solitons, when present, constitute the main part of the long-time asymptotics of initial-value problems for localised initial data, and this is the reason why solitons proved to be a very important part of the physical world we live in, across all scales. SS ; NM ; Hoefer ; Dreiden15

Initially developed as a purely analytical technique, AS ; ZS in recent years the IST formed the basis for the development of efficient numerical approaches to the analysis of nonlinear problems, most notably within the framework of another famous integrable model, the Nonlinear Schrödinger (NLS) equation. ZS ; Osborne ; SCDA An efficient IST-based numerical approach to solving the KdV equation was also developed. TOD

Recently, 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. KS ; KT This condensed matter problem can be viewed as an analogue to the fluid mechanics problem of calculating the reflected and transmitted waves when a surface or internal soliton passes through an area of rapid depth variation. Pelinovsky ; Tappert ; Johnson ; Djordjevic ; Grimshaw08 The theoretical predictions KS agreed well with experimental studies, Dreiden10 and were also confirmed by numerical simulations. KT Long longitudinal bulk strain solitary waves were experimentally observed in various elastic waveguides, including rods, bars, plates and shells, and modelled using Boussinesq-type equations. Samsonov_book ; Porubov_book ; Dreiden10 ; Dreiden15 The exceptional stability of bulk strain solitons Samsonov_book ; SDS_2008 ; SSB_2016 makes them an attractive candidate for the introscopy of layered structures, in addition to the existing methods. Dreiden12 ; SBDPS

The 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. KSZ The radiating strain solitary wave has recently been observed in laboratory experiments. Dreiden12 More 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 SG ; CTMK 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 KSZ (in non-dimensional and scaled form):

[TABLE]

Here, $f$ and $g$ denote the longitudinal strains in the layers, while the coefficients $c,\alpha,\beta,\delta,\gamma$ are defined by the physical and geometrical parameters of the problem KSZ (see Section II for details).

In the symmetric case ($c=\alpha=\beta=1$) system (1) admits the reduction $g=f$, where $f$ satisfies the equation

[TABLE]

The Boussinesq equation (2) has particular solitary wave solutions: ${\displaystyle f=3(v^{2}-1)\ \mathrm{sech}^{2}\left(\frac{x-vt}{\Lambda}\right),\quad\Lambda=\frac{2v}{\sqrt{v^{2}-1}},}$ 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. $c$ is close to 1), these pure solitary wave solutions are replaced with radiating solitary waves, KSZ ; KM that is solitary waves radiating a co-propagating one-sided oscillatory tail. BGK ; Shrira ; Bona Figure 1 illustrates the pure solitary wave solutions of system (1) with $\delta=\gamma=0$ (for a fixed value of $v$), and subsequent evolution of this initial condition into a radiating solitary wave.

Radiating solitary waves have been extensively studied in the context of perturbed KdV equations, coupled KdV systems, perturbed NLS equations, and coupled NLS systems. Boyd ; Grimshaw ; GI ; GJ ; Karpman ; Lombardi ; VB 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. KSZ The dispersion relation has the form

[TABLE]

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\to 0$, while approaching zero as $k\to\infty$.

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, Ostrovsky ; review 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. KM We 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.

II Problem Formulation

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. Dreiden12 Two 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 bi-layer are assumed to have close properties, leading to the generation of a radiating solitary wave in the bonded section. KSZ We study the scattering of this wave by the subsequent delaminated region.

The mathematical problem formulation consists of the following sets of scaled regularised non-dimensional equations in the respective sections of the complex waveguide: KS ; KSZ ; KT

[TABLE]

for $x<x_{a}$;

[TABLE]

for $x_{a}<x<x_{b}$; and

[TABLE]

for $x>x_{b}$. Here, the functions $u^{(i)}(x,t)$ and $w^{(i)}(x,t)$ describe longitudinal displacements in the ‘top’ and ‘bottom’ layers of the three sections of the waveguide, respectively. The values of the constants $\alpha$, $\beta$ and $c$ depend on the physical and geometrical properties of the waveguide, while the constants $\delta$ and $\gamma$ depend on the properties of the soft bonding layer, and $\epsilon$ is a small amplitude parameter. KSZ

These equations are complemented with continuity conditions at the interfaces between the sections: continuity of longitudinal displacement

[TABLE]

and continuity of normal stress

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

as well as some natural initial and boundary conditions, which will be imposed later.

III Weakly Nonlinear Solution

Differentiating the governing equations with respect to $x$ and denoting $f^{(i)}=u_{x}^{(i)}$ and $g^{(i)}=w_{x}^{(i)}$, we obtain the equations ‘in strains’:

[TABLE]

for $x<x_{a}$;

[TABLE]

for $x_{a}<x<x_{b}$; and

[TABLE]

for $x>x_{b}$.

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 time-averaging method used for the homogeneous initial-value problem. KM All 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. KS ; KT Therefore, 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.

III.1 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

[TABLE]

where the characteristic variables are given by $\xi=x-t$, $\eta=x+t$, and the slow variable $X=\epsilon 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 $\mathcal{O}\left(\epsilon\right)$ we have

[TABLE]

and a similar equation can be obtained for the second layer. We average equation (14) with respect to the fast space variable $x$ using $\displaystyle\lim_{\chi\rightarrow-\infty}\frac{1}{x_{a}-\chi}\int_{\chi}^{x_{a}}\dots\;\;\mathrm{d}x,$ in the reference frame moving with the linear speed of right- and left-propagating waves (at constant $\xi$ or $\eta$). 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 $\xi$,

[TABLE]

The same result can be obtained for $P^{(1)}$ at constant $\eta$. Therefore, we can average (14) at constant $\xi$ to obtain

[TABLE]

Similarly, averaging (14) with respect to $x$ at constant $\eta$ gives

[TABLE]

In each case, we can integrate with respect to the relevant characteristic variable and, recalling that there are no waves at infinity, we obtain

[TABLE]

Substituting (18) and (19) into (14) and integrating with respect to the characteristic variables, we obtain an expression for $P^{(1)}$ of the form

[TABLE]

where $\phi_{1}^{(1)}$, $\psi_{1}^{(1)}$ are arbitrary functions.

Similarly, we obtain the equations

[TABLE]

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. KS ; KT For the case of an incident solitary wave this implies that $\phi_{1}^{(1)}=0$ and $\phi_{2}^{(1)}=0$.

III.2 Second region: bi-layer with soft bonding

We assume that the layers have close properties, so that $c-1=\mathcal{O}\left(\epsilon\right)$. In this case the cRB equations admit solutions in the form of radiating solitary waves. KSZ ; KM Thus, we assume that

[TABLE]

and therefore we can make the rearrangement

[TABLE]

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

[TABLE]

The characteristic variables $\xi$, $\eta$ 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 $\displaystyle\lim_{x_{b}\to\infty}\frac{1}{x_{b}-x_{a}}\int_{x_{a}}^{x_{b}}\dots\ dx.$ Substituting the asymptotic expansion into the equation for $f^{(2)}$ in (12) the equation is satisfied at leading order, while at $\mathcal{O}\left(\epsilon\right)$ we have

[TABLE]

For the equation governing $g^{(2)}$ we have

[TABLE]

We average equations (24) and (25) with respect to the fast space variable $x$ as defined earlier. In each case, we average at constant $\xi$ or $\eta$ and note that $P^{(2)}$ and $Q^{(2)}$ are both zero when averaged at either constant $\xi$ or constant $\eta$. Averaging (24) and (25) at constant $\xi$ gives

[TABLE]

and therefore (26) and (27) form a system of coupled Ostrovsky equations. KM We note that the Ostrovsky equation was initially derived to describe long surface and internal waves in a rotating ocean. Ostrovsky ; review Coupled Ostrovsky equations appear naturally in the description of nonlinear waves in layered waveguides, both solid and fluid. KM ; AGK14

Similarly, averaging (24) and (25) at constant $\eta$ gives

[TABLE]

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. KM

Substituting (26) and (28) into (24) and integrating with respect to the characteristic variables, we obtain

[TABLE]

where $\phi_{1}^{(2)}$, $\psi_{1}^{(2)}$ are arbitrary functions.

Similarly, substituting (27) and (29) into (25) and integrating with respect to the characteristic variables, we obtain

[TABLE]

where $\phi_{2}^{(2)}$, $\psi_{2}^{(2)}$ are arbitrary functions.

III.3 Third region: delamination

We now consider the third region, where the same bi-layered 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

[TABLE]

where we now use two sets of characteristic variables $\xi=x-t$, $\eta=x+t$, and $\nu=x-ct$, $\zeta=x+ct$, while $X=\epsilon x$. Substituting this into system (13) the equation is satisfied at leading order, while at $\mathcal{O}\left(\epsilon\right)$ we have

[TABLE]

We define the averaging in this region as $\displaystyle\lim_{\chi\to\infty}\frac{1}{\chi-x_{b}}\int_{x_{b}}^{\chi}\dots\ dx$. Averaging at constant $\xi$ and $\nu$ and integrating with respect to the appropriate characteristic variable, we obtain two KdV equations of the form

[TABLE]

Substituting the results for (34) into (32) and integrating with respect to the characteristic variables, we obtain

[TABLE]

where $\phi_{1}^{(3)}$, $\psi_{1}^{(3)}$ are arbitrary functions. The second radiation condition states that if the incident wave is right-propagating, then there should be no left-propagating waves in the transmitted wave field. KS ; KT Thus, $\psi_{1}^{(3)}=0$.

Similarly, substituting (35) into (33) and integrating with respect to the characteristic variables, we obtain

[TABLE]

where $\phi_{2}^{(3)}$, $\psi_{2}^{(3)}$ are arbitrary functions. By the same argument as above $\psi_{2}^{(3)}=0$.

III.4 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 $f^{(i)}=u^{(i)}_{x}$, $g^{(i)}=w^{(i)}_{x}$, we obtain the following conditions in terms of the strain rates:

[TABLE]

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

[TABLE]

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=-\infty$ or $x=x_{b}$ the expression will be zero. Therefore, from (40) we obtain

[TABLE]

Similarly, from (39) we obtain

[TABLE]

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 $f^{(i)}=u^{(i)}_{x}$, $g^{(i)}=w^{(i)}_{x}$, we obtain the following conditions in terms of the strain rates:

[TABLE]

Then, similarly to the previous considerations, we obtain

[TABLE]

We now make use of the continuity conditions for normal stress and, substituting the relevant weakly nonlinear solution into (9) (noting that $u_{x}^{(i)}=f^{(i)}$ and $w_{x}^{(i)}=g^{(i)}$) we obtain to leading order

[TABLE]

Similarly, substituting the relevant weakly nonlinear solutions into (10) we have, to leading order,

[TABLE]

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

[TABLE]

where $C_{R}^{(i)}=0$ and $C_{T}^{(i)}=1$ for all $i$.

Similarly for the bottom layer we have

[TABLE]

where ${\displaystyle D_{R}^{(1)}=\frac{c-1}{c+1},D_{T}^{(1)}=\frac{2}{c\left(1+c\right)},D_{R}^{(2)}=0,D_{T}^{(2)}=1.}$ These coefficients agree with previous work for a perfectly bonded waveguide KS ; KT 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. KM ; KMP14 However, this is beyond the scope of our present paper.

IV Numerical modelling

In a case study we compare the results of the semi-analytical 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 $\Delta x=\Delta t=0.01$ and, for the pseudospectral method, we use $\Delta\xi=0.3$ (the same step size is used for all characteristic variables) and $\Delta X=0.001$. In all cases, we assume $\alpha=1.05$, $\beta=1.05$, $c=1+\epsilon/2$ and $\epsilon=0.05$.

We note that, similarly to the single Ostrovsky equation, the coupled Ostrovsky equations (26) - (27) and (28) - (29) imply zero mass constraints:

[TABLE]

Therefore we first use initial conditions for the incident strain solitary wave which include a pedestal term, KMP14 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)

[TABLE]

where we have $A=-\frac{v\sqrt{v^{2}-1}}{\sqrt{2\epsilon}}$, $\Lambda=\frac{2\sqrt{2\epsilon}v}{\sqrt{v^{2}-1}}$ and

[TABLE]

$\sigma$ 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, $\kappa=\Delta t$, and the corresponding strain has zero mass (for $\sigma=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,\kappa)$ 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),

[TABLE]

where $\tilde{A}=-\frac{v_{1}}{2}$, $\tilde{\Lambda}=\frac{2}{\sqrt{v_{1}}}$ and

[TABLE]

where $v_{1}$ is related to $v$ by the approximate relation $v=1+\epsilon v_{1}+\mathcal{O}\left(\epsilon^{2}\right)$. For the initial conditions in other sections of the bi-layer, we make use of the relations in Section III.4 to obtain the initial conditions in terms of (52).

IV.1 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.3. Firstly we recall that the solution of an initial-value problem for the KdV equation

[TABLE]

on the inifinite line, for a sufficiently rapidly decaying initial condition $U_{0}(\chi)$, is related to the spectral problem for the Schrödinger equation

[TABLE]

where $\lambda$ is the spectral parameter. GGKM In 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. KS ; KT However, 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. SB We consider the potential $U_{0}(\chi)$ 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. LL Since the potential $U_{0}(\chi)$ is localised, the eigenfunctions have the asymptotic behaviour

[TABLE]

where $\lambda=-r^{2}$. We rewrite the Schrödinger equation (54) in the form $\Psi_{\chi}=\Phi$, $\Phi_{\chi}=\left[U_{0}\left(\chi\right)-\lambda\right]\Psi$, and solve this system from the boundary $\chi=a$ to $\chi=b$. We normalise the solution by setting $\Psi(a)=1$ and $\Phi(a)=r$. The system is then solved using the standard Runge-Kutta 4th order method. The ratio of the values of these two functions is tested at the right boundary against the relation $\Phi(b)/\Psi(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. KS ; KT

IV.2 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 $\sigma=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 $\epsilon$, and this has been tested for a number of values.

If we take the same initial conditions with non-zero mass, i.e. $\sigma=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.1, 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

[TABLE]

where $r$ is defined by the eigenvalue $\lambda=-r^{2}$, and $\chi_{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 1.

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 separated 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.

IV.3 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 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}\approx 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.

We note that the results presented in Figure 9 are obtained using the semi-analytical method. The finite-difference 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 double-humped 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.

IV.4 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.3. 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 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 two-layered waveguide with a soft bonding layer (‘imperfect interface’) from a layered lattice model with all essential degrees of freedom of a layered elastic waveguide KSZ . For a single layer, the model leads to a ‘doubly dispersive equation’ (DDE) Samsonov_book ; Porubov_book , earlier derived for the long longitudinal waves in a bar of rectangular cross-section using the nonlinear elasticity approach KS . In dimensional variables, the DDE for a bar of rectangular cross-section $\sigma=2a\times 2b$ has the form

[TABLE]

where $c=\sqrt{E/\rho},c_{1}=c/\sqrt{2(1+\nu)},\beta=3E+2l(1-2\nu)^{3}+4m(1+\nu)^{2}(1-2\nu)+6n\nu^{2},J=4ab(a^{2}+b^{2})/3$, $\rho$ is the density, $E$ is the Young’s modulus, $\nu$ is the Poisson’s ratio, while $l,m,n$ are the Murnaghan’s moduli. Non-dimensionalisation, regularisation of the dispersive terms and scaling bring the equation to the form (2).

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 multiple-scale expansions and averaging with respect to the fast space variable. The developed approach is an extension of our earlier work, KT where we considered a simple bi-layer 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. KM ; KMP14 We 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 ($\geq 25$ FWHM), then there is a significant reduction in the amplitude of the transmitted radiating solitary wave ($\geq 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 ($\geq 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\times 10$ mm cross-section have been reported previously. Dreiden12 ; Dreiden14 The typical amplitude of the strain for the observed compression solitary waves is $\mathcal{O}\left(10^{-4}\right)$, and the soliton velocity is about $5-7\%$ greater than the linear longitudinal wave speed. SSB_2016 It has been reported that, in PMMA and PS, solitons can propagate to distances tens of times greater than their width without significant decay. SDS_2008 ; SSB_2016

The 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. Dreiden12 Our 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. SG ; CTMK

Acknowledgments

We thank Alexander Strohmaier and Greg Roddick for useful discussions of numerical approaches to the calculation of the spectrum of the Schrödinger equation. The authors would like to acknowledge the support of the Engineering and Physical Sciences Research Council (EPSRC). M.R.T. is supported by an EPSRC studentship.

Appendix A Numerical formulation for Boussinesq equations

The problem (4) - (10) is treated as a set of boundary value problems (BVPs) with the continuity conditions defining the interaction between the equations. In order to determine the solution at a given time in a section of the bar, we solve by a method similar to that presented in our earlier work.KT The formulation presented there was for two sections of a bi-layer, rather than three. While the formulation can be extended to a bi-layer with three sections, it will significantly increase the complexity of the numerical scheme and hence the computation time.

To alleviate this, we consider a modified approach where we compute the solution in two sections of the bi-layer at a time. Referring to Figure 3 as an example, we would first calculate for the region $x<x_{a}$ and $x_{a}<x<x_{b}$ with two constraints: that the generated radiating solitary wave has not reached the boundary $x=x_{b}$, and that the waves reflected from $x=x_{a}$ have not reflected from the left-hand boundary of the problem and into the domain $x_{a}<x<x_{b}$. The speed of the incident solitary wave and the radiating solitary wave will be close to the characteristic speed in each section of the bar, and therefore both constraints can be satisfied by choosing an appropriate time interval for the calculation. A similar approach is then followed for the region $x_{a}<x<x_{b}$ and $x>x_{b}$, where we choose an appropriate time so that the waves in this section have not reflected from the boundary $x=x_{b}$ and back again from the boundary $x=x_{a}$.

We summarise the method for, say, the second and third sections. Discretising the domain $\left[x_{a},L\right]\times\left[0,T\right]$ into a grid with equal spacings $h=\Delta x$ and $\kappa=\Delta t$, the analytical solutions $u^{(n)}\left(x,t\right)$, $w^{(n)}\left(x,t\right)$ are approximated by the exact solution of the finite difference schemes $u^{(n)}\left(ih,j\kappa\right)$, $w^{(n)}\left(ih,j\kappa\right)$, denoted $u_{i,j}^{(n)}$ and $w_{i,j}^{(n)}$ respectively. We make use of first-order and second-order central difference approximations in the main equations for $u_{x}^{(n)}$, $u_{xx}^{(n)}$, $u_{tt}^{(n)}$ and the approximation for $u_{ttxx}^{(n)}$ can be obtained iteratively using the approximations for $u_{tt}^{(n)}$ and $u_{xx}^{(n)}$. Similar approximations are used for $w$. To simplify the obtained expressions, we introduce the notation

[TABLE]

Using the finite-difference approximations in system (4) and making use of the notation above, we obtain a tri-diagonal system of the form

[TABLE]

and

[TABLE]

for $x_{a}<x<x_{b}$ and a similar system for the third section of the bi-layer. Assuming the domain can be discretised, we denote the central point as $N=\frac{x_{b}-x_{a}}{h}$ and therefore conditions (7) and (8) translates directly to

[TABLE]

In the continuity conditions (9) and (10), we make use of the central difference approximations presented above, and introduce ‘ghost points’ of the form $u_{N+1,j+1}^{(2)}$, $u_{N-1,j+1}^{(3)}$, $w_{N+1,j+1}^{(2)}$ and $w_{N-1,j+1}^{(3)}$. The continuity condition relates the ghost points to each other (omitted for brevity, but the expressions are of the same form as the ones in our earlier workKT ).

As we are considering localised initial data for strains, if we take $L$ large enough then we can enforce zero boundary conditions for the strains, i.e. $u_{x}=0$, $w_{x}=0$. We note that (57) - (58) and the associated system for $x>x_{b}$ are similar to the system considered in our previous work and the same method can be applied here. This is summarised as follows: solve each tridiagonal system (57) - (58) (and the systems for $x>x_{b}$) in terms of the ghost points at the central boundary, then use the expressions for the ghost points in conditions (59) and the discretisation of (9) and (10) to obtain a solvable algebraic system for the ghost points. This is then substituted back into the tridiagonal systems to obtain the solution.

Appendix B Numerical formulation for coupled Ostrovsky equations

In previous studies, the uncoupled Ostrovsky equationOS12 and coupled Ostrovsky equationsAGK13 have been solved using finite-difference and pseudospectralF96 techniques respectively. We implement a pseudospectral method using the Runge-Kutta 4th-order method for ‘time’ stepping in the Fourier space, and the nonlinear terms are calculated in the real domain and transformed back to the Fourier space for use in the calculation.

Firstly we transform the solution interval from $\left[-L,L\right]$ to $\left[0,2\pi\right]$ via the transform $\tilde{\xi}=s\xi+\pi$ and $s=\pi/L$. We denote the dependent variables by $u(\xi,X)$ and $w(\xi,X)$ and, assuming unity coefficients for simplicity, we obtain (omitting tildes)

[TABLE]

The nonlinear term is calculated in the real domain, then transformed to the Fourier space. We can rewrite the nonlinear terms by introducing the notation

[TABLE]

The solution interval is discretised by $N$ equidistant points with spacing $\Delta\xi=2\pi/N$, where $N$ is a power of 2. This allows us to use the Discrete Fourier Transform

[TABLE]

and a similar transform for $w$. Here $k$ is an integer which represents the discretised (and scaled) wavenumber. In what follows we first exclude the zero mode, i.e. we let $k\neq 0$, assuming zero mass initial conditions for $u$ and $w$.

The inverse transform is

[TABLE]

and again a similar transform for $w$. We make use of the Fast Fourier Transform (FFT) algorithm to implement these transforms efficiently. The discrete Fourier transform of equations (60) with respect to $\xi$ gives

[TABLE]

We now introduce the Runge-Kutta 4$\mathrm{th}$-order method. Assume that the solution at $X$ is given by $\hat{u}_{j}=u(k,j\varkappa)$ and $\hat{w}_{j}=u(k,j\varkappa)$, where $\varkappa=\Delta X$. Then the solution at $X=(j+1)\Delta X$ is given by

[TABLE]

where $a_{i},b_{i},c_{i},d_{i}$ are functions of $\xi$ at a given moment in time, $X$, and are defined as

[TABLE]

for $i=1,2$. The functions $F_{i}$ are found as a rearrangement of (61) to contain all non-time derivatives. Explicitly we have

[TABLE]

For the cases in Section IV.3, 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 re-entering the region of interest and interfering with the main wave structure.AGK13 The sponge layer is defined as

[TABLE]

for some constants $\nu,K$. We choose $K$ so that $KL=12$ and $\nu$ is chosen so that damping occurs quickly. The sponge layer is incorporated as

[TABLE]

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.Boyd01 This 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 $\mathcal{O}\left(10^{-2}\right)$, defined by the initial condition. The maximum of the zero mode of $u-w$ is $\mathcal{O}\left(10^{-3}\right)$ for the cases considered in the paper, and therefore this approximation introduces only a small error, approximately satisfying the equations (61) (multiplied by $k$) for zero modes. In the real space, the average $\frac{1}{2L}\int_{-L}^{L}(u-w)d\xi$ is $\mathcal{O}\left(10^{-5}\right)$, showing that the zero mass constraint is approximately satisfied. A more accurate approach is required for general initial conditions with non-zero mass.KMP14