Reflection of waves in a waveguide from a boundary with nonlinear stiffness: application to axial and flexural vibrations

The reflection of time-harmonic waves in a waveguide with a nonlinear boundary stiffness is considered with applications to rods and beams. Incident waves at frequencies that are multiples of a fundamental frequency give rise to reflected propagating and nearfield waves at the same frequencies. An infinite set of equations is developed for the reflection coefficients, which depend on the amplitudes and phases of the incident waves. Nonlinear boundary conditions are applied, and equations is truncated by using the harmonic balance method and solved numerically. The case of zero linear boundary stiffness, i.e. essential nonlinearity, is studied. First, the case where there is only one incident wave is considered. An approximate solution is found when retaining only two reflected waves. Numerical examples are presented, energy being seen to leak into the higher harmonics. The minimum magnitudes of the reflection coefficients of axial and flexural vibrational waves at the fundamental frequency and the maximum energy that can leak into the higher harmonics are determined. Accuracy and convergence when retaining different numbers of reflected harmonics are illustrated. The case of two incident waves is then considered. Multiple incident waves affect the leakage of energy to higher harmonics and can have a significant effect on the reflection coefficient for the fundamental harmonic. With some parameters, a much lower reflection coefficient is obtained for the wave at the fundamental frequency as compared to the case of one incident wave. It is seen that with two incident flexural waves, the reflection coefficients can be multi-valued for certain values of the system parameters. A numerical study is performed to show the region of multiple solutions.


Introduction
Exploiting nonlinearity to control vibrations in elastic structures such as rods and beams is of interest to academia and industry. This paper concerns the reflection of time-harmonic axial and flexural waves in a thin rod and an Euler-Bernoulli beam from a boundary with nonlinear stiffness. The nonlinearity typically scatters energy into higher frequencies; the energy then decays more rapidly.
Many previous papers have concerned the effects of a linear boundary condition on wave reflection in uniform and non-uniform continuous structures such as rods, plates, Euler-Bernoulli beams and Timoshenko beams (e.g. [1][2][3][4][5][6][7]). The wave propagation approach has also been widely used to study vibrations of nanorods [8] and nanobeams [9] with cracks. In these studies, the crack has been modelled as a linear spring element.
Some studies concerned initial value problems and free wave propagation and developed solutions based on the d'Alembert's method [10][11][12][13][14]. This paper concerns steady-state, time-harmonic axial and flexural waves in a rod or a beam with a nonlinear boundary stiffness using a wave propagation approach. Nonlinear dynamic systems typically have more involved behaviour than their linear counterparts. The source of nonlinearity might be in the material properties of the system or due to large deformations (geometrical nonlinearities) or at boundaries. Many studies concern vibrations of elastic structures with material and geometrical nonlinearities using a wave approach, e.g. [15].
Different analytical and numerical approaches can be used to study nonlinear vibrations in rods and beams. Santo et al. [16] studied the vibrational behaviour of a rod with multiple springs with both linear and nonlinear stiffness using perturbation and polynomial methods.
Wave propagation approaches have also been used to model the dynamics. Nayfeh et al. [17] and Vakakis [18] studied the axial vibration of a rod with a Duffingtype nonlinear boundary stiffness through wave scattering. A wave incident on a boundary can be scattered into an infinite number of reflected waves with higher frequencies. It was shown that the physical behaviour of the system can be described by considering only the lower harmonics. In these studies, the nonlinearity was weak, and the strength of nonlinearity was defined in terms of the linear force of the boundary spring compared to the nonlinear force, so that the case of essential nonlinearity, i.e. when the linear spring force is zero, was not investigated. Although the power transferring from the first harmonic to higher harmonics was studied, the maximum energy that can leak into higher harmonics was not investigated. Furthermore, the results were shown for the case when only one wave is incident on the boundary. The purpose of this paper is to extend these works for the cases of a beam or a rod in a more general form, i.e. by including incident and reflected waves with higher harmonics and studying the case of essential nonlinearity. The higher harmonic incident waves might, for example, arise from the reflection of the negative-going waves from a boundary or discontinuity elsewhere on the waveguide.
Wave propagation in beams is more complex than in rods since nearfield waves also exist. Nearfield or evanescent waves have amplitudes that decay exponentially with distance and do not propagate energy along the waveguide [19]. Following the analysis presented in [17,18], Brennan et al.
[20] studied hardening and softening nonlinear boundaries for rods and beams. However, in this study, the effects of incident and reflected waves of higher harmonics were neglected: only a single incident and reflected wave was considered. It was shown that, for hardening boundary stiffness, when retaining only the first harmonic, the magnitude of the reflected wave is identical to that of the incident wave, but there is a phase change in the reflected wave. This conclusion also follows from the conservation of energy. However, with softening stiffness and when a nearfield wave is incident upon the boundary, the displacement at the end of the beam was seen to jump at certain frequencies, and the solution becomes multi-valued. No jump was found in the displacement at the end of the rod and beam with hardening boundary stiffness.
Tang et al. [21] followed the previous study [20] and calculated the natural frequencies of a rod and a beam by using a phase closure principle. It was shown that softening-type stiffness substantially affects the phases of the reflection coefficients with instabilities in a certain range of frequencies. In [21,22], the nonlinearity was also weak, and the case of essential nonlinearity was not considered. The present paper can also be considered as an extension in the methodology and results of works [20,21] by retaining incident and reflected waves with higher harmonics and investigating the case of essential nonlinearity.
Other wave approaches have also been used to study the vibration of nonlinear systems. In [22,23] a wave finite element method was used to analyse composite structures with damage. The damage was modelled by a nonlinear coupling element, and only the first incident wave was considered; however, reflected waves with higher frequencies were retained. Chouvion [24] used a different wave approach and worked on nonlinear vibrations of beam-like structures. The principle behind this matrix-based approach was the ray-tracing method. In Chouvion's approach, reflected and incident waves matrices were coupled and connected at discontinuities. The method was applied to a rod and a turbine blade (beam-like network) with damping. Essential nonlinear stiffness boundaries have a wide range of physical applications as they can absorb the system's vibrational energy in a one-way and irreversible process. The nonlinear energy sink (NES) is a well-known application of such systems [25]. Chouvion [26] applied the same approach as in [24] to a beam with an NES in a subsequent study. The beam was connected to a linear torsional spring at each end, and the nonlinearity of the spring component of the NES was assumed to be cubic. The frequency response functions of the system were determined by the retention of waves with higher harmonics. However, the focus in [24] and [26] was on the time-harmonic forced response of the systems, rather than the calculation of the reflection coefficients, and the effects of the magnitude and phase of incident waves with higher harmonics on the reflection coefficients were not studied. Furthermore, the maximum leakage of energy from the first harmonic to the higher harmonics was not investigated.
In the present paper, the same wave approach as in [17,18] and [20,21] is used to study the nonlinear vibration of a semi-infinite, uniform, thin rod and an Euler-Bernoulli beam. Incident, reflected and reflected nearfield waves with multiple frequencies are considered. Under the assumption of cubic nonlinearity, a set of equations is developed for the reflection coefficients. The infinite set of equations is then truncated at an arbitrary number of harmonics and solved by using the harmonic balance method. The case of essential nonlinearity, for which the effects of nonlinearity are strongest, is investigated. For the case of one incident wave and essential nonlinearity, reflected waves with higher frequencies are retained to study energy transfer and the maximum energy that can leak from the first harmonic to higher harmonics.
It is also shown that multiple incident waves have significant effects on the reflection coefficients and can decrease or increase their magnitudes. Numerical examples are presented, and the conditions for which the reflection coefficient of the first harmonic is minimum are indicated. Furthermore, it is seen that, for flexural vibration, with two incident waves, at certain frequencies the response becomes multivalued, and there are at least two solutions for the reflection coefficients. The region in which jumps occur in the magnitudes of the reflection coefficients is determined.
The paper is organized as follows: In Sect. 2, the governing equations are derived, and reflection coefficients are defined. The nondimensional parameters and strength of nonlinearity are defined regarding a reference frequency, and a set of equations is developed in terms of the reflection coefficients for the case of cubic nonlinearity. In Sect. 3, equations are truncated using the harmonic balance method for the case when there is only one incident wave and an arbitrary number of reflected waves. In Sect. 4, an analytical solution based on the method of expansion in a small parameter is found, and expressions are derived for the reflection coefficients when one incident and two reflected waves are retained. Section 5 presents the equations for the case of multiple incident and reflected waves. Numerical results are presented in Sect. 6, and finally, Sect. 7 is the conclusion. The focus is primarily on flexural waves, but results for the simpler case of axial waves are also presented.

