## Abstract

Integrity of layered structures, extensively used in modern industry, strongly depends on the quality of their interfaces; poor adhesion or delamination can lead to a failure of the structure. Can nonlinear waves help us to control the quality of layered structures? In this paper, we numerically model the dynamics of a long longitudinal strain solitary wave in a split, symmetric layered bar. The recently developed analytical approach, based on matching two asymptotic multiple-scales expansions and the integrability theory of the Korteweg–de Vries equation by the inverse scattering transform, is used to develop an effective semi-analytical numerical approach for these types of problems. We also employ a direct finite-difference method and compare the numerical results with each other, and with the analytical predictions. The numerical modelling confirms that delamination causes fission of an incident solitary wave and, thus, can be used to detect the defect.

## 1. Introduction

The birth of the concept of a *soliton* is ultimately linked with the Boussinesq equation appearing in the continuum approximation of the famous Fermi–Pasta–Ulam problem [1]. Although the vast majority of the following studies of solitons were devoted to waves in fluids and nonlinear optics (see [2–4] and references therein) it was also shown, within the framework of the nonlinear dynamic elasticity theory, that Boussinesq-type equations can be used to model the propagation of long nonlinear longitudinal bulk strain waves in rods and plates [5,6]. The existence of longitudinal bulk strain solitons in these solid waveguides, predicted by the model equations, was confirmed by experiments [7–9]. Very recently the theoretical and experimental studies were extended to some types of adhesively bonded layered bars [10–14] and thin-walled cylindrical shells [15].

Unusually slow decay of the observed solitons in some polymeric waveguides makes them an attractive candidate for the introscopy of the waveguides, in addition to the existing methods. Recently, we analytically studied scattering of the longitudinal bulk strain soliton by the inhomogeneity modelling delamination in a symmetric layered bar with perfect bonding [10]. The developed theory was supported by experiments [16,17]. In order to extend the research to other, more complicated types of layered structures, the analytical and experimental studies need to be complemented with numerical modelling. The direct numerical modelling of even the simplest problem proved to be rather expensive. The development of numerical schemes based on analytical results is an actively developing direction of research (e.g. asymptotic techniques have been used to improve the efficiency of numerics for tsunamis in [18], the inverse scattering transform (IST) has been used for the Korteweg–de Vries (KdV) in [19], Padé approximants have been used to reduce the Gibbs phenomenon in [20], series expansions in the vertical excursions of the interface and bottom topography are used in the three-dimensional modelling of tidally driven internal wave formation in [21]).

The aim of this paper is to develop an efficient semi-analytical numerical approach based on the weakly nonlinear analysis of the problem, which could be used to model the scattering of nonlinear waves on the defects. The developed approach can be applied to other problems of this type. For example, a system of coupled Boussinesq-type equations was derived to describe longitudinal waves in a layered waveguide with soft bonding [11]. In [22], we constructed a hierarchy of weakly nonlinear solutions of the initial-value problem for localized initial data in terms of solutions of Ostrovsky-type equations [23] for unidirectional waves, which opens the way for the application of the method developed in this paper to the study of the scattering of nonlinear waves in other delaminated layered structures.

The problem of calculating the reflected and transmitted waves from a wave incident on an interface between two media has been discussed in several areas. In fluids, it is known that when a surface soliton passes through an area of rapid depth variation, an incident solitary wave evolves into a group of solitons [24–26]. There is also a similar effect for internal waves [27–29]. The problem of collimated light beams incident on an interface separating two linear dielectric media has been extensively studied [30–32]. This can be extended to nonlinear dielectric media and theoretical results have been found [33]. Interface Sturm–Liouville systems have been studied in [34] and applications of these systems have been discussed in [35,36]. It is worth noting that the systems considered in these examples often have simple boundary conditions, while our problem takes account of continuity of stress.

