Wave-based analysis of jointed elastic bars: stability of nonlinear solutions

In this paper we develop two new approaches for directly assessing stability of nonlinear wave-based solutions, with application to jointed elastic bars. In the first stability approach, we strain a stiffness parameter and construct analytical stability boundaries using a wave-based method. Not only does this accurately determine stability of the periodic solutions found in the example case of two bars connected by a nonlinear joint, but it directly governs the response and stability of parametrically forced continuous systems without resorting to discretization, a new development in of itself. In the second stability approach, we pose a perturbation eigenproblem residue (PER) and show that changes in the sign of the PER locate critical points where stability changes from stable to unstable, and vice-versa. Lastly, we discuss follow-on research using the developed stability approaches. In particular, we identify an opportunity to study stability around internal resonance, and then identify a need to further develop and interpret the PER approach to directly predict stability.


Introduction
In a companion paper [1], we developed two wavebased methods for predicting the nonlinear vibration response of jointed, continuous elastic bars. These methods consist of a nonlinear wave-based vibration approach (WBVA) informed by re-usable, perturbationderived scattering functions, and a numerical plane wave expansion (PWE) approach exploiting wave solutions as expansion quantities. Both approaches were applied to finding periodic solutions in example systems in which continuous bars, connected by nonlinear joints, are constrained at their ends and harmonically excited. As documented in [1], the two approaches exhibit very good accuracy while significantly reducing the computation time required to obtain periodic solutions. They are also notable for having fixed problem size, independent of frequency. Arguably, stability of the periodic solutions predicted is as important as the solutions themselves. In [1] we report stability results for the examples studied using two new stability methods, which we deferred discussion of until the present paper.
Direct methods for assessing stability of wave-based solutions have not been presented in the literature, the necessity of which has been stated in recent studies (see, for instance [2,3]). When stability has been treated, previous studies, owing to a lack of alternatives, have chosen to project the wave-based solution onto a finite element mesh and then use time-domain monodromy matrix calculations to estimate stability. Although this approach is very attractive owing to its simplicity, it leads to inaccuracies at larger amplitudes, which we document in Appendix A. In this paper, we assess stability of the multi-harmonic periodic solutions found in [1] using two new approaches: a strained parameter, analytical method applied to the governing partial differential equations based loosely on the technique as applied to parametrically excited ordinary differential equations [4], and a computational method employing the residue of the perturbed system's eigenvalue problem. Both are informed by Floquet theory. Floquet theory (see [4]) provides the theoretical basis for the fact that stability transitions of a system perturbed around a periodic solution occur when the perturbed system has a purely imaginary eigenvalue with magnitude equaling an integer multiple of the frequency of the response (among other cases, see Sect. 3 below). Detection of such points along a forced response curve, as opposed to a complete spectral decomposition of the perturbed system for each response on the curve (such as the HBM perturbation approach in [5]), underpins both proposed approaches. In the strained parameter approach we find stability tongues in the parameter space distinguishing stable and unstable solutions, while in the perturbed eigenproblem residue approach we employ Floquet theory to find changes in the computed residues which indicate a stability transition.
The paper is organized as follows: we first develop an analytical, wave-based stability approach based on straining a single stiffness parameter in Sect. 2. This strained parameter approach yields stability curves (i.e., Arnold tongues) and surfaces associated with a parametrically excited continuous system. Next, in Sect. 3 we develop a computational approach termed the perturbation eigenproblem residue (PER) for predicting stability changes in the forced problem. In Sect. 4 we apply both methods to determine the stability of the two example systems studied in [1]. The paper concludes with discussions of the results and avenues for future work.

