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][3][4] and references therein) it was also shown, within the framework of the nonlinear dynamic elasticity theory, that Boussinesqtype 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][8][9].

Weakly nonlinear solution
We consider the scattering of a long longitudinal strain solitary wave in a perfectly bonded twolayered 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 u ± (x, 0) = F ± (x), (2.2) and associated continuity conditions , (2.4) 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 β = n 2 + k 2 n 2 (1 + k 2 ) , (2.5) 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.
To leading order the right-propagating incident wave is defined by the solution of the KdV equatioñ Similarly, the reflected wave satisfies the KdV equationR 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 We are looking for the leading order transmitted wave satisfying 14) and higher order corrections are given by so to leading order we have (2.16) and at O( ) Similarly for (2.4) we have, to leading order 18) 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. and , (2.22) 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 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 (2.23) and  (2.25) 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 ψ(η − , 2.27) and

(b) 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][25][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 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.
U ∼ − N n=1 2k 2 n sech 2 (k n (χ − 4k 2 n τ − χ n )) + radiation, (2.34) where χ 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 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 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 A → B, and l → m.
If c = 1 (the structure is of one and the same material), then C R = 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).

(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 u ± i,j . We make use of first-order and second-order central difference approximations in the main equations, namely and the approximation for u ± ttxx can be obtained iteratively using the approximations for u ± tt and u ± xx . To simplify the obtained expressions, we introduce the notation Substituting these into system (2.1) gives 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 u − N+1, j+1 and u + N−1, j+1 . 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 x = 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 righthand 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, namelỹ The equation system for i = 0, . . . , N can be written in matrix form as 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 φ ± 1 , ψ ± 1 , φ ± 2 and ψ ± 2 from (3.8) and the associated system for the second equation. To solve this system in terms of the ghost point u − N+1, j+1 , we firstly express the variables u + N−1, j+1 and u + N+1, j+1 in terms of this ghost point. Therefore, denoting u = u −

N+1, j+1
for brevity, we make use of ( 3.4) and obtain where we have and 12) 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 ψ − 1 = ψ + 1 and ψ − 2 = ψ + 2 (this has been confirmed by numerical calculation) and therefore h 0 ≡ 0. This leads to the result and therefore we can determine u − N,j+1 , u + N,j+1 and similarly u − N−1, j+1 and u + N+1, j+1 . 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 where 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 e ± i,j ), where e − i,j is the sum of the numerical solutions to (2.8) and (2.10) and e + i,j is the numerical solution to (2.14).
Using this scheme in equation (2.8) (which is identical to the equation (2.10) in characteristic variables), we obtain ( 3.16) and for (2.14) we obtain (3.17) 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 v 1 = √ 1 + 2 v. 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 f 0 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 e − = u − x and e + = u + x .  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.

(a) Test cases
Another test case is for a value of c such that c < 1. Considering the expressions for C R and C T , that is and noting that 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 C R = − 3 5 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). 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 C R = 2 3 and C T = 1 15 , 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.
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

(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 We

(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 semianalytical 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 rspa.royalsocietypublishing.org Proc. R. Soc  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; 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 semianalytical 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.
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].