Modelling of nonlinear wave scattering in a delaminated elastic bar.

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.

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

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] with appropriate initial conditions and associated continuity conditions and where 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 where 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 2a and the height of each layer is 2b/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.

Figure 1.

Two-layered symmetric bar with delamination at x>0. (Online version in colour.)

Derivation of KdV equations

Let us consider a weakly nonlinear solution to (2.1) of the form [10] where the characteristic variables are given by Here, the functions 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 To leading order the right-propagating incident wave is defined by the solution of the KdV equation Similarly, the reflected wave satisfies the KdV equation Substituting these conditions into (2.6) and integrating with respect to both characteristic variables, we obtain an expression for higher order terms satisfying where ϕ(ξ−,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 where and higher order corrections are given by where 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 so to leading order we have and at O(ϵ) Similarly for (2.4) we have, to leading order and at O(ϵ) Using (2.16) and (2.18), we can derive initial conditions for the previously derived KdV equations, defining both reflected and transmitted ‘strain’ waves at x=0 in terms of the incident wave, i.e. where and are the reflection and transmission coefficients, respectively. It is worth noting that, if the materials are the same in the bonded and delaminated areas (c=1), then CT=1, CR=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 and

Considering (2.17) and (2.19), we have and The functions 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 and The constants of integration in the KdV equations and (2.27), (2.28) should be found from additional physical conditions.

Fission

We can rewrite the transmitted wave equation in the canonical form via the change of variables The method of the IST [38] can be used to determine the solution for the KdV equation. We make a similar approximation to what was done in [26], for a soliton moving into a region with different properties, neglecting some short waves as the soliton moves over the x=0 cross section. We define the transmitted wave field by the spectrum of the Schrödinger equation where the potential is given by The sign of 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]) The number of solitary waves produced in the delaminated area, N, is given by the largest integer satisfying the inequality where In the above, the parameters β 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 , the solution will evolve into a train of solitary waves, ordered by their heights, propagating to the right and some dispersive radiation (a dispersive wave train) propagating to the left (in the moving reference frame), i.e. where χ 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 This gives rise to a series of predictions based upon the value of β 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 where we have used (2.21). It is clear that the sign of B is dependent upon the sign of the reflection coefficient CR. 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 , and .

If c=1 (the structure is of one and the same material), then CR=0 and there is no leading order reflected wave.

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

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 . We make use of first-order and second-order central difference approximations in the main equations, namely and the approximation for can be obtained iteratively using the approximations for and . To simplify the obtained expressions, we introduce the notation Substituting these into system (2.1) gives and Assuming the domain can be discretized, we denote and therefore continuity condition (2.3) becomes In the continuity condition (2.4), we make use of the central difference approximations presented above, and introduce ‘ghost points’ of the form and . Therefore, (2.4) becomes

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=0. Therefore, applying a central difference approximation to this condition, we have and we now solve for the boundary points using these relations. We note that (3.1) and (3.2) form tridiagonal systems, specifically for the values i=0,…,N and i=N,…,2N. Denoting the right-hand side of (3.1) as f and the right-hand side of (3.2) as g, we make a rearrangement around the central boundary, namely The equation system for i=0,…,N can be written in matrix form as with a similar system for i=N,…,2N. This system is tridiagonal and therefore we can use the Thomas algorithm (e.g. [43]) to solve both tridiagonal systems in terms of and , respectively. This intermediary solution is then substituted into (3.5) to obtain a nonlinear equation in terms of one of the ghost points, which can then be used to determine the solution for this time step. To obtain this expression, we denote where we calculate the coefficients , , and from (3.8) and the associated system for the second equation. To solve this system in terms of the ghost point , we firstly express the variables and in terms of this ghost point. Therefore, denoting for brevity, we make use of (3.4) and obtain Substituting (3.9) and (3.10) into (3.5), we can find a quadratic equation for u of the form where we have In order to simplify the above equations, we need to find expressions for the coefficients in (3.9), so we consider the intermediary steps in the Thomas algorithm to obtain explicit representations. By considering the intermediary steps, we can deduce that, for N large enough, we can set and (this has been confirmed by numerical calculation) and therefore h0≡0. This leads to the result and therefore we can determine , and similarly and . These values can be substituted into the implicit solution of the tridiagonal system to determine the solution at each time step. It is worth noting that this scheme is second order.