The paper is organized as follows. In §2, we briefly overview the problem formulation for the scattering of an incident soliton in a symmetric layered bar with delamination, and the weakly nonlinear solution of the problem [10]. The approach gives rise to three KdV equations describing the behaviour of the leading order incident, transmitted and reflected waves, and terms for the higher order corrections. We describe fission of the transmitted and reflected strain solitary waves and establish predictions for the number of solitary waves present in each section of the bar. In §3, we present two numerical schemes. In §3b we describe the finite-difference approximations for the problem, with a ‘kink’ as the initial condition for the displacement (corresponding to the strain soliton). The system is treated as two boundary value problems and the finite-difference approximation of the system produces two tridiagonal matrices with a nonlinear condition at *x*=0 (the common point in both boundary value problems). In §3b, we introduce a semi-analytical numerical scheme, which is based upon the weakly nonlinear solution discussed in §2. The three KdV equations are solved numerically using a strong stability preserving Runge–Kutta (SSPRK) scheme proposed in [37]. In §4, we compare the results of the numerical schemes from §3a,b, against each other and the predicted results, and the agreement is checked for various configurations of the delaminated bar. Finally, we conclude our findings in §5 and discuss potential extensions of this work, pertaining to other types of bonding of layered structures, described by coupled regularized Boussinesq (crB) equations [11].

## 2. Weakly nonlinear solution

We consider the scattering of a long longitudinal strain solitary wave in a perfectly bonded two-layered bar with delamination at *x*>0 (figure 1). The material of the layers is assumed to be the same (symmetric bar), while the material to the left and to the right of the *x*=0 cross-section can be different. The problem is described by the following set of regularized non-dimensional equations [10]
*c*,*α* and *β* are constants defined by the geometrical and physical parameters of the structure, while *ϵ* is the small wave amplitude parameter. The functions *u*^{−}(*x*,*t*) and *u*^{+}(*x*,*t*) describe displacements in the bonded and delaminated areas of the structure, respectively. Condition (2.3) is continuity of longitudinal displacement, while condition (2.4) is the continuity of stress. In order to reduce the number of parameters in our numerical experiments, in what follows we will use the values presented in [10], namely that *α*=1 and
*n* represents the number of layers in the structure and *k* is defined by the geometry of the waveguide. Referring to figure 1, the cross section *x*=0 has width 2*a* and the height of each layer is 2*b*/*n*. In terms of these values, *k*=*b*/*a* and, as there are two layers in this example, *n*=2. In §4, we consider various configurations of the bar.

### (a) Derivation of KdV equations

Let us consider a weakly nonlinear solution to (2.1) of the form [10]
*I*(*ξ*_{−},*X*), *R*(*η*_{−},*X*) and *T*(*ξ*_{+},*X*) describe the leading order incident, reflected and transmitted waves, respectively, while the functions *P*(*ξ*_{−},*η*_{−},*X*) and *Q*(*ξ*_{+},*η*_{+},*X*) describe the higher order corrections. Substituting this weakly nonlinear solution into (2.1) the system is satisfied at leading order, while at *O*(*ϵ*) we have
*ϕ*(*ξ*_{−},*X*), *ψ*(*η*_{−},*X*) are arbitrary functions. The first radiation condition requires that the solution to the problem should not change the incident wave (see [10] for the discussion of the radiation conditions). In our problem, the incident wave is the solitary wave solution of the KdV equation (2.8), with corrections at *O*(*ϵ*^{2}). Therefore, the radiation condition implies that *ϕ*(*ξ*_{−},*X*)=0. We find *ψ*(*η*_{−},*X*) from conditions (2.3) and (2.4) later. A similar calculation for the second equation in (2.1) gives

We are looking for the leading order transmitted wave satisfying
*q*(*ξ*_{+},*X*), *r*(*η*_{+},*X*) 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. Thus, *r*(*η*_{+},*X*)=0. From (2.3), we see
*O*(*ϵ*)
*O*(*ϵ*)
*x*=0 in terms of the incident wave, i.e.
*c*=1), then *C*_{T}=1, *C*_{R}=0 and there is no leading order reflected wave.