Problem description and equations of motion
We consider the dynamic behaviour of a semi-infinite waveguide lying along the x-axis. The boundary at x = 0 is attached to a spring with nonlinear stiffness. Figure 1 shows two examples that will be considered in detail: Fig. 1a shows a thin rod (an example of a system that obeys the wave equation) undergoing axial vibrations with displacement u(x,t), while Fig. 1b shows an Euler-Bernoulli beam, which vibrates transversely in bending with displacement w(x,t). The boundary applies a nonlinear axial or transverse force on the waveguide (the beam is assumed to have zero moment at the end for simplicity), which depends only on the displacement at x = 0. It is assumed that the stiffness and the force can be written as a polynomial in the spring extension.
We now assume that a group of time harmonic waves of amplitudes a þ n at frequencies nx (n integer) are incident upon the boundary at x = 0. The source of the incident waves might be either excitation or the reflection of the negative-going waves from a boundary or discontinuity at some point x\0. The incident waves produce a group of reflected waves of amplitudes a À n and, for the beam, reflected nearfield waves of amplitudes a À N:n at the same frequencies. The response in the waveguide is assumed to be the superposition of the positive-and negative-going waves.
The case of flexural vibrations of an Euler-Bernoulli beam with a nonlinear boundary is now considered. The equations for the simpler and less interesting system of a thin rod in axial vibrations can be found using the same approach and are merely stated at the end of this section.