Stability via the strained parameter approach
We now formulate an analytical approach to determine stability of the multiple periodic solutions found in [1] for the single-jointed system shown in Fig. 1. The approach shares many similarities with the "strained parameter", or Lindstedt-Poincaré, approach employed in the study of the damped Mathieu equation [4]. While the Mathieu equation is a single ordinary differential equation, the strained parameter approach developed herein assesses stability of continuous systems governed by partial differential equations. In addition, the motivating example studied consists of multiple continuous domains connected by a nonlinear joint.
The equations governing vibration response in the linearly elastic bars are given by, where the left and right bars have displacement fields denoted by u L (X, t) and u R (X, t), while f L (X, t) and f R (X, t) denote forcing per unit length in the left and right subdomains, respectively. For simplicity of resulting expressions, the two bars considered share the same Young's modulus E y , density per unit volume ρ, and cross-sectional area A. The three values for α i , i = 1, 2, 3, form a partition of unity. Four boundary conditions are required to completely specify the vibration response of two jointed bars. The two clamped ends provide two boundary conditions, u L (0, t) = 0 and u R (l, t) = 0, while two others arise at the joint, where X J = (α 1 + α 2 )l. The joint considered has linear and cubic restoring stiffness characterized by coef- Fig. 1 Two linearly elastic bars joined by a nonlinear spring and a linear dashpot. Each bar has one fixed end. Singlefrequency plane waves are injected at the location of the forcing. Note that each wave coefficient includes itself plus its complex conjugate-this is only shown for a + , but assumed for all other coefficients. The global coordinate X denotes position starting from the left end. Periodic solutions to this problem have been found in [1] ficients K and ε , respectively, and linear damping characterized by C. Note that ε denotes a small bookkeeping parameter later set to one [4]. We identify any periodic solution [1] with left and right displacement fields denoted by u * L (X, t) and u * R (X, t), respectively. Next, we assess the local stability of these solutions by introducing small perturbations δu L (X, t) and δu R (X, t) such that the total displacement fields are now given by, Substituting the above into the governing field equations, Eqs. (1), and (2), the joint conditions, Eqs. (3), and (4), and the zero displacement boundary conditions, we retain only linear terms in δu L (X, t) and δu R (X, t), yielding where u * denotes the joint displacement. Retaining only terms up to and including O(ε) in the stability analysis (i.e., only retaining the fundamental harmonic in u * J (t)), we rewrite the bracketed term on the right-hand side of Eqs. (9), and (10) as, where we have used u * J (t) ≡ a * cos (ωt + φ * ), φ * denotes the joint displacement phase relative to the forcing, K ≡ K + 3ε a * 2 2 , and P ≡ 3 a * 2 2 . We note that both the amplitude, a * , and phase, φ * , are known from the solutions obtained in [1], and these fully distinguish one solution from another.
Using Eq. (13) in Eqs. (7-12), we recover a parametrically excited problem as illustrated in Fig. 2(a). While this problem arises from studying the stability of a directly excited system, it may be of general interest. For example, similar problems will arise when two domains are connected by time-varying stiffness, such as in parametrically driven micromechanical oscillators [6]. For the forced problem, of interest herein, stability of the solution is guaranteed when the parametrically excited problem is stable. Fig. 2 a Damped, parametrically excited system arising from the local stability analysis of the corresponding forced system; solution approach at b zeroth-order, c first-order, and d second-order We next proceed to find wave-based solutions to the linear problem illustrated in Fig. 2(a). To do so, we employ a strained parameter approach. We first expand δu L and δu R , and subsequently expand the static component of the spring stiffness, which amounts to "straining" the stiffness in the problem. We also place the damping at first order to avoid decay in the ensuing zeroth-order problem: C → Ĉ . Substitution of the expanded and ordered quantities into Eqs. (7-13) results in three ordered problems, as illustrated in Fig. 2(b-d), governed by, where the forcing arising at the joints at O(ε 1 ) and O(ε 2 ), respectively, is given by: Similar to the wave-based solution procedure for the nonlinear forced problem [1], we seek exact wave solutions for the three linear problems shown in Fig. 2 Starting with the zeroth-order problem in Fig. 2 we assign wave coefficients a + 0 , a − 0 , . . ., d − 0 and their complex conjugates. Note that the subscript notation here is different from that used in [1] where subscripts were used to denote the harmonic index while they are used here to denote the ordered wave components of the perturbation analysis. These coefficients are related by propagation relationships, joint conditions, and clamped boundary conditions, where k denotes the wavenumber associated with the system's j th natural frequency, ω j . We point the reader to more details on wave-based approaches in [1], particularly as concerns joint conditions for the wave coefficients. Note that the zeroth-order stability problem is an eigenvalue problem. Next, we assemble the coefficient relationships into the form, where z (0) s denotes an eigenvector holding the eight zeroth-order wave coefficients (see Fig. 2