We can now simplify (2.17) and (2.19) using the KdV equations (2.8), (2.10), (2.14) and relations (2.20) to obtain

Considering (2.17) and (2.19), we have
*f*,*g* are now completely defined in terms of the leading order incident, reflected and transmitted waves. We restore the dependence of *f* and *g* on their respective characteristic variables to obtain

### (b) Fission

We can rewrite the transmitted wave equation in the canonical form
*x*=0 cross section. We define the transmitted wave field by the spectrum of the Schrödinger equation
*A* will determine whether any solitary waves are present in the transmitted wave field. If *A*<0, the transmitted wave field will not contain any solitons and the initial pulse will degenerate into a dispersive wave train. However when *A*>0, there will always be at least one discrete eigenvalue, corresponding to at least one solitary wave in the transmitted wave field, and accompanying radiation determined by the continuous spectrum. It is worth noting that, in some cases, more than one secondary soliton can be produced from the single initial soliton, referred to as fission of a soliton [24–26].

The discrete eigenvalues of (2.30) take the form (e.g. [39])
*N*, is given by the largest integer satisfying the inequality
*β* and *c* depend on the properties of the material and the geometry of the waveguide. We can see from (2.33) that, for small *δ*, there will always be one soliton while, as *δ* increases, more solitons will emerge. As *χ*_{n} is the phase shift.