Flexural vibrations of an Euler-Bernoulli beam
The governing equation of motion for flexural vibrations of the beam is [1]: where w x; t ð Þ is the transverse displacement and q, S, E, I are the density, cross-sectional area, elastic modulus and second moment of area of the crosssection. For time-harmonic behaviour at frequency x free waves of the form exp i xt Ç jx ð Þcan propagate in the positive and negative x-directions, respectively, together with nearfield or evanescent (non-propagating) waves of the form exp ixt Ç jx ð Þ , whose Here is the wavenumber, while the phase velocity c p ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi EI=qS is a function of frequency, the flexural waves being dispersive. For the wave groups in Fig. 1b, the wavenumber of the n th harmonic at frequency nx is thus ffiffi ffi n p j for flexural waves. The spring applies a force, which is assumed to be a polynomial function of the displacement at x = 0 with constant coefficients k m so that the boundary conditions, assuming zero bending moment M, can be written as where V is the shear force. The displacement of the beam is assumed to be the superposition of the positive-and negative-going waves shown in Fig. 1b and thus can be written as: where CC denotes the complex conjugate of the preceding terms. The ratios of the amplitude of the n th incident wave, reflected propagating and reflected nearfield wave relative to that of the first incident wave are defined as: Note that the reflection coefficients r n and r N;n are functions of a þ 1 and a n , in contrast to the linear case. Substituting (5) in (4), the flexural displacement of the beam in terms of a n ,r n and r N;n becomes: a n e i nxtÀ ffiffi r n e i nxtþ ffiffi 2.1.1 Effect of the phase of the fundamental incident wave The wave amplitudes are complex and can be expressed in terms of magnitude and phase. The phases of the n th incident, reflected and reflected nearfield waves are defined as U n ; / n and w n , respectively. First, assume a þ 1 is pure real, i.e.
Þand reflected nearfield waves a À N;n exp i nxt þ w n ð Þ ð Þ . It follows that Now suppose that a þ 1 is complex, i.e. U 1 6 ¼ 0, so that the first incident wave at Defining t 0 ¼ t þ U 1 =x, the relations between the phases of the reflected waves are given by: in which U 0 n ¼ U n À nU 1 . For simplicity, henceforth it is assumed that U 1 ¼ 0 unless otherwise stated.