(b)) and A (0)
s ω j ; K 0 denotes the zeroth-order stability coefficient matrix as a function of the j th natural frequency ω j and parameterized by stiffness K 0 . Note that the dependence of A (0) s on frequency is nonlinear, and as a result, there are an infinite number of natural frequencies satisfying Eq. (40) consistent with the fact that the jointed system is continuous. While it is not apparent until the next order analysis, we also note that K 0 will be chosen based on requiring the eigenfrequency ω j to be equal to the forcing frequency ω. The zeroth-order problem is completed by finding the eigenvalues and eigenvectors associated with Eq. (40) [7].
Following completion of the zeroth-order problem, we turn our attention to the first-order problem. We note that the zeroth-order spring stretch required to evaluate Eq. (35) is now known and given by, where explicit dependence of ω j on K 0 is noted. Since the zeroth-order solution is an eigenfunction, we can freely choose the harmonic amplitudes a (0) and b (0) . Substitution of Eq. (41) into Eq. (35) yields an updated first-order forcing, We note secular, or resonant, terms on the right-hand side of Eq. (42) having frequency dependence equal to one of the zeroth-order system's eigenfrequencies; namely the underlined terms. These must be eliminated in order for the response to be bounded and to satisfy the assumed ordering of the asymptotic expansions. In addition, for particular choices of K 0 such that 2ω − ω j (K 0 ) = ω j (K 0 ), additional secular terms appear as indicated by the doubly underlined terms. In particular, the latter condition holds for a certain value of K 0 close, but not equal, to K from the forced problem. When this occurs, the excitation frequency ω equals the j th eigenfrequency ω j (K 0 ), as remarked during the zerothorder problem.
Proceeding with the choice of K 0 such that ω = ω j (K 0 ), we remove secular terms by requiring, Zeroing the determinant of the coefficient matrix in Eq. (43) yields K 1 , We note that the two values of K 1 in Eqs. (44) begin to form two sides of a stability tongue in the P versus K plane, as seen later. With known values for K 1 , it is also possible to find the eigenvectors, but these are not needed from this point forward.
With secular terms removed, we must still find the forced response associated with the non-secular terms in Eq. (42). As is the usual practice in perturbation approaches [4], we neglect the homogeneous (or unforced) solution in all orders above the zerothorder without loss of generality. With the choice of K 0 established, the non-secular forcing occurs at frequency 3ω j (K 0 ). We find the forced response of the O(ε 1 ) problem using another exact wave approach, as depicted by the wave coefficients in Fig. 2(c) and their complex conjugates. These first-order wave coefficients are related by propagation relationships, joint conditions, and clamped boundary conditions, where the joint conditions have been developed similar to the discussion in [1] with the exclusion of damping and the inclusion of imposed forces . Similar to the zeroth-order problem, we assemble the coefficient relationships into matrix form, where z (1) s holds the eight first-order wave coefficients, A (1) s 3ω j (K 0 ) denotes the first-order stability coefficient matrix, and f (1) s holds forcing terms proportional to P 4 a (0) and i P 4 b (0) . Unlike the eigenvalue problem at zeroth-order, Eq. (48) admits a single solution for the coefficients via inversion of A (1) s 3ω j (K 0 ) , when this inverse exists. For cases where this is not possible (see discussion of one such case in Sect. 4), it is implied that the strained system is resonant at 3ω j (K 0 ) in addition to having ω j (K 0 ) as a resonance. We will, however, not dwell on this aspect presently and proceed with the general formulation of the method assuming that A Following the solution of Eq. (48), the spring displacement required in the second-order problem can be written as, where we note that a (1) depends only on a (0) (and not b (0) ), and b (1) depends only on b (0) (and not a (0) ). In fact, a (1) . This results from the similarity in how a (0) and b (0) appear in the first-order joint relationships.
With the first-order problem complete, we turn final attention to the second-order problem. Substitution of ω = ω j (K 0 ) and Eqs. (41), (49) into Eq. (36) yields an updated second-order forcing, of which only the secular portion is of interest, Elimination of the two secular terms yields, where we introduce r (1)/(0) to represent the ratio of first-to zeroth-order spring displacement components. Lastly, we reconstitute the strained stiffness, recalling that K 0 is found by enforcing that the excitation frequency ω equals the j th natural frequency of the zeroth-order stability problem, ω = ω j (K 0 ). It is apparent from Eq. (52) that non-zero damping lifts the stability curves off of the K axis, and thus has a stabilizing effect. To return fully to assessing stability of the forced problem, we recall that K ≡ K + 3ε a * 2 2 and P ≡ 3 a * 2 2 , where a * and φ * denote the joint displacement's amplitude and phase for the periodic solution under consideration. We also note thatĈ and always appear with ε, and thus, no generality is lost in setting ε to one. Inverting Eq. (52) yields two curves defining the stability tongue [8] separating stable and unstable solutions in the P versus K plane (depicted in Fig. 7 for the single-jointed case).
We note that the strained parameter approach may also be developed for the two-jointed system studied in [1], where it is only necessary to strain one of the two joint stiffnesses, regardless of the fact that these two joints have the same stiffness-i.e., if they differed, it would still be only necessary to strain one or the other, and not both. This is due to the fact that the stability problem has codimension-1, i.e., only one resonant condition must be met, and thus, only one parameter must be strained [9].