Is soliton fission possible if the waveguide is made of one and the same material? This requires setting *c*=1. We find that fission can occur and is dependent upon the geometry of the waveguide. Recall the expression for *β* given in (2.5), where *n* is the number of layers and *k*=*b*/*a*. The number of solitary waves produced is the largest integer satisfying
*β* and these will be checked in §4.

A similar description can be made for the reflected wave, which is already written in canonical form in (2.10), making use of the ‘initial condition’ and reflection coefficient as presented in (2.20) and (2.21), respectively. The wave field is defined by the spectrum of the Schrödinger equation (2.30), where the potential *U* is given by
*B* is dependent upon the sign of the reflection coefficient *C*_{R}. If *c*<1, then *B* is negative and the reflected wave field does not contain any solitary waves. The initial pulse will degenerate into a dispersive wave train. For *c*>1, *B* is positive and there will be at least one solitary wave present in the reflected wave field, accompanied by radiation. The solitary waves can be described using formulae (2.32) and (2.33) and making the change

If *c*=1 (the structure is of one and the same material), then *C*_{R}=0 and there is no leading order reflected wave.

## 3. Numerical schemes

To numerically solve equations (2.1), with continuity conditions (2.3) and (2.4), we consider two numerical schemes. The first builds on the scheme for a single Boussinesq equation [22], treating the problem as two boundary value problems (BVPs) with the continuity conditions defining the interaction between the two problems. The second approach is a semi-analytical method that solves the equations derived in §2 using an SSPRK scheme (see [37,40,41] for method and [42] for examples).

### (a) Direct finite-difference scheme

Firstly, we consider the Boussinesq-type equations posed in §2. Discretizing the domain [−*L*,*L*]×[0,*T*] into a grid with equal spacings *h*=Δ*x* and *κ*=Δ*t*, the analytical solution *u*^{±}(*x*,*t*) is approximated by the exact solution of the finite difference scheme *u*^{±}(*ih*,*jκ*), denoted

As we are considering localized initial data for strains, if we take *L* large enough then we can enforce zero boundary conditions for the strain, i.e. *u*_{x}=0. Therefore, applying a central difference approximation to this condition, we have
*i*=0,…,*N* and *i*=*N*,…,2*N*. Denoting the right-hand side of (3.1) as *f*_{i} and the right-hand side of (3.2) as *g*_{i}, we make a rearrangement around the central boundary, namely
*i*=0,…,*N* can be written in matrix form as
*i*=*N*,…,2*N*. This system is tridiagonal and therefore we can use the Thomas algorithm (e.g. [43]) to solve both tridiagonal systems in terms of *u* of the form
*N* large enough, we can set *h*_{0}≡0. This leads to the result

To determine the initial condition, we differentiate (2.1) and denote *ε*(*x*,*t*), and therefore the exact ‘kink’ solution for (2.1). Explicitly we have
*v*_{1} is the velocity of the wave. 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).

### (b) Semi-analytical approach

To solve equations (2.8), (2.10) and (2.14) numerically, we use the SSPRK(5,4) scheme as described in [37] (see [42] for example of this method). A hybrid Runge–Kutta scheme was used initially (e.g. [44] and the methods used for the finite-difference terms in the equation were based upon [45,46]). However, it was found the SSPRK(5,4) scheme had a higher accuracy, owing to its stability preserving properties, and the computational time was similar for both schemes. In addition, the SSPRK(5,4) scheme can be adapted more easily.

We discretize the domain [−*L*,*L*]×[0,*T*] into a grid with equal spacings *h*=Δ*ξ*, *κ*=Δ*X*, and the analytical solution *e*^{±}(*x*,*t*) is approximated by the exact solution of the numerical scheme, *e*^{±}(*ih*,*jκ*) (denoted

The SSPRK(5,4) scheme is as follows. Let *e*(*x*,*t*) denote a solution of the KdV equation. Given that the solution at time *t*_{j}=*jκ* is given by
*t*_{j+1}=(*j*+1)*κ* is given by
*F* is the finite-differenced form of all the terms in the KdV equation involving spatial derivatives. Note that the coefficients here are chosen in such a way to optimize the time step at each point and, due to the complexity of each coefficient, are presented to 15 decimal places.

To obtain an expression for *F*, central difference approximations are applied for the first and third derivatives and an average is taken for the nonlinear term (as was performed by Zabusky & Kruskal, see [1,47]). We also assume boundary conditions of the form
*v* is the velocity of the soliton and can be related to the velocity of the solution (3.14) by the formula

## 4. Results

In what follows, we let *ϵ*=0.01,*v*=0.35 and take a step size of *h*=*κ*=0.01 in the finite-difference scheme and a step size of *h*=0.1, *κ*=0.001 in the SSPRK(5,4) scheme. The value of *h* and *κ* chosen for the SSPRK(5,4) scheme was chosen as the same step size taken for the hybrid Runge–Kutta scheme in [44]. Subsequent numerical experiments have shown that the SSPRK(5,4) scheme is not sensitive to considerable changes in the step sizes, and the results for these chosen step sizes and much larger values of *h*=0.25, *κ*=0.01 produce comparable results with an error of 6.5×10^{−3}. As shown in [22], the finite-difference method for the Boussinesq equation is linearly stable (using a von Neumann linear stability analysis) for values of *κ* satisfying
*f*_{0} is the constant used in the linearized scheme. In practice, the stricter condition of
*h* and *κ* chosen above satisfy this relation.

The results are plotted for the strain waves, to correspond with the prescription of the initial condition. Therefore, we denote

### (a) Test cases

Firstly, we consider the case where *β*=*c*=1. The initial condition for the Boussinesq equation (2.1) is given by (3.14). We solve the equation for the displacement *u*^{+}(*x*,*t*) and plot the strain *u*^{+}_{x}(*x*,*t*). The initial condition for the KdV equation (2.14) is provided by (3.19). Viewing the solution of the first scheme as exact, the solution of the second scheme will be accurate to leading order, with a small correction to the wave at *O*(*ϵ*^{2}). The solution for *e*^{−} at the initial time and for *e*^{+} at a sufficiently large time should describe the same right-propagating solitary wave. The result of the calculation for this case, for both numerical schemes, is presented in figure 2 for the transmitted wave. As *c*=1, there is no leading order reflected wave. Figure 2 shows a good agreement between the solutions, with a small phase shift between the two cases. This can be remedied by adding higher order terms to the solution of the second scheme.

Another test case is for a value of *c* such that *c*<1. Considering the expressions for *C*_{R} and *C*_{T}, that is
*c*>0 (a physical requirement), we observe that *C*_{T}>0 for all *c*, while *C*_{R}<0 for *c*<1 and *C*_{R}>0 for *c*>1. Recalling our observations from §2 we expect that, for a value of *c*<1, a dispersive wave train will be present in the reflected wave field. For *c*=0.25, we note that

The next case to consider is for large values of *c*. Using the results of §2, specifically (2.21) and (2.22), we would expect the leading order reflected wave to be closer in amplitude to the initial wave and the transmitted wave to be of a much smaller amplitude. An example is presented in figures 4 and 5 for *c*=5, where we have

### (b) Predictions

Following the scheme outlined in §2b, we define *c*=1 and recall that
*n* and *k* to obtain different cases. Table 1 in [10] presents results for nine different configurations and all of these configurations have been checked. For brevity, we only present results for four cases. The expected number of solitary waves in the transmitted wave field in these cases are summarized in table 1. The results are presented in figures 6–9. In each of the following cases, we note that there is no leading order reflected wave, as *c*=1, and therefore only the transmitted wave field is presented. We observe that the number of predicted solitons is in agreement with the numerics in all cases. In the cases where the smaller solitons are close in amplitude to the radiation, a further simulation was run for a larger time to ensure that the smaller soliton moves away from the radiation generated by the continuous spectrum. The figures are all presented at the same time value of *t*=1000 for comparison.

### (c) Predicted heights

A prediction for the amplitude of the lead soliton is provided in [10], namely that the ratio of its amplitude to that of the incident soliton is
*κ*=1,2,3 and *n*=2,3,4 for the following results. Taking *h*=*k*=0.01 in the finite-difference scheme, we have a maximum error of 5.877×10^{−3} and for step sizes *h*=0.1,*k*=0.001 in the SSPRK(5,4) scheme we have a maximum error of 1.461×10^{−4}.

### (d) Comparison of schemes

Reviewing the results presented here, it can be seen that the semi-analytical approach produces results comparable to the direct finite-difference scheme for many cases, particularly for long waves. As the model is a long-wave model, this is the desirable behaviour for the schemes. Crucially, the semi-analytical approach requires the solving of, at most, two equations in each section of a bar (reflected and transmitted) while the direct finite-difference method requires the solution of multiple tridiagonal equations systems, and the solution of the nonlinear equation at *x*=0 has to be substituted into the implicit solution at all other points. Therefore, the semi-analytical approach is more desirable as additional sections are included in the bar. In addition, an analysis of the time required to calculate a solution by each scheme shows that the semi-analytical approach is faster.^{1} The accuracy of this method can be improved further by adding higher order terms [49,50]. It must be noted that, in these papers, we have developed a semi-analytical numerical approach to the solution of the initial value problem for Boussinesq-type equations; however the advantages of the semi-analytical method compared to the standard methods were not obvious. They became apparent for the type of problems discussed in this paper.

## 5. Conclusion

In this paper, we considered two numerical schemes for solving two boundary-value problems matched at *x*=0. This problem represents the propagation of a strain wave in a bar with a perfect bonding between the two layers for the first boundary-value problem, and a bar with complete delamination in the second boundary-value problem.

At this point, we looked for a weakly nonlinear solution by taking a multiple-scales expansion in terms of the appropriate set of fast and slow variables. This produced three KdV equations satisfying initial-value problems, describing the leading order incident, reflected and transmitted waves. The initial values in these equations are fully described in terms of the leading order incident wave, with reflection and transmission coefficients being derived to describe these initial conditions. It was noted that a bar made of one and the same material will not generate a leading order reflected wave. In addition, expressions were found for the higher order corrections in terms of the leading order incident and reflected waves, and the reflection and transmission coefficients.

Following [10], the process of fission of an incident solitary wave was then discussed and expressions for the number of solitons in the delaminated section and their eigenvalues were found. This was applied to the case when the materials in the perfectly bonded and delaminated sections of the bar are one and the same, and a formula for the number of solitons present in the delaminated section was found in terms of the geometry of the waveguide. This allowed for predictions to be made to the number of solitons for a given configuration of the bar. A similar description was given for the reflected wave and a similar prediction can be made in this case.

Two numerical schemes were suggested to describe this behaviour, taking account of the original derivation and the weakly nonlinear approach. For the original equations, we made use of finite-difference techniques, which resulted in two tridiagonal systems and a nonlinear difference equation linking the systems. These systems were solved, in terms of ‘ghost points’ in both systems, using a Thomas algorithm [43] and the result of this calculation was used to find the solution of the nonlinear equation for one of the ghost points. The solution for this ghost point is then substituted back into the implicit solution of the tridiagonal system to determine the solution at a given time value. The result allowed us to analyse the leading order transmitted and reflected waves in their respective domains and compare the results to known analytical predictions. The numerical modelling has confirmed that the incident soliton fissions in the delaminated area, which can be used to detect the defect.

The second scheme solves the derived KdV equations in terms of the incident strain wave. This semi-analytical approach makes use of a SSPRK(5,4) scheme, with the finite-differenced form of the spatial derivatives used as the function in the scheme [37]. This scheme was used to calculate the incident, transmitted and reflected waves from their respective equations. The reflection and transmission coefficients determine the magnitude of the initial condition in the equations for the reflected and transmitted waves, respectively.

These results were then compared to each other and to known analytical results. It was found that the schemes were in agreement to leading order, with a slight phase shift between the two solutions representing the *O*(*ϵ*^{2}) difference between their propagation speeds. Furthermore, the numerics were in agreement with the predicted theory for the number of solitons in the delaminated section for several choices of waveguide geometry. In addition, it was found that the semi-analytical approach has a smaller computation time than the direct method, and is simpler to implement in comparison to the direct finite-difference approach. This effect will be more dramatic as the number of equations is increased (corresponding to more sections in the bar). Finally, the methods used in this paper can be generalized to a bar with a soft bonding layer, which leads to a system of coupled Boussinesq equations for *x*<0 and uncoupled equations for *x*>0 [11].

## Data accessibility

This work does not have any experimental data. Calculations were performed using the C++ compiler, v. 15.0.0.

## Authors' contributions

All authors contributed to the analysis, computations and the writing of this paper. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

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. The paper is published under the Creative Commons Attribution licence.

## Acknowledgements

The authors would like to thank Noel Smyth and the referees for useful discussions of the finite-difference approach.

## Footnotes

↵1 Both methods were programmed in C and compiled using the Intel compiler. The spatial domain for the semi-analytical method used 75 000 points, while the spatial domain for the finite-difference method used 400 000 points. When run on a 2.66 GHz Intel Core i5 machine, the approximate CPU time for the direct finite-difference scheme was 12 h, compared to 4 h for the semi-analytical method. Each scheme was parallelized and therefore the actual calculation time was reduced by a factor based upon the number of cores in the processor. Subsequent numerical experiments have shown that in the simple settings considered in this paper the domain size can be further reduced to improve the calculation time, so we can consider a spatial domain containing 15 000 points for the semi-analytical method, and 150 000 points for the finite-difference method. This reduced the calculation time to approximately 20 min for the finite-difference scheme and 15 min for the semi-analytical scheme. However, as the semi-analytical scheme is not sensitive to considerable increases in the step sizes, we can take step sizes of

*h*=0.25 and*κ*=0.01 to obtain a calculation time of 34 s. Therefore, the resulting calculation in the semi-analytical scheme was 35 times faster than the finite-difference scheme.

- Received August 21, 2015.
- Accepted October 12, 2015.

- © 2015 The Authors.