To determine the initial condition, we differentiate (2.1) and denote . The initial condition is then taken as the exact solitary wave solution of this equation for ε(x,t), and therefore the exact ‘kink’ solution for (2.1). Explicitly we have where v1 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).

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 ), where is the sum of the numerical solutions to (2.8) and (2.10) and is the numerical solution to (2.14).

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κ is given by then the solution at t=(j+1)κ is given by where the function 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 Using this scheme in equation (2.8) (which is identical to the equation (2.10) in characteristic variables), we obtain and for (2.14) we obtain The initial condition is taken as the exact solution of the KdV equation (2.8), namely where v is the velocity of the soliton and can be related to the velocity of the solution (3.14) by the formula . For equation (2.14), a multiplicative factor is applied to the initial condition in (3.18), namely that which we derived in (2.22). Explicitly, we have Similarly, for equation (2.10) a different multiplicative factor is applied to the initial condition (3.18) of the form (2.21) and therefore we have

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 where f0 is the constant used in the linearized scheme. In practice, the stricter condition of is imposed, to help accommodate for the nonlinearity effects. The values 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 and .

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

Figure 2.

The solution e+ at t=1000 in equations (3.2) and (3.17) with exact initial conditions and initial position x=−50. (Online version in colour.)

Another test case is for a value of c such that c<1. Considering the expressions for CR and CT, that is and noting that c>0 (a physical requirement), we observe that CT>0 for all c, while CR<0 for c<1 and CR>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 and the result for this calculation is presented in figure 3. It can be seen that a dispersive wave train is present in both numerical schemes. The KdV approximation overestimates the amplitude of the wave train, but correctly resolves all its main features. It is worth noting that this solution looks like the Airy function solution of the linearized KdV equation, which is the small amplitude limit of the similarity solution of the KdV equation related to the Painleve II equation (see [4,48] and references therein for more details).

Figure 3.

The solution e− at t=1000 in equations (3.1) and (3.16) with c=0.25 and an initial position of x=−50. (Online version in colour.)

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 and we can see from the coefficients that we would expect the reflected wave to have a larger amplitude than the transmitted wave, by approximately one order of magnitude in this case.

Figure 4.

The solution e− at t=1000 in equations (3.1) and (3.16) with c=5 and initial position x=−50. (Online version in colour.)

Figure 5.

The solution e+ at t=1000 in equations (3.2) and (3.17) with c=5 and initial position x=−50. (Online version in colour.)

Predictions

Following the scheme outlined in §2b, we define c=1 and recall that We now consider the behaviour in our problem formulation using these parameters. The initial condition remains the same and we change the value of 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.

Table 1.

Predictions on the number of solitons present in the transmitted wave field in the delaminated section of the bar, for various choices of n and k.

no. layers, n	ratio of height and half width, k	value of dispersion coefficient, β	no. solitons
2	1	58	2
3	1	59	2
3	3	15	3
4	2	14	3

Figure 6.

The solution e+ at t=1000 in equations (3.2) and (3.17) with , corresponding to n=2 and k=1, and initial position x=−50. (Online version in colour.)

Figure 9.

The solution e+ at t=1000 in equations (3.2) and (3.17) with , corresponding to n=4 and k=2, and initial position x=−50. (Online version in colour.)

The solution e+ at t=1000 in equations (3.2) and (3.17) with , corresponding to n=3 and k=1, and initial position x=−50. (Online version in colour.)

The solution e+ at t=1000 in equations (3.2) and (3.17) with , corresponding to n=3 and k=3, and initial position x=−50. (Online version in colour.)

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 We consider the combination of values κ=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.

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.

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