Stability via the perturbed eigenproblem residue approach
The second stability approach is based on a PWE formulation of the perturbed system of the original problem. We consider a generic single jointed-bar vibration problem of the form and formulate the approach using this case as an example. Application to cases with multiple joints is straightforward since the method does not make assumptions on the number or types of nonlinearities present. In the above, f L (X, t) and f R (X, t) denote general forcing functions (independent of the solution u L ,R ). We note that this approach does not demand the presence of a small parameter, i.e., this approach is applicable for arbitrarily strong nonlinearities and not just weak nonlinearities (as was the case for the strained parameter approach presented in Sect. 2). Following the notation in Sect. 2, we use starred superscripts to denote the periodic solutions u * L ,R of the original forced vibration problem and δu L ,R to denote infinitesimal perturbations of the Substituting these expressions into Eqs. (53-58) and dropping the higher-order terms yields a homogeneous set of governing equations for δu L ,R . These can be expressed as Since u * L ,R (X, t) (and therefore u * , u * , . . . ) is fixed from the solution of Eqs. (53-58), the above is a linear homogeneous PDE with time-varying linear stiffness and damping acting at the joint, and the solution can be obtained using a wave-based approach. An important aspect in the homogeneous case, however, is that the frequency of the response is complex and unknown, unlike the inhomogeneous case where the frequency is real and fixed by the excitation. Resorting to a multi-harmonic wave-based representation of the solution, the solution δu L can be represented, for instance, around the joint as where k in the above comes from the dispersion relationship λ = √ E/ρk. Here, λ denotes the unknown complex frequency that the perturbed system responds at. Assembling this system and using the FFT to obtain the wave-component stiffness coefficient contributions for the periodic stiffness terms in Eqs. (63) and (64) above (which are 2π/ω-periodic), we obtain a linear set of equations of the form where z ∈ C N w N h is a complex vector consisting of the N h complex harmonics of the N w wave components; A ∈ C N r N h ×N r N h is a complex matrix that is a function of λ alone (but parameterized by the periodic solution u * , ω). Note that no AFT (Alternating Frequency -Time) procedure needs to be carried out in terms of the complex frequency λ here; we use the appropriate entries of the frequency-domain Jacobian of the PWE residue function (see [1]) to get the periodic stiffness components once. As expected from Floquet theory, this is a Nonlinear EigenProblem (NEP) in λ. However, since A contains exponentials of λ, it is not, in general, possible to spectrally decompose this problem quite easily and, furthermore, there exist an infinite set of eigenpairs. We unsuccessfully attempted the use of the Padé approximant and interpolatory Krylov techniques 1 for partial spectral decomposition of this NEP, but still feel that such approaches must be explored in detail in the future to enable quantitative stability analysis. Qualitative analysis of stability can, however, be carried out merely by seeking out the points along a given forced response curve where transitions of stability occur in a manner that completely avoids the challenges associated with spectral decomposition. Floquet theory asserts [4] that stability transitions occur when the perturbed system has as eigenvalue imξ , a purely imaginary quantity with integer m > 0, with three possibilities for ξ : ξ = 0; Perturbation remains unchanging in time; ξ = ω; Perturbation varies periodically with the same fundamental period as the solution; and ξ = ω/2; Perturbation varies periodically with twice the fundamental period as the solution.
Choosing ξ = 0 gives rise to the trivial solution and need not be considered in the context of continuum dynamics with homogeneous essential boundary conditions. In a multi-harmonic approach such as the wavebased strategies considered presently, m can be set to 1 without loss of generality since the higher harmonics will be intrinsically included for each ξ . In our context, therefore, the cases corresponding to ξ = ω and ξ = ω/2 will have to be considered. If ξ is an eigenvalue of the system in Eq. (68), the determinant of A λ=iξ will be exactly zero. Checking for points along any given forced response curve where this happens allows one to obtain the stability transition or critical points directly. Since A is complex, the determinant is, in general (for the cases when ξ is not an eigenvalue), also complex, making the task of detecting the zero point on a forced response curve challenging. We first therefore create a fully real counterpart of the Jacobian matrix (this process was previously alluded to in [1]). Denoting by {z} and {z} the real and imaginary parts of the vector of wave component harmonics z, Eq. (68) can be rewritten as Symbols denoting the functional and parametric dependences of A are dropped in the above for brevity. Decomposing the complex quantity Az into its real and imaginary parts denoted by {Az} and {Az}, respectively, yields the following fully real form of Eq. (68): The Jacobian matrix of this system, is fully real and its determinant will be real always. Along a forced response curve, a zero-crossing of the determinant of A] λ=iξ will be the same as a sign change in the determinant ofÂ λ=iξ . The quantity det (Â) λ=iξ will be referred to as the Perturbation Eigenproblem Residue (PER) henceforth. The signchange property will be clear if elaborated upon in the following manner: suppose the true eigenvalues of Eq. (70) are denoted by γ j with j being the index of the eigenvalue, we can write the PER as a characteristic polynomial in the following manner: Here |(.)| and (.) denote the absolute value and angle of their complex arguments respectively, and α is some real constant. A sufficient condition for the existence and validity of such a polynomial representation is the analyticity of the PER for all points λ on the complex domain. The existence of a Taylor series representation is guaranteed in this case, justifying the use of the polynomial form above. This is clearly the case here since λ only appears as smooth functions in the matrix A(λ) (and therefore in its determinant). This analyticity condition could be violated in some problems with nonsmooth nonlinearities (e.g., Filippov dynamical systems, see, for instance [10]). In these cases, classical Floquet theory is not strictly applicable and further considerations become necessary. We will presently not concern ourselves with such issues and proceed with the discussion making note of the fact that this framework is invalid for such cases. We also note here that there are no bounded values of λ for which the PER becomes unbounded; i.e., the characteristic polynomial does not have any poles on the complex plane. It, however, has an infinite number of zeros (despite the finite size ofÂ), referred interchangeably here as eigenvalues of the NEP, consistent with that expected in continuum dynamics.
Since the PER is always real, the effective phase ∞ j=1 λ − γ j is either 0 • or 180 • . Consider two points right before and after a transition in stability such as shown in Fig. 3. Just before the critical point (point A in Fig. 3), the zero γ (A) m closest to iξ will lie just to its left (since the solution is known to be stable) such that the quantity iξ − γ (A) m equals 0 • . Equivalently, the vector in the complex plane joining γ (A) m to iξ points to the right since this solution is known to be stable. The total angle of the PER at this point can be expressed as where the zero term is canceled in the final expression, which itself has to be either 180 • or 0 • since the PER is a real quantity. Looking at a point just after the transition (point B in Fig. 3), the zero γ (B) m has crossed through iξ and will now lie just to its right. Consequently, the quantity now takes a value of 180 • , and we can express the total angle of the PER at this point as the sum, Around these transition points, due to the arbitrarily small change in the excitation frequency ω and ensu-ing response amplitude, the zeros away from iξ are assumed to undergo only negligible changes and therefore γ is assumed for j = m. The first term in Eq. (74) is approximated as P E R A (see Eq. (73) above) and can be simplified as This shows that the difference in complex angle between the PER values across a critical point is exactly 180 • . Since the PER is always a real number, this amounts to a sign change as mentioned earlier. Similar reasoning can be employed to observe that a sign change will also happen at the other critical point in Fig. 3 where the solutions transition from unstable to stable. One drawback of the method presented is the fact that it can only predict stability changes and it is not possible, in general, to determine the stability of a solution from just the sign or the value of the PER. Exact stability determination would instead require knowledge of the location in the left or right half-plane for all infinite number of γ j , of which little can be said at this point. An exploration of the necessary conditions for the PER sign change (stability transition is just a sufficient condition) is deferred for future research. But the efficacy of the approach is empirically demonstrated by the fact that it provides satisfactory estimates of the stable and unstable branches on the response curve for both the numerical examples considered in [1].

Stability results for an example system
In this section we present stability results for the single-jointed example problem illustrated in Fig. 1 and parameterized by quantities appearing in Table 1. Periodic solutions to this problem are presented in [1], which also contains stability results for a second example problem consisting of three bars connected by two nonlinear joints.
The stability boundary predicted by the strained parameter approach for a fixed excitation frequency ω is a parabolic curve (known as a stability or Arnold tongue), in the P − K space, where P denotes the periodic stiffness and K the constant stiffness of the parametrically forced system (recall set to 1)-such a curve will be presented later. If, in addition, we vary the frequency, a stability surface results separating regions of stable and unstable parametrically forced systems. Table 1 Mechanical and geometric parameters used for the single-jointed bar example (see Fig. 1

for accompanying schematic)
Parameter Value E y 262 GPa ρ 1280 kg/m 3 K 10 9 N/m C 320 N/(ms −1 ) We present this surface in the subplots of Fig. 4 near the second natural frequency of the single-jointed system shown in Fig. 1, and label it as the Stability Boundarynote that this boundary is the same surface in all six sub-figures. We provide a narrated video in the supplementary material which includes revolving views of the figures and further explanations. The stability analysis is carried-out up to, and including, first-order terms. When we return to the directly forced system, specifying values for the joint linear and nonlinear stiffness, K and , respectively, together with the excitation frequency ω and response amplitude a * for a given periodic solution, yields a point in the K − P − ω space; see Eq. (13). If this point falls inside of the stability boundary, the solution is unstable, whereas as if it falls outside, the solution is stable. If in addition we sweep frequency, response curves result which lie in the same space -these appear in Fig. 4 as curves whose stable portions are indicated by solid black lines, and whose unstable portions by dashed green lines. A firstorder WBVA approach is used to obtain the solutions (see [1] for details) and a first-order strained parameter approach is used to obtain the stability boundaries. We present solution curves at forcing levels of 7.5 MN, 15 MN, and 30 MN, and we provide two views for each forcing level to aid visualization. It is important to note that while K is fixed along these curves, K is not, and thus the curves bend towards increasing K as resonance is approached and the response amplitude a * increases. If we also vary the linear stiffness K in the directly forced problem, we can generate a solution surface, shaded blue in Fig. 4 and labeled as the Solution Manifold. The solution curves and manifold occupy more space inside of the stability boundary as the forcing increases, as depicted in the progression of sub-figures. For example, at 7.5 MN forcing, the solution manifold never intersects the stability boundary, indicating all solutions are stable. In fact, the solutions for 7.5 MN forcing closely resemble the linear solutions such that no appreciable bending of the frequency response curve occurs, and thus multiple solutions near resonance have yet to appear. At 15 MN and 30 MN forcing, appreciable bending occurs together with multiple solutions and the solution manifolds intersect with the stability boundaries. For these forcing levels, the portions of the response curves associated with mid-amplitude response are now unstable. Lastly, we note that the intersection of the surface manifold with the stability boundary marks the Transition Boundary, as indicated by the solid magenta curves in the sub-figures.
We study the convergence of the stability tongues as the analysis order increases in Fig. 5. The sub-figure on the left presents the frequency response curve for the single-jointed system studied, evaluated with a forcing of 45 MN. At such a large forcing level, the ratio of nonlinear to linear restoring force at resonance is 0.25, well-above that which can be expected for solution accuracy. However, this forcing level yields three solutions whose stability can be used to test the accuracy of first-order versus second-order stability analyses. The three solutions are indicated in the left subfigure at 81.6 k rad/s using three markers. Subsequently, values for P and K for each are computed and indicated with the same markers in the right sub-figure, together with the stability tongue corresponding to 81.6 k rad/s. We note from the right sub-figure that the stability tongues match closely at low values of P, and thus, only first-order stability is necessary to assess the stability of solutions at low forcing, justifying the use of first-order analysis in Fig. 4. At higher forcing amplitudes, the mid-and high-amplitude solutions correspond to larger values of P. If first-order stability analysis is employed, the stability determination can be inaccurate, and thus second-order stability must be employed. This is clearly documented in Fig. 5(b) where the high-amplitude solution appears inside of the first-order stability tongue, erroneously suggesting unstable response. However, using second-order stability analysis, this solution subsequently lies outside of the stability tongue as the tongue appears rotated away from the first-order tongue. In our investigations of solution response and stability using the WBVA and the strained parameter approach, focusing specifically on the resonance peak where issues in stability pre-diction might arise, we have found first-order analysis is always accurate for responses in which the ratio of nonlinear to linear restoring force is less than 0.1.
Next we investigate the PER stability approach in the context of the single-jointed system. For all fre- quency responses considered, we set ξ equal to the excitation frequency ω in detecting critical points for which the sign of the PER changes. No critical points were found by setting ξ to ω/2, implying that stability loss occurs only through the ξ = ω case for the present system, consistent with the strained parameter approach. Figure 6 plots the frequency response and the perturbation eigenproblem residue corresponding to the PWE solution of the single-jointed system harmonically forced at an amplitude of 30 MN using both  First-and second-order stability boundaries from the strained parameter approach for the single-jointed system Fig. 8 Comparison of the forced response stability predictions for the single-jointed system a using the FE-HB approach with frequency domain Hill's coefficients [5] for stability determination, and b using the PWE frequency response and the PER stability approach developed in Sect. 3. In all figures, the stable and unstable branches are plotted using solid and dashed lines, respectively a linear and log scale. Specifically, we plot the PER, det (Â) λ=iξ , evaluated using λ = iω. The left subfigure shows the PER clearly changing signs at the two bifurcation points in the frequency response curve, correctly indicating stability transitions in agreement with the strained parameter approach. In the right sub-figure, the log value of the PER drops sharply at the bifurcation points, but interestingly, shows a drop and local minimum at approximately 79.8 k rad/s. The PER approach does not provide insight into the mechanism for the log value exhibiting a local minimum at approximately 79.8 k rad/s. Revisiting the stability boundaries predicted by the strained parameter approach, however, provides the requisite insight. Figure 7 presents the stability boundaries when we carryout the strained parameter approach to the first-and second-orders. At the second order, which also necessitates the inclusion of third harmonic terms, K 2 (and thus K ) goes unbounded at the same point in frequency where the PER exhibits a local minimum. Examining the coefficient matrix A(3ω; K 0 ) evaluated near the local minimum frequency, we find that it has a zero eigenvalue at 79.7636 k rad/s. Since the determinant of a matrix equals the product of its eigenvalues, it is apparent that 3ω at this frequency satisfies an internal resonance condition. Further investigation of the internal resonance mechanism lies beyond the scope of this paper, but may offer an interesting direction for followon efforts.
Lastly, Fig. 8 presents the finite element validation and PWE forced response solutions for the singlejointed system, where we assess stability of the PWE solutions using the PER approach. We note strong agreement between the two approaches in both the predicted response and accompanying stability. Notably, at a forcing level of 45 MN, the ratio of nonlinear to linear restoring force at resonance is over 0.7, which is well-beyond the validity region for the WBVA and the strained parameter stability approach. However, the PWE combined with the PER stability approach accurately captures both the frequency response and the multiple solution stability at all amplitudes considered. This can be further contrasted with the poor performance of the finite-element substitution approach discussed earlier and documented in Fig. 9(a).