Cubic nonlinear stiffness
Generally, the reflection coefficients can be found by substituting Eq. (6) for the response into the boundary conditions (3), separating terms at each frequency and solving. Given that the nonlinear boundary spring stiffness can potentially include an infinite number of terms (i.e. M ¼ 1 in Eq. 3), this can be a formidable task. However, for simplicity, it is now assumed that the spring possesses cubic nonlinearity and that k 2 ¼ 0. Under this assumption, an incident wave at one frequency produces reflected waves with only odd multipliers of the fundamental frequency [17,18]. Consequently, the second part of Eq. (3) becomes.
The bending moment is zero and it follows that This is substituted into Eq. (6) Putting x = 0 in Eq. (11) results in The shear force at x = 0 is given by: To put the equations in a non-dimensional form, a reference frequency is introduced. For a rectangular cross-section of width b and thickness h, then jh ¼ 1 at this frequency. For a beam with an arbitrary cross-section, ffiffiffiffiffi 12 p jr g ¼ 1 at this frequency, where r g is the radius of gyration of the cross-section. The characteristic impedance (the ratio of the shear force to transverse velocity for a propagating wave) of the beam at x 0 is are introduced, where X represents the non-dimensional frequency ratio, K is the relative linear stiffness of the boundary spring (k w 0 ð Þ) to that of the beam at frequency x 0 and represents linear effects, while A is the ratio of the nonlinear force in the spring (due only to the first incident wave) and that of the beam at frequency x 0 and represents the nonlinear effects.
Equation (14) can now be rewritten in terms of the dimensionless parameters as X 3=2 X 1 n¼1;3;5;::: n¼1;3;5;::: n¼1;3;5;::: It can be seen that the reflection coefficients depend on the three dimensionless parameters together with the magnitudes and phases of a n . In this paper, weak nonlinearity is defined to be the case when the dimensionless nonlinear force is small compared to the linear forces of the spring and beam together, i.e.
. This is in contrast to previous studies [17,18] and [20,21] where weak nonlinearity was defined to be the case when the relative forces of the linear and nonlinear boundary springs are small, i.e. A ( K or k 3 a þ2 Our definition allows us to study the case of essential nonlinearity (K = 0 or k 1 = 0) at which the nonlinearity can be stronger. However, it should be noted that K = 0 does not necessarily mean strong nonlinearity since the linear dynamic stiffness of the beam can be large compared to A, i.e. A . X 3=2 may be small and the nonlinearity is weak.

Power and conservation of energy
The time average power in a time-harmonic propagating wave of amplitude a AE at frequency x is is the group velocity of flexural waves. Due to energy conservation at the boundary, the net power, i.e. the difference between the powers of the positive-and negative-going waves, must be zero, the nonlinearity causing energy to leak into higher harmonics. However, evanescent (nearfield) flexural waves carry no power. Hence, the net time average power (P) of the incident and reflected flexural waves is given by: When there is only one incident wave, Eq. (17) in terms of the reflection coefficients becomes:

Axial waves in a thin rod
The analysis for the case of a rod in axial vibrations follows the same steps as above, although there are of course no nearfields. Only the main results are stated here. Axial vibrations of a rod are governed by equation where u is the axial displacement and P is the axial force. The boundary condition at x = 0 for cubic stiffness is For time-harmonic motion, the wave speed c g , wavenumber j and characteristic impedance Z of the rod at frequency x are [1] Axial waves are non-dispersive so that the wavenumber of the n th harmonic wave at frequency nx is nj.
The displacement of the rod can be written as: Substituting (5) into equation (22) gives The characteristic impedance Z is given by Eq. (21), and the dimensionless parameters of equation (15)  n a n À r n ð Þe inxt n¼1;3;5;::: n¼1;3;5;::: Weak nonlinearity is again defined to be the case when A=ðK þ XÞ ( 1. As can be seen, Eq. (24) is very similar to Eq. (16); hence, in Sects. 3, 4 and 5 equations for flexural waves are analysed in detail, and equivalent equations for axial waves are presented in Appendix B. Note that the time average power of an axial wave at frequency x and of amplitude a þ equals 1 2 Zx 2 a þ j j 2 .
Hence, the net time average power or, for the case of one incident wave in terms of reflection coefficients, 3 The case of only one incident wave and higher harmonic reflected waves Consider the case when there is only one incident wave, i.e. a n ¼ 0; ðn [ 1Þ, so that Eq. (16) becomes n¼3;5;::: n¼1;3;5;::: n¼1;3;5;:::

Truncation and solution
The infinite set of equations (27) can be truncated at some finite value of n = N and solved through the harmonic balance method. The accuracy of the solution depends on the value chosen for N and is discussed in Sect. 6.2. The terms of e inxt are collected and the coefficients are equated to zero. For example, taking N = 9 and collecting the terms in e ixt , the equation for the first harmonic is: where the overbar represents the complex conjugate. Similarly, collecting and equating the coefficients of time-dependent terms at frequency 3x result in In a similar way, by collecting terms of higher harmonics, equations for n = 5, 7 and 9 can be found as presented in Appendix A. The analysis for large values of N can become extremely cumbersome.

Case of one reflected wave (N = 1)
In previous studies [20,21], it was shown that if all higher harmonics are neglected, i.e. truncating the system of equations at N = 1, the nonlinearity has no effect on the magnitude of the reflection coefficient r 1 , only a change in its phase. Thus, the magnitude of the reflection coefficient is similar to that of the linear case, i.e. r 1 j j ¼ 1. Since there is no energy dissipation at the boundary, all the incident energy is necessarily reflected into the first harmonic. It should be noted that putting N = 1 in Eq. (16) results in the same equations as those developed in [20,21].

A solution for the case of one incident wave and two reflected waves
In order to give insight into the underlying behaviour of the system, first the case when there is only one incident wave and only the first and third reflected harmonics are retained is discussed. Truncating Eq. (27) at N = 3 results in This set of equations can be solved for various values of parameters. Special cases are when A or K is zero. When A = 0, then the system is linear with In this case, the denominator is the complex conjugate of the numerator so that |r 1 |= 1, and since energy is conserved, r 3 = 0. At low-frequency ratios ðX 3=2 . K ( 1Þ, the behaviour of the beam is similar to one with a simply supported boundary and r 1 ¼ r SS ¼ À1, while, at high frequencies (X 3=2 . K ) 1), the beam behaves the same as one with a free end and r 1 ¼ r free ¼ Ài.

Essential nonlinearity
When the linear stiffness K = 0, the wave reflection is governed by essential nonlinearity. Note that even in this case the nonlinearity can be weak when A\X 3=2 . A new parameter is now introduced. Equations (30 and 31) can thus be rewritten in terms of c as It can be seen that the reflection coefficients depend only on c. The case c ( 1 denotes weak nonlinearity and the beam behaviour effectively corresponds to that of a beam with a free end, i.e. r 1 ¼ r free ¼ Ài; r 3 ¼ 0. On the other hand, if c ) 1 the beam behaves similar to a beam with a simply supported boundary and r 1 ¼ r SS ¼ À1; r 3 ¼ 0. Hence, at some value of c, |r 1 | reaches a minimum and |r 3 | reaches a maximum, which implies that the effects of nonlinearity are strongest.

Analytical solution using small parameter expansion for weak nonlinearity
The method of expansion in a small parameter is used to find an approximate solution to Eqs. ( (30 and 31). Assuming e to be a small parameter, eA denotes weak nonlinearity. Hence, Eqs. (30 and 31) become The solution to Eqs. (36 and 37) in terms of the small parameter e is r 1 ¼ r 01 þ er 11 þ e 2 r 21 þ Oðe 3 Þ; ð38Þ Since the difference between the amplitudes of the first incident and reflected waves is at least of order Oðe 2 Þ, we consider a second-order expansion. Substituting (38 and 39) in (36 and 37) and equating the coefficients of e n , the terms in the expansion of r 1 are: Similarly, for r 3 r 03 ¼0; As can be seen from the coefficients of e 0 , r 03 ¼ 0 and r 01 ¼ r lin and the reflection coefficients depend on the linear parameters of the system (K and X). A appears in the expressions for the higher-order terms of the series, so that r 11 and r 21 are functions of can be noted that r 01 is of order 1, The perturbation method is accurate when the nonlinearity is weak, i.e. b 1 j j is small.

Case of two incident and reflected waves
It is now assumed that in Eq. (16) a 3 6 ¼ 0, i.e. the second incident wave with frequency 3x is considered. Taking N = 3 in Eq. (16) and collecting the first and third harmonic terms result in: For the case of essential nonlinearity, these equations can be written in terms of c as  The case c ( 1 denotes weak nonlinearity and the beam behaviour effectively corresponds to that of a beam with a free end, i.e. r 1 = -i. We are interested in investigating the effects of stronger nonlinearity on the reflection coefficients with increasing c. It can be seen from the figures that, first, when nonlinearity increases |r 1 | decreases. However, when A or c is very large, the beam behaves similar to a beam with a simply supported boundary, as expected, and r 1 = -1. Since energy is conserved, the magnitude of the reflection coefficient of the first harmonic for the case of one incident wave cannot be larger than one and energy leaks from the first harmonic to the third harmonic. Consequently, there is a minimum value of |r 1 | when the effects of nonlinearity are strongest. The leakage of energy depends on the value of c. With increasing c, first |r 1 | decreases while |r 3 | increases. At c ' 0.073, the minimum magnitude of r 1 is about 0.8, and the maximum magnitude of r 3 is about 0.15. From this point, as c increases, the end of the beam will be effectively pinned and r 1 asymptotes to -1. Figures 2c and d show the changes in the phases of reflection coefficients. When c\\1, the phases of the first and third harmonics are -90°. As c increases, / 1 and / 3 asymptote to 180°and 90°, respectively.  Figure 3 shows the magnitudes of the reflection coefficients of the first and third harmonics as a function of frequency ratio X for A = 0.01 and three different values of K. The minimum magnitude of the reflection coefficient of the first harmonic and the maximum magnitude of the reflection coefficient of the third harmonics are defined as |r 1 | min and |r 3 | max . As can be seen in Fig. 3, when K = 0, |r 1 | min is minimum and |r 3 | max is maximum. As K increases, |r 1 | min increases since the linear stiffness becomes more pronounced and decreases the effect of nonlinearity. From energy conservation r 1 j j 2 min þ9 ffiffi ffi 3 p r 3 j j 2 max ¼ 1, and hence, in Fig. 3, as |r 1 | min increases, |r 3 | max necessarily decreases.
In Fig. 4, the effect of the nonlinear stiffness on the magnitudes of the reflection coefficients is examined. Results are presented for K = 0.1 in Fig. 4a and b as a function of the frequency ratio X. As can be noted from Fig. 4a, increasing A changes the frequency at which |r 1 | min and |r 3 | max occur. Furthermore, as A increases |r 1 | min decreases and asymptotes to about 0.8 and |r 3 | max increases and asymptotes to about 0.15.
Examining the results in Figs. 3 and 4 indicates that generally, as the nonlinear effects become more pronounced, i.e. the linear component (K) decreases or nonlinear component (A) increases, |r 1 | min and |r 3 | max asymptote to those values in the case of essential nonlinearity in Fig. 2. Results are shown only for the case of essential nonlinearity in the rest of the paper since this case is of primary interest and effect of nonlinearity is more pronounced.

Flexural waves: one incident wave, multiple reflected waves and convergence
In this section, the solutions to the system of equations truncated at a higher number of harmonics (N = 5, N = 7, and N = 9) are presented and compared to the results with N = 3 for the case of essential nonlinearity. Results show that the underlying behaviour of the system is the same as that seen when retaining only the first and third reflected waves and that the solutions converge as N increases. The set of equations are truncated for various N, and the magnitudes and phases of the reflection coefficients are illustrated in Fig. 5 as a function of c. The qualitative behaviour is the same as when N = 3. There is a significant difference between |r 1 | min and |r 3 | max for N = 3 and 5; however, these two parameters are almost identical for N = 5, 7 and 9. Comparing the magnitudes of r 1 to r 9 shows that, as the order of the odd harmonic increases, the magnitude of the reflection coefficient decreases, i.e. r 1 j j[ r 3 j j[ r 5 j j[ r 7 j j[ r 9 j j. Retaining a higher number of harmonics has less effect on the phases of the reflected waves than on their magnitudes.
To study the accuracy of the solution, the relative error r n j j m N of the maximum or minimum magnitudes of the n th reflected wave retaining N harmonics is defined as the error with respect to r n j j m 9 . Table 1 shows the errors together with |r n | m when N = 9. Truncating the equations at N = 9 results in very good accuracy for the reflection coefficients of lower harmonics (r 1 and r 3 ), good accuracy for r 5 and a reasonable accuracy for r 7 .
The parameter P n , which denotes the normalized power n 2 ffiffi ffi n p r n j j 2 of the n th harmonic, is shown in Fig. 6 as a function of c with essential nonlinearity (K = 0). When the system is linear P 1 ¼ 1, however, as illustrated in Fig. 6, with nonlinearity the power leaks from the first harmonic to the higher harmonics. Comparing the powers shows that P 1 [ P 3 [ P 5 [ P 7 [ P 9 . As expected, the incident energy is mostly reflected into the first and third harmonics and only about 5% reflected to the higher harmonics.
Furthermore, there is substantial leakage of energy (' 33%) from the first harmonic to the third harmonic.
From the net power equation for flexural waves (Eq. (18)), the accuracy of the harmonic balance solution can be defined as the difference between the total power of the reflected waves and unity (the power of the incident wave). For example, for flexural waves with N = 9, the error is With small enough solver tolerances that used in MATLAB, this error is less than 10 -8 .

Flexural waves: case of two incident and reflected waves
The results presented in this section concern the effect of multiple incident waves on the wave reflections at the end of the beam for the case of essential nonlinearity. The reflection coefficients now depend on c together with the magnitude and phase of a 3 . The source of the multiple incident waves might be either multi-harmonic excitation or the reflection of the negative-going waves with different frequencies and phases from a boundary or discontinuity at some point x\0.
6.3.1 Flexural waves: the effect of U 3 on the magnitudes of the reflection coefficients Figures 7 and 8 show the effect of the phase difference U 3 between the incident waves with a 3 j j ¼ 1 as a function of c. (Fig. 7 is a contour plot, while Fig. 8 is a 2D plot for the same parameters). It can be seen that the reflection coefficients depend strongly on the phase of the incident third harmonic. For example, when U 3 ¼ p and c ( 1, r 1 j j ' 1, then, with increasing c, |r 1 | increases and reaches a local maximum and then decreases and reaches a minimum, again increases and reaches a local maximum and finally asymptotes to 1 (|r 1 | has two maximums and one minimum). However, when U 3 ¼ 0 and c ( 1, r 1 j j ' 1, then as c increases the magnitude of r 1 decreases and reaches a minimum, then increases and reaches a maximum, again decreases and reaches a minimum and for c ) 1 asymptotes to 1 (|r 1 | has one maximum and two minimums). As illustrated in Figs. 7a and 8a, the minimum value |r 1 | min occurs when U 3 ' p, while |r 1 | max occurs at U 3 ' 0. Since energy is conserved r 1 j j 2 þ9 ffiffi ffi 3 p r 3 j j 2 ¼ 1 þ 9 ffiffi ffi 3 p a 3 j j 2 , and hence when r 1 j j is maximum r 3 j j is minimum as shown in Figs. 7b and 8b. Figures 7a and 8a show that, in contrast to the case of one incident wave, r 1 j j can be significantly larger than 1 if U 3 ' 0. Furthermore, with two incident waves, |r 1 | min can be smaller than when a 3 ¼ 0. When U 3 ' 0, both incident waves are in phase when t ¼ np=x and produce a larger displacement. On the other hand, when U 3 ' p the incident waves are out of phase and produce a smaller displacement. As illustrated in Figs. 7b and 8b, although U 3 has a significant effect on r 1 j j, it does not substantially affect r 3 j j due to dominance of the direct reflection of the incident third harmonic. Figures 9 and 10 show the magnitudes of the reflection coefficients as functions of c and a 3 j j for (i) (j) Fig. 9 is a contour plot, while Fig. 10 is a 2D plot for the same parameters). The figures illustrate the effect of a 3 j j on the magnitudes of the reflection coefficients and identify the value of a 3 j j at which |r 1 | min is minimum (when U 3 ' p) and |r 1 | max is maximum (when U 3 ' 0). It can be seen from Fig. 9a or 10a that when both incident waves are pure real, |r 1 | can be greater than 1 and reaches a maximum of about 2.5 at a 3 j j ' 1:5. Furthermore, in Fig. 9c or 10c, |r 1 | reaches a minimum of about 0.4 at a 3 j j ' 0:5, which is significantly smaller than the case with one incident wave ( r 1 j j min ' 0:8). Figures 9b or 10b and 9d or 10d show that the reflected third harmonic, a À 3 , is dominated by the incident wave, a þ 3 .

Flexural waves: multiple solutions
For some parameter values, there are multiple solutions for the reflection coefficients. This behaviour is not uncommon in nonlinear systems, e.g. the Duffing oscillator. These are accompanied by jumps in the values of |r| as c increases or decreases, i.e. by exciting the beam at certain frequencies there are two possible reflection coefficients for the beam. When c ( 1, there are no jumps in the magnitudes of the reflection coefficients; however, as c increases, for some values of U 3 and a 3 j j the effects of nonlinearity are strong and the magnitudes of r 1 and r 3 are multivalued, undergoing a rapid change or jump as c changes. When both incident waves are pure real, there is no jump, as shown in Figs. 9a, b, 10a and b. In Figs. 7 and 8, multiple solutions occur when 0:0039\c\0:0071 and p=3 U 3 5p=6, while in Figs. 9c, d, 10c and d, there are two solutions for the reflection coefficients in the range 0:0126\c\0:0215 and 0:5 a 3 j j 0:8. As demonstrated in Figs. 11 and 12, in these parameter ranges, the magnitudes of the reflection coefficients become multi-valued for some values of c and there are at least two possible solutions for |r 1 | and |r 3 |. Results for the reflection coefficients are shown for both increasing and decreasing c. As can be seen in Figs. 11 and 12, due to energy conservation, when |r 1 | jumps up from the lower branch to the upper branch, |r 3 | jumps down. Table 1 Flexural waves: reflection coefficients, reflected power and the percentage error of the minimum or maximum magnitude of the reflection coefficients with respect to the value when N = 9

Flexural waves
The minimum or maximum magnitude of the reflection coefficients for each harmonic r n j j m9 n 2 ffiffi ffi n p r n j j m9 2 Percentage relative error  For example, in Fig. 12, when a 3 j j ¼ 0:6, the magnitudes of r 1 and r 3 change substantially in the region 0:0175\c\0:018. The difference between the two solutions is largest at U 3 ¼ 2p=3 in Fig. 11 and at a 3 j j ¼ 0:6, in Fig. 12.
Note that to have an accurate solution, the tolerances should be small enough. The error of the numerical method for multiple combinations of initial conditions is less than 10 -8 .
Numerical studies have been performed to determine the magnitude and phase of a 3 for which there are multiple solutions. The parameter region in which the solution is multi-valued is shown in Fig. 13. When 40 0 \U 3 \193 0 or 0:45\ a 3 j j\1:27, there is more than one solution for the magnitudes of the reflection coefficients for some values of c.

Axial waves
In this section, various numerical examples are presented for axial waves with essential nonlinearity. First, results with only one incident wave and two reflected waves being retained (i.e. N = 3) are presented. Then, the convergence of the solution is illustrated. Finally, the case with two incident and reflected waves is considered.

Axial waves: the case of one incident wave and two reflected waves
The magnitude and phase of the reflection coefficients of the first and third harmonics are shown in Fig. 14 as a function of c ¼ ðA=XÞ with K = 0. Numerical results are shown, together with the analytical solutions for weak nonlinearity presented in Appendix B. As for the case of flexural waves, results are in good agreement for small values of c. However, as c increases, the nonlinear effects become more pronounced, and the analytical results become less accurate: more terms need to be included in the expansion to give better accuracy, but this becomes increasingly intractable, especially when more harmonics are retained. As illustrated in Fig. 14, the qualitative behaviour of the rod is similar to that of the beam. When c ( 1 the end of the rod is in effect free and r 1 ¼ 1. The magnitude of the reflection coefficient of the first harmonic has a minimum value of about 0.9193 at c ' 0:29, while the magnitude of r 3 has a maximum value of 0.1311. At this point, the effects of the nonlinearity are strongest and the maximum leakage of energy from the first to the third harmonic occurs. As c increases further, r 1 increases while r 3 decreases. When c [ [ 1, the end of the rod effectively becomes fixed and r 1 ¼ À1. In Fig. 14c, it can be seen that the Fig. 9 Flexural waves: Retaining two incident and two reflected waves: the magnitudes of the reflection coefficients of the first and third harmonics for the case of essential nonlinearity as functions of c and a 3 j j: a |r 1 |, 10 Flexural waves: Retaining two incident and two reflected waves: the magnitudes of the reflection coefficients of the first and third harmonics for the case of essential nonlinearity with various a 3 j j as a function of c: a |r 1 |, 6.4.2 Axial waves: higher harmonics for one incident wave Table 2 shows the results of a convergence study similar to that in Sect. 6 When only the first and third harmonics are retained, |r 1 | min = 0.9193. When the fifth harmonic is retained |r 1 | min decreases to 0.9141; however, it does not change significantly when higher harmonics are retained. When the fifth harmonic is retained, |r 3 | max decreases noticeably from 0.1312 (when only the first and third harmonics are retained) to 0.1267. When the seventh harmonic is taken into account, |r 3 | max decreases slightly to 0.1264. Furthermore, retaining the ninth harmonic has little effect on the magnitude of r 3 . As for the case of the beam, as higher harmonics are retained, the maximum magnitudes of the higher reflection coefficients decrease and r 3 j j max [ r 5 j j max [ r 7 j j max [ r 9 j j max . The normalized power in the n th reflected harmonic for axial waves, n 2 r n j j 2 is also shown in Table 2. As can be seen, n 2 r n j j 2 decreases rapidly with increasing n. Retaining N harmonics gives moderate accuracy for r N , good accuracy for r N-2 and very good accuracy for the reflection coefficients for lower harmonics.
Similar to the case of flexural waves, a 3 j j affects the reflection coefficients substantially, but the behaviour is a little different from those of flexural waves. This is due to the interaction of waves with different frequencies, which makes the vibration of the system more complex. When U 3 ' p, r 1 j j is minimum and hence r 3 j j is maximum. Furthermore, comparing Figs. 14 and 17 it is seen that the minimum magnitude of r 1 can be smaller than when there is just one incident wave. Additionally, the magnitude of r 1 can be larger than one. As can be seen in Figs. 15a or 16a, the maximum value of r 1 j j is about 1.4 and occurs when U 3 ' 0. In contrast to the first harmonic, as illustrated in Fig. 15b or 16b, U 3 has little effect on r 3 due to the dominance of the third incident wave. When a 3 j j ' 0:7 and U 3 ¼ 0, r 1 j j is maximum and r 3 j j is minimum; however, when a 3 j j ' 0:4 and U 3 ¼ p, r 1 j j is minimum and r 3 j j is maximum. Examining Figs. 17 and 18, it can be noted that when a 3 j j[ 1, a þ 3 is dominant and, r 3 j j ' a 3 j j.

Conclusion
This paper presented analytical and numerical studies of the reflection coefficients for time-harmonic waves in a waveguide from a boundary with nonlinear stiffness. A wave solution was applied to a semiinfinite, thin rod and a semi-infinite Euler-Bernoulli beam. Incident, reflected propagating and nearfield waves comprising a series of harmonics at integer multiples of the frequency of the first incident wave were included. For simplicity, the nonlinearity was assumed to be cubic, with only odd harmonics incident and reflected at the boundary. An infinite set of equations was derived, and solutions were developed using the harmonic balance method by truncating these and including a finite number of harmonics. The conclusions obtained are as follows.
(1) Dimensionless parameters It was shown that the reflection coefficients depend on three dimensionless parameters: the frequency ratio, the linear stiffness and the nonlinear stiffness of the spring, together with the relative amplitudes of the incident waves a n . (2) Strong nonlinearity Weak nonlinearity was defined to be the case when A= K þ X ð Þ(1 for axial waves and A . K þ X 3=2 ( 1 for Fig. 13 Flexural waves: the region in which the solution is multi-valued Table 2 Axial waves: Reflection coefficients, reflected power and the percentage error of the minimum or maximum magnitude of the reflection coefficients with respect to the value when N = 9

Axial Waves
The minimum or maximum magnitude of the reflection coefficients for each harmonic r n j j m9 n 2 r n j j m9 2 Percentage relative error  flexural waves. So that the case of essential nonlinearity was studied. For essential nonlinearity (K = 0), the reflection coefficients depend only on the nondimensional parameter c ¼ A=X and a n . (3) Minimum reflection with one incident wave For the case with one incident wave and when two reflected waves are retained (N = 3), the method of expansion in a small parameter was used to find analytical expressions for the reflection coefficients of the first and third harmonics.
Results were compared with numerical solutions for essential nonlinearity. It was seen that the magnitude of the reflection coefficient of the first harmonic for axial and flexural waves has a minimum and energy leaks from the first harmonic to the third harmonic. Since energy is conserved, with increasing c, the magnitude of the reflection coefficient of the first harmonic decreases and reaches a minimum, while the (a) (b) Fig. 15 Axial waves: Retaining two incident and two reflected waves; the magnitudes of the reflection coefficients of the first and third harmonics for the case of essential nonlinearity with a 3 j j ¼ 1 as functions of c and U 3 a |r 1 |, b |r 3 | Fig. 16 Axial waves: Retaining two incident and two reflected waves; the magnitudes of the reflection coefficients of the first and third harmonics for the case of essential nonlinearity with a 3 j j ¼ 1 as functions of c and U 3 a |r 1 |, b |r 3 | magnitude of r 3 increases and reaches a maximum. The minimum magnitude of r 1 was about 0.92 for axial waves and about 0.8 for flexural waves. For a single incident wave, the maximum energy that leaks from the first harmonic wave is about 16% of the incident energy for axial waves and 36% for flexural waves. The qualitative behaviour was similar for axial and flexural vibrations. (4) Linear and nonlinear effects The effects of K and A on |r 1 | min and |r 3 | max were studied. It was shown that as we increase nonlinearity or decrease linearity of the system, the minimum or maximum magnitudes of the reflection coefficients asymptote to those values for essential nonlinearity, i.e. for large values of A/K, r 1 j j min ' 0:8 for flexural waves and r 1 j j min ' 0:92 for axial waves. (5) Higher harmonics and convergence For the case with one incident wave and truncating the system of equations at N = 5, 7 and 9, the reflection coefficients and powers in each harmonic were obtained and compared with the results for N = 3. The results showed that the qualitative behaviour of the system with a larger number of retained reflected waves is the same as that with N = 3 and that the solution was seen to converge for increasing N. Retaining the fifth harmonic gives very good accuracy for the first and third harmonics. (6) Multiple incident waves The case of two incident and two reflected waves was studied. The Fig. 17 Axial waves: The magnitudes of the reflection coefficients of the first and third harmonics as functions of c and a 3 j j, retaining only the first and third harmonics: a |r 1 |, U 3 ¼ 0, b |r 3 |, U 3 ¼ 0, c |r 1 |, U 3 ¼ p, d |r 3 |, U 3 ¼ p Fig. 18 Axial waves: The magnitudes of the reflection coefficients of the first and third harmonics for various values of a 3 j j as a function of c and retaining only the first and third harmonics: a |r 1 |, U 3 ¼ 0, b |r 3 |, U 3 ¼ 0, c |r 1 |, U 3 ¼ p, d |r 3 |, U 3 ¼ p reflection coefficients now depend on the magnitude and phase of the ratio of the amplitudes of the third to first harmonic incident waves a 3 . It was seen that the reflection coefficients depend significantly on the incident third harmonic, as described in Sects. 6.3 and 6.4.3. There are qualitative differences in the dependencies of the reflection coefficients on the dimensionless frequency for the rod and beam. For example, for the beam, there can be two local maximums and one minimum (or one maximum and two local minimums), while for the rod only one maximum and one minimum were observed. A numerical study to show the minimum and maximum reflection coefficients of the first harmonic with two incident waves was done. When both incident waves are pure real, r 1 j j has a maximum value significantly larger than 1, while when the phase difference of incident waves is p, r 1 j j min is smaller than the case with only one incident wave. (7) Multiple solutions for reflection coefficients It was seen that for certain values of a 3 j j and U 3 there can be two solutions for the reflection coefficients of flexural waves: jump phenomena occur as c increases or decreases. Because energy is conserved, when |r 1 | jumps up from the lower branch to the upper branch, |r 3 | jumps down. Numerical studies were performed and the region in which the system becomes multivalued was determined. No jump phenomenon was found for the reflection of axial waves.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions. This article was supported by University of Auckland.
Data availability Enquiries about data availability should be directed to the authors.

Declarations
Conflict of interest The authors have not disclosed any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.