Concluding remarks
As pointed out by previous researchers, a direct analysis of stability has been lacking in the literature for nonlinear wave-based solution approaches. We address this deficiency by developing two direct approaches, and then apply the methods to determine solution stability of an harmonically forced, example two-bar system connected by a nonlinear joint. The first stability approach expands (or strains) a stiffness parameter in the problem, which allows one to construct analytical stability boundaries in the parameter space of the perturbed problem using a wave-based method. Not only does this determine stability of periodic solutions arising from direct excitation of a continuous system, but it also directly governs the response and stability of parametrically forced continuous systems, which is a new development in and of itself. The second stability approach is primarily a computational tool which searches for sign changes in the residue of the per-turbed eigenproblem. The strained parameter approach holds for only moderate nonlinear response due to its reliance on an asymptotic expansion; the PER avoids this limitation and in fact is used to determine stability of periodic solutions at large nonlinear response.
We now list some avenues for future research that were identified during the course of this study. One issue with the PER approach, as presented herein, is that it can only determine locations of stability change, and does not assess stability directly. Therefore, a PER approach that can assess the latter warrants investigation. Further, the PER method described in this paper only involves a sufficient condition for the sign change. An investigation of the necessary conditions could potentially yield greater insights. Both stability approaches (strained parameter as well as PER) rely on Floquet theory, whose applicability for non-smooth problems is not very clear at this point. This will have to be studied in detail (see [10], for instance) to understand the theoretical limitations of the methods more clearly. The internal resonance observed in the singlejointed problem (see Fig. 6 and related discussions) should also be given additional attention.

Data availability
The datasets generated during and/or analyzed during the current study are not publicly available due to limitations on archiving data long-term, but are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.

Appendix A
In this appendix, the PWE approach [1] is used to obtain the frequency responses of the single jointed system (Fig. 1) around its first resonance. We then project these solutions onto a 90-element finite element mesh, for use as initial conditions, and estimate its monodromy matrix using transient simulations over a single period. The predicted stability using this hybrid approach is indicated in Fig. 9(a). Comparing these stability predictions with pure finite element stability predictions ( Fig. 9(b), which use frequency domain Hill's coef- Fig. 9 Comparison of the forced response stability predictions for the single-jointed system a using the PWE-computed frequency responses with stability assessed using transient Floquet multipliers estimated from the FE model, and b the FE-HB-computed frequency responses using the frequency domain Hill's coefficient approach [5] (reproduced from Fig. 8). In all figures, the stable and unstable branches are plotted using solid and dashed lines, respectively ficients [5]), shows that although the response curves are nearly the same, the stability predictions contradict, specifically around the peaks at large forcing levels. This shows that the hybrid approach can lead to incorrect conclusions about solution stability, and motivates the need for direct wave-based stability approaches.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.