Wave-based analysis of jointed elastic bars: nonlinear periodic response

In this paper, we develop two wave-based approaches for predicting the nonlinear periodic response of jointed elastic bars. First, we present a nonlinear wave-based vibration approach (WBVA) for studying jointed systems informed by re-usable, perturbation-derived scattering functions. This analytical approach can be used to predict the steady-state, forced response of jointed elastic bar structures incorporating any number and variety of nonlinear joints. As a second method, we present a nonlinear Plane-Wave Expansion (PWE) approach for analyzing periodic response in the same jointed bar structures. Both wave-based approaches have advantages and disadvantages when compared side-by-side. The WBVA results in a minimal set of equations and is re-usable following determination of the reflection and transmission functions, while the PWE formulation can be easily applied to new joint models and maintains solution accuracy to higher levels of nonlinearity. For example cases of two and three bars connected by linearly damped joints with linear and cubic stiffness, the two wave-based approaches accurately predict the expected Duffing-like behavior in which multiple periodic responses occur in the near-resonant regime, in close agreement with reference finite element simulations. Lastly, we discuss extensions of the work to jointed structures composed of beam-like members, and propose follow-on studies addressing opportunities identified in the application of the methods presented.


Introduction
Interest in the dynamics of jointed structures, in which mechanical joints connect elastic members, has grown during the past decade [1,2].This interest has been spurred in part by new insights into the influence that joint dynamics and tribology have on system-wide dynamic response and hysteretic damping [2].Owing to the complexity of such problems, nearly all predictive methods employ high-fidelity, discretized models (e.g., finite element models) to compute dynamic response -see, for instance, [3,4].Elastic members away from the joint (trusses, beams, etc.) are accurately modeled using linear representations, while approximate nonlinear representations are typically employed to represent frictional contact behavior at the mechanical joints.Due to the small domain size of the joints in comparison to the members, these methods allocate the majority of their computational effort to linear components, while allocating comparatively smaller computational effort to the intricate nonlinear dynamics of the mechanical joint.Herein, we address this apparent paradox through the development of a nonlinear wave-based analysis approach in which known wave solutions represent elastic members, effectively eliminating their computational burden, while mechanical joints are represented by discrete, nonlinear models whose frequency-and amplitudedependent reflection and transmission behavior is first found using perturbation approaches.The present research represents a first step in this direction, and as such, we consider the elastic members to be simple bars, and we represent the joints by nonlinear springs with linear damping.To place the paper in its proper context, we first review relevant literature on joint dynamics and damping, wave-based vibration analysis techniques, and lastly, scattering of elastic waves from nonlinear interfaces Some of the earliest modern efforts on characterizing the nonlinear dynamic effects of mechanical joints in jointed structures can be found in investigations of the dynamics of the so-called Gaul resonator [5,6], which connects two masses using a bolted lap joint.The dynamic response of such jointed structures, at the interface level, are characterized by intermittent periods of slipping, where the interfaces slide over one another with non-conservative force-displacement characteristics, and sticking, characterized by conservative force-displacement characteristics (instantaneously).Under harmonic operation this manifests as the well-known nonlinear hysteretic characteristics of the jointed structure.Practically, these structures demonstrate an amplitude-dependent softening (reduction in stiffness/resonant frequency with increasing amplitude) and dampening (increase in the effective damping factor with increasing amplitude) as have been observed in many experimental studies (see, for instance [7][8][9][10]).Mathematical modeling of hysteresis has been an active area of research over the past century (see reviews such as [11,12]).These range from simple saturating springs (like the Jenkins model [13]), to models based on structural properties of hysteresis (like the Masing hypothesis [14]), and more involved mathematical formulations (such as the Bouc-Wen model [15]).The authors in [12] present a recent review detailing such models and their influence in a steady-state setting.
Researchers in the structural dynamics community have recently applied known propagating wave solutions (as opposed to traditional standing wave solution techniques) to analyze exact steady-state dynamics of linear elastic structures, which requires knowledge of wave scattering at discontinuities [16][17][18][19].In frame-like and/or lattice structures composed of one-dimensional members (e.g., bars and beams), the wave-based vibration approach (WBVA) [17,18,[20][21][22][23][24][25][26][27][28] exploits known wave solutions propagating in the forward and backward directions along the elastic members, together with characterization of reflection and transmission at connections and boundaries, to obtain exact system eigenfrequencies and eigenfunctions under free motion conditions, or steady-state response under harmonic excitation.The present-day wave-based vibration approach owes its existence to the culmination of studies that first focused on developing reflectivity and transmissibility relationships for a variety of boundary conditions and interfaces using bar and beam (e.g., Euler-Bernoulli and Timoshenko) theories [17,21,22,29,30].These studies related the incident wave coefficients to their reflected and transmitted counterparts, typically using matrix representations.Armed with the knowledge of reflection, transmission, and excitation-generated wave coefficients, researchers then posed closed, linear algebraic conditions governing exact solution of the wave coefficients as applied to two-and three-dimensional structures (e.g., frames) [20,24,31,32], for both free and harmonically-forced motions.It is important to note that in these papers the connections between members are commonly referred to as 'joints' -however, these are more identifiable as rigid connections since they lack both compliance and damping.Herein, we exclusively refer to this procedure, of determining vibration response using an assemblage of wave coefficient conditions, as the wave-based vibration approach.Related, but fundamentally different, methods not grouped under the WBVA designation include the wave finite element approach [33] and spectral wave analysis [34], among others.Extensions of the exact WBVA have been explored by Leamy for the study of Bloch waves in periodic two-dimensional lattices [27], and by Lv and Leamy to study boundary-damped systems [28].Lv and Leamy also describe an automated assembly procedure and implementation in Matlab, similar to that used in finite element assembly procedures, which generates the problem solution simply by user specification of the material properties and geometry, connectivity of the members, and forcing properties.Implementation of the WBVA in structures incorporating two-dimensional elastic members, such as membranes and plates, is an open research problem.However, one basis for such a technique can be found in [35], in which wave solutions are used to find exact and nearly-exact eigenfunctions of rectangular plates supported by a variety of boundary conditions.
Wave scattering at nonlinear joints has received sparse attention in comparison to linear connections, and as such a WBVA incorporating such joint models has yet to be presented.Gaul [16] studied wave transmission and energy dissipation through joints using a variety of discrete nonlinear representations, but simplified the analysis using equivalent elements (i.e., linear springs and dashpots) for representing the constitutive features of the joints.Improving upon this, Vakakis [19] studied nonlinear scattering of longitudinal waves in a rod with a weakly nonlinear elastic joint at its end using a perturbation approach, and applied the phase-closure principle to compute vibration backbones.Later, Brennan et al. [36] derived the reflection and transmission coefficients for both rods and beams terminated by a cubic spring using harmonic balance.Other researchers have developed a multi-harmonic formulation for the description of the interaction of waves in linear waveguides through nonlinear coupling elements [37], and extended it using a numerical implementation of the harmonic balance method to estimate the forced steady-state response of linear frame structures with nonlinear supports [38].Most recently, [39,40] presented a general perturbation approach for determining scattering in a beam due to an array of attached nonlinear oscillators.
We first aim to extend the WBVA to analyze weakly nonlinear systems, while finding the requisite reflection and transmission behavior at nonlinear joints using a perturbation approach.Importantly, the perturbation approach considers coupling of incident waves and their anticipated higher harmonics, which first arise at the joint and then reappear incident after downstream reflections.This is followed by a critical reconstitution step before usage in the nonlinear WBVA.Although demonstrated for longitudinal waves in multiple jointed bars, the nonlinear WBVA is general and can be extended in a straight-forward manner to beams and frames joined by a variety of discretely represented, nonlinear joints, which will be the subject of future work.Concurrently, we implement an alternating frequency-time domain plane wave expansion (PWE) approach for analyzing the same jointed bar systems, and then make comparisons between the two approaches for harmonic loading.We validate the results using finite element models evaluated by the harmonic balance method.While we present stability results for the example systems studied, development of the stability methods used can be found in a companion paper [41].
The paper is organized as follows: in Sec. 2 we first develop requisite scattering relationships for a damped, nonlinear joint between two semi-infinite bars using perturbation analysis.We then formulate the nonlinear WBVA for finding system dynamic response using the example of two clamped bars connected by the damped, nonlinear joint; later, we consider three bars connected by two joints.In Sec. 3 we detail the PWE approach and its application to the same two-bar, clamped system.This is followed by Sec. 4 where we present forced response results from the developed approaches together with finite element-harmonic balance results for validation purposes.The paper concludes with discussion of the results and avenues for future work.

Nonlinear Wave-Based Vibration Approach
In this section we develop the nonlinear wave-based vibration approach through regular perturbation theory.We start by developing scattering relationships corresponding to single and multiple wave incidence on a joint of a system of linear-elastic bars.We then apply the scattering relationships to a wave-based vibration problem to predict system response resulting from harmonic loading.

Scattering of a Single Incident Wave
As depicted in Fig. 1, we consider two semi-infinite, linearly elastic bars joined by a nonlinear spring in parallel with a linear damper.The spring is characterized by linear and cubic stiffness coefficients K and εΓ, respectively.The bars support longitudinal motions only such that the left and right bars have displacement fields denoted by u L (x, t) and u R (x, t).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. Generalization to two bars of unequal material and geometric properties is straight-forward.We choose to consider this system as it represents the simplest setting to study wave propagation inside a jointed structure in which the joint is modeled using a nonlinear stiffness element and a linear damper.Consideration of hysteretic joints is deferred for future work.The equations governing wave propagation in the linearly-elastic bars are given by, where f L (x, t) and f R (x, t) denote forcing per unit length in the left and right subdomains, respectively.We will consider forcing later; in the absence of forcing, the bars admit freely propagating waves.We model a single-frequency incident plane wave radiating from the left boundary with complex amplitude a, I (x, t) = ae ikx e −iωt + c.c., which then reflects off, and/or transmits through, the joint located at x = 0. Note that we denote the complex conjugate of all preceding terms by +c.c.such that the total incident field is real (as required for nonlinear analysis).Since the bars carrying all waves are linear and governed by Eqs. ( 1) and ( 2), the wavenumber k is related to the frequency ω by ω = c p k, where c p = Ey ρ denotes the wave speed.The reflected and transmitted waves are both a function of space x and time t, which we denote by R (x, t) and T (x, t), respectively.For a linear problem, these fields would simply be of the form [Rae −ikx e −iωt + c.c.] and Tae −ikx e −iωt + c.c. , respectively, where R and T denote complex reflection and transmission coefficients to be determined by suitable joint conditions.However, due to the nonlinear spring stiffness, these forms do not suffice.Four boundary conditions are required to completely specify the wave propagation problem consisting of two jointed bars.Radiation and causality conditions provide two conditions: the presence of a leftward propagating incident wave in the first bar with complex amplitude a, and the absence of a rightward propagating wave from the right in the second bar.Two other boundary conditions arise at the joint, which state that the internal force at x = 0 must equal the action of the spring and dashpot, Satisfaction of these boundary conditions requires an asymptotic approach predicated on the smallness of the parameter ε.Assuming ε ≪ 1 and positive, we decompose the displacement fields into orders of ε, similar to the approach described in [40], where 6) and (7) into Eqs.(1), ( 2), ( 4) and ( 5), followed by collection of terms with the same powers of ε, yields the ordered field equations and boundary conditions, O ε 1 : In addition to the boundary conditions appearing above, the incident radiation condition I (x, t) is present in u L0 (and not u L1 or u L2 ) since Eq.(3) appears at O ε 0 , and causality dictates that no leftward-moving waves exist in either u R0 or u R1 .

Zeroth-order solution
The O ε 0 problem recovers a linear problem whose solution we seek in the form, where the assumed solution forms incorporate the radiation and causality conditions.Substitution of Eqs.(20) and (21) into Eqs.(10) and (11) provides algebraic equations for solution of the reflection and transmission coefficients at O ε 0 , yielding With the O ε 0 solution complete, we next seek solutions to the O ε 1 problem.

First-order solution
Substituting the known solutions for u L0 (x, t) and u R0 (x, t) into the O ε 1 problem yields updated boundary conditions, where coefficients c 1 and c 3 are given by, with overbars denoting complex conjugate and where the identity R 0 + T 0 = 1 has been used to simplify the expressions 1 .Due to the presence of the third harmonic (i.e., e i3kx e −i3ωt ) on the righthand side of the updated boundary conditions, we seek solutions to Eqs. ( 12) and ( 13) decomposed into the fundamental harmonic and third harmonic components, 1 ae ikx e −iωt + T (3) where R (1) 1 and T (1) 1 capture higher-order corrections to reflection and transmission of the fundamental harmonic, and R  28) and ( 29) into Eqs.( 24) and ( 25), followed by separation of the orthogonal harmonics (i.e., ω and 3ω), yields the necessary conditions to find the four coefficients.We then find the resulting amplitude-dependent expressions, T (1) 1 For non-zero damping we have R0 + T 0 = 1 and

Second-order solution
The procedure for the second-order and higher follow that detailed for the first-order; namely, the boundary conditions are updated with all known quantities from lower orders, solutions are assumed capturing the harmonics present in these updated equations, and reflection and transmission coefficients are then determined.For the second-order, this results in, 2 ae ikx e −iωt + T 2 ae i3kx e −i3ωt + T (5) where

Solution reconstitution
Based on the expansions in Eqs. ( 6) and ( 7) and the solutions found at each order, we can reconstitute the total displacement fields in the left and right bars, u R (x, t) = T (x, t) = T (1) ae ikx e −iωt + T (3) ae i3kx e −i3ωt + T (5) ae i5kx e −i5ωt + c.c.
where the total reflection and transmission coefficients are given by, We note that ε and Γ always appear together, and thus there is no loss in generality in setting ε to one once the final expressions are obtained.

Scattering of Multiple, Ordered Incident Waves
Vibration problems in which a single frequency excitation acts on two finite-length jointed bars motivates the study of scattering depicted in Fig. 2. In such a vibration problem, as we will study in Sec.2.3, application of the excitation source will first generate an incident wave at the excitation frequency ω, which upon encountering the joint, will result in reflected and transmitted waves at the excitation frequency and its multiples, as found in Sec.2.1.Since the vibration problem consists of finite-length bars, each of these reflected and transmitted waves will once again return after they reflect from the opposite ends of the bars.Thus, as in Fig. 2, in the steady-state there will be incident waves on the joint from both sides composed of zeroth-order waves at the fundamental (or excitation) frequency, and higher-order waves at the higher harmonics.We note that the solutions for each incident wave acting separately cannot be superimposed due to the nonlinear joint, and thus a combined treatment is necessary.We model multi-harmonic incident plane waves radiating from the left and right boundary, which then reflect off, and/or transmit through, the joint located at x = 0. Note that in this case we cannot distinguish a sole reflected wave in one half, and a sole transmitted wave in the other.Instead, we recognize that both reflected and transmitted waves exist in both halves, which we term RT L (x, t) and RT R (x, t) for the left and right halves, respectively.For the cubic spring shown in Fig. 2, based on the analysis of Sec.2.1, the higher-harmonic coefficients we need to retain for an analysis up to and including the second order are a 3 , a 5 , b 3 , and b 5 such that the incident left and right expressions simplify to, The ordered field equations and boundary conditions remain those studied earlier, namely Eqs. ( 8) to (19).Due to the linearity of the zeroth-order problem, its solution consists of a super-position of two incident wave solutions, where R 0 and T 0 are again given by Eqs. ( 22) and (23).For the first-order problem, an incident wave now appears in each bar at the third harmonic, and the solution again contains the fundamental and its third harmonic, but reflection and transmission behavior cannot be separated.Thus we seek solutions of the form, 1R e ikx e −iωt + RT 1R e i3kx e −i3ωt + c.c., where RT 1L , RT 1R , RT 1L , and RT 1R denote functions of the wave coefficients a, b, a 3 , and b 3 capturing the aforementioned combined reflection and transmission.Note that due to coupling, these functions cannot be written as coefficients multiplying any one of the wave coefficients, as done in Eqs. ( 54) and (55).Similar to the solution procedure detailed in Sec.2.1, which requires updating first the O(ε 1 ) field equations and boundary conditions with the O(ε 0 ) solutions, we find the requisite functions in closed form, where denote R 0 , T 0 evaluated with ω → nω and k → nk, and we indicate explicit functional dependence on the wave coefficients in Eqs.(58) to (61).The wave coefficients a and b appear together and cannot be separated, as expected.
We then follow the procedure once more to obtain the second-order solutions.We thus seek solutions of the form, 2R e ikx e −iωt + RT 2R e i3kx e −i3ωt + RT 2R e i5kx e −i5ωt + c.c., where 1R (a, b, a 3 , b 3 ) , (64) With assistance from computer algebra, the procedure can be carried-out to arbitrarily high orders.However, for the purpose of this paper, we will restrict our analysis up to the second order accurate perturbations (up to and including the O(ε 2 ) terms).Lastly, we reconstitute the solutions to obtain final expressions for the total reflection/transmission and displacement fields, where I L and I R are given by Eqs. ( 52) and (53) and the total reflected/transmitted fields are given by, 2R e ikx e −iωt + εRT where we suppress dependence on wave coefficients in the reflection/transmission functions for sake of brevity.We note that in the final expressions ε always appears together with Γ, a 3 , a 5 , b 3 , and b 5 , and thus no loss in generality results from setting ε to one.

Nonlinear WBVA
In this section we pose and solve a nonlinear wave-based vibration approach capable of finding the forced response of two finite, jointed bars excited by single-frequency forcing of the form 1 2 F 0 e −iωt + c.c..The system of interest is depicted in Fig. 3.The three values for α i , i = 1, 2, 3, form a partition of unity.Note that ε no longer appears in the problem, and instead we consider only values of F 0 , K, and Γ such that the ratio of cubic restoring force to linear restoring force is small -we check adherence to this condition after computing the response (see Secs. 4.1 and 4.2).
To date, the wave-based vibration approach has been limited strictly to linear systems.We extend the approach to structures in which nonlinearities arise solely at discontinuities (e.g., joints and other connections).While the members considered are modeled as bars and not beams, the approach can be extended in a straight-forward manner to Euler-Bernoulli and Timoshenko beams.This will be the focus of follow-on work.The wave-based vibration approach takes advantage of known, exact wave solutions in the structural members, together with wave conditions at discontinuities (i.e., boundaries and loading locations), to yield exact analytical solutions to complex, multi-connected systems.The utility of such an approach is that it (i) yields solutions that would otherwise require discretization to obtain, and (ii) enables significant re-use after formulating the scattering behavior of the joint.Fig. 3 depicts the jointed-bar forced vibration problem whose solution we seek using a nonlinear WBVA.We identify wave coefficient vectors (i.e., a + , a − , b + , b + , ..., f − ) at the locations of discontinuities, which here include boundaries, the forcing location, and the jointed connection.The coefficient vectors hold individual coefficients corresponding to the fundamental, third, and fifth harmonic frequency -e.g., a + = [a + 1 a + 3 a + 5 ] T .As per Secs.2.1 and 2.2, this is consistent with carrying-out the perturbation approach up to, but not including, O(ε 3 ).Each individual wave coefficient holds the magnitude and phase of either a leftward or rightward moving wave at each location identified.These coefficients can then be used to represent the bar displacement along the bar until the next set of coefficients is reached.For example, for all points to the right of the forcing, and up to the joint, the displacement is given by, Note from Eq. ( 74) that the local origin for the displacement field begins with x equal to zero at the location identified for the coefficients -i.e., immediately to the right of the forcing for coefficients c 1 , c 3 , c 5 , and their complex conjugates.We next list the equations required to find the vibration response of the nonlinear problem illustrated in Fig. 3.In regions of the bars between discontinuities, coefficients bounding the continuous domain are related by propagation relationships, where j = 1, 3, 5.Note that we suppress expressions for complex conjugate coefficients since it is only necessary to solve the problem for the original quantities and then later append to the solution its complex conjugate.For clamped ends, the zero-displacement boundary conditions yield, At the forcing location, a jump discontinuity occurs for the fundamental frequency coefficients, but not for the third-and fifth-harmonic coefficients, documented in Secs.4.1-4.2.We note further that if the scattering problem in Sec.2.2 is carriedout to third-order (and higher), the fifth harmonic (and higher) would also appear in the analysis of the fundamental harmonic.
Eqs. (75) to (84) hold 36 nonlinear algebraic equations for 36 unknown wave coefficients, and represent the nonlinear WBVA problem formulation.A large subset of the equations, totaling 30, are in fact linear.As such, Eqs.(75) to (78) can be written as, where A (z red ; P, ω) denotes a sparce coefficient matrix parameterized by material and geometric properties from the set P = {E y , ρ, A, α 1 , α 2 , α 3 , l}, frequency ω, and a function of individual wave coefficients from the reduced coefficient vector 4EyAk .Due to the sparsity of A (z red ; P, ω) and the small problem size, Eq. ( 85) returns relatively simple expressions for the elements of z in terms of z red .We find these expressions using computer algebra -see Appendix A. Following substitution of the now-eliminated wave coefficients z into Eqs.(79) to (84) results in six nonlinear algebraic equations for the six unknown wave coefficients stored by the reduced coefficient vector z red , which represents a significant reduction in problem size.After specification of the parameter set P, and frequency ω, we numerically compute solutions for z red and then compose the total displacement field from the coefficients, j e ijkX e −ijωt + a − j e −ijkX e −ijωt + c.c., 0 < X < α 1 l c + j e ijk(X−α 1 l) e −ijωt + c − j e −ijk(X−α 1 l) e −ijωt + c.c., α 1 l < X < (α 1 + α 2 )l e + j e ijk(X−(α 1 +α 2 )l) e −ijωt +e − j e −ijk(X−(α 1 +α 2 )l) e −ijωt + c.c., (α where X denotes a global position coordinate with zero location at the leftmost edge of the domain.

Nonlinear Plane-Wave Expansion Approach
We now present an approach for the forced response estimation problem that is closely related to, but fundamentally different from, the nonlinear wave-based vibration approach presented above.Eq. ( 86) expresses the solution of the problem shown in Fig. 3 through the individual directed wave components.The same solution form is assumed for the second approach, referred to henceforth as the Plane Wave Expansion (PWE) approach.The point of departure of this approach from the above nonlinear WBVA is the treatment of the nonlinearities at the joint.While the WBVA treats it from a perturbation standpoint, leading up to the derivation of re-usable, nonlinear multiharmonic scattering relationships, the PWE approach employs a numerical Galerkin projection approach applied to the entire problem in time, while retaining the traveling wave formalism in space.An approach with some similarities was proposed and demonstrated recently in [38,42].However, herein, we do not discretize with traditional interpolation functions and we use only waveform basis functions.Dropping ε, we write the joint force relationship in Eq. ( 4) as, , or more generically, In the above, the nonlinear internal force generated from the joint is represented as f nl , with the relative displacement at the joint location written as ∆u, while ( ) denotes a partial derivative with respect to time.Since the PWE approach is sufficiently generalizeable, we will use this form of the force balance at the joint for the remaining treatment.
Restricting the present discussion to modeling the joined members as bars, all the terms in the solution expansion (Eq.(86)) are harmonic; i.e., we can decompose the solution into complex harmonic wave components via a Fourier series.Evaluating this at any location (here attention is only on the joint), will result in a time-periodic function with fundamental period 2π/ω.Using a similar notation as in previous sections, we write u L , ∂u L /∂x (local to the joint), and ∆u(t) as where coefficients δ j are known functions of d + j , d − j , e + j , and e − j .We note that this is an assumed form of the solution, while Eq. ( 86) is a derived solution form found asymptotically.Notably, the asymptotic solution reveals that only odd harmonics participate in the response (due to the cubic nonlinearity), which is not strictly known a priori for the PWE.
The nonlinear function f nl evaluated with the time-periodic ∆u(t) will result in a time-periodic restoring force, whose harmonic coefficients we denote as f A Galerkin projection of Eq. (91) onto the Fourier basis leads to nonlinear, coupled algebraic equations in terms of the wave coefficients, The nonlinear force harmonic coefficients {f (nl) j } j=0,1... can be numerically calculated using the Alternating Fourier-Time (AFT) approach, schematically represented in Fig. 4. The displacement harmonics {δ j } j=0,1... (derived from the wave components) are first transformed to the time domain (using an inverse FFT), yielding ∆u(t), since the nonlinearity is usually defined in the time domain in a more convenient form.Using the definition of the nonlinearity, f nl (∆u, . . . ) is evaluated (along with its gradients).Following this, the nonlinear force harmonics can be calculated by an FFT operation, transforming the force (and its gradients) from the time domain to the frequency domain.
The boundary conditions may be treated in a manner similar to the nonlinear WBVA.For a clamped end at the left end of the structure in Fig. 3, where the local displacement field is given as the zero displacement condition (at x = 0) simplifies in the following fashion (also see Eq. ( 76)): In other words, the boundary conditions are written at each harmonic level as algebraic relationships between the wave-components.Given a vibration problem such as the one in Fig. 3, the wave component coefficients are related by linear propagation relationships at each harmonic order (j in the above, see Eqs. ( 75) to (78)).Appending the above algebraic relationships (force balance at the joint(s), Eq. ( 92), and clamped boundary condition(s), Eq. ( 94)) to the wave propagation relationships yields a nonlinear algebraic system of the form Supposing the number of wave components in the problem is N w and the number of harmonic components is N h , z ∈ C NwN h is the vector of wave components and their harmonics; F ∈ C NwN h is the constant vector with the inhomogeneous terms, and F is a vector consisting of the homogeneous terms in the system (linear propagation and joint nonlinear force balance).Setting the residual r ∈ C NwN h , using z, to zero provides the PWE solution to the vibration problem.The usual practice to study steady state response is to set all the elements of F to zero except for the index corresponding to the first harmonic of the excitation location, which will be set to the forcing amplitude.This allows us to capture all the super-harmonics that are generated from such an excitation.If one is interested in the subharmonics, for instance the 1/3 rd subharmonic, the force will be applied at the third harmonic (so the non-zero element of F will be the index of the third harmonic of the excitation point).Wave components at the first and second harmonic terms from such an analysis corresponds to the 1/3 rd and 2/3 rd subharmonics respectively.The former approach will be taken for all the numerical studies in Sec. 4 for simplicity, but the validity of neglecting subharmonic terms for the numerical examples considered is verified through time integration (see App. C).
The development of the above approach is closely related to the traditional Harmonic Balance Method (HBM) which, in recent years, has become extremely popular in the nonlinear structural dynamics community (see [43] for a detailed presentation and review).As a consequence, estimation of the analytical Jacobian of the residue in the above (Eq.( 92)) is a straight-forward exercise, enabling the use of any gradient-based solver to obtain the PWE solution.It must be noted that this approach is only suitable for the study of periodic responses.Although extensions are possible for quasi-periodic cases (see [44], for instance), it is not applicable when the responses are aperiodic (e.g., chaotic).

Periodic Response and Validation
We now apply the two wave-based solution techniques to two examples of one-dimensional bars connected with nonlinear joints and forced harmonically.The indicated stability is found using the techniques developed in [41].Additionally, a discretized Finite Element (FE) model is developed for the two examples for validation purposes.We use the traditional nonlinear multi-harmonic balance method [43] to predict the forced response of the FE models, with stability estimated using the frequency domain Hill's coefficient estimation approach [45].Unless otherwise indicated, we retain the first five harmonics for all the harmonic balance simulations (FE and PWE).The sufficiency of five harmonics is verified by conducting simulations around each resonance with seven harmonics and comparing the relative magnitude of the seventh harmonic with respect to the first and ensuring it is below a threshold (10 −3 used here).Time-domain validations are also conducted around the resonances as a secondary check.Similar steps are undertaken to ensure that the second-order perturbation is sufficient to capture the response.App.C presents numerical results and discussions for one particular case from the examples below.No viscous dissipation is assumed in the linear bars; it is, however, straightforward to include the effects of a linear viscous term in the wave-based formulation using the developments above with complex dispersion relationships (see, for instance [46]).

Bars Connected with a Single Nonlinear Joint
The first example we consider is the single-jointed bar system already used for the development of the nonlinear WBVA in Sec.2.3 (see Fig. 3 for a schematic).Table 1 presents the numerical values of the parameters used for the analysis.A linear finite element model with 90 equal-length C 0 elements is used for validation.The nonlinear joint is represented as a nonlinear phenomenological relationship between appropriate nodal degrees-of-freedom (DoFs).
Before considering the nonlinear forced response, we first determine the eigenmodes of the linearized system (Γ = 0), which aids in identifying near-resonant frequencies in the nonlinear model.Applying the traditional WBVA to this problem in the absence of forcing yields an eigenvalue problem similar in form to Eq. (85), A(ω; P)z = 0.
Since A(ω; P) has nonlinear dependence on ω, this is mathematically a Nonlinear Eigen-Value Problem (NEP or NEVP).Further, since matrix A is, in general, complex, its fully real counterpart, Table 1: Mechanical and geometric parameters used for the single-jointed bar example (see Fig. 3 for accompanying schematic).
is constructed (see [41] for a derivation) and the determinant of Â is taken as the real-valued characteristic polynomial that takes on zero when ω equals an eigenvalue.The characteristic polynomial of this structure is shown in Fig. 5a, which also indicates the natural frequencies using markers.In addition, Fig. 5b-g show the deflection shapes corresponding to the resonant modes.Since this problem is linear, the WBVA yields exact eigenvalues and eigenfunctions.
With the eigenmodes identified, we use the developed approaches to calculate the nonlinear response at frequencies close to the first three natural frequencies.Fig. 6 documents the forced response results (first, third and fifth harmonic amplitudes and first harmonic phase) corresponding to forcing the system harmonically with an amplitude of 7.5 MN near the first eigenfrequency.The response at a point immediately to the right of the joint is used for all frequency response plots.The unstable regimes of the solutions are plotted using dashed lines.Stability for the finite element solution is determined using the frequency-domain Hill's coefficient estimation method [45].For the wave-based approaches we determine stability using the approaches developed in [41].Specifically, the strained parameter approach is used for the nonlinear WBVA solutions, and the perturbation eigenproblem residue approach is used for the PWE solutions (although either can be used for both with no appreciable differences).
At the considered excitation level, the ratio of the peak nonlinear force to linear force at the joint is approximately 0.032, indicating nonlinear response well within what is often considered weakly nonlinear (i.e., 0.1).The zeroth-order WBVA solution recovers the linearized response, and as such, does not contain higher harmonic content or multiple solutions.The first-order WBVA solution already provides a very good approximation of the response until the third harmonic (recall, fifth harmonic terms are not generated at the first-order).At this order, the response exhibits classical frequency response bending and multiple solutions associated with Duffing [47] and other nonlinear ordinary differential equations.We note that the WBVA captures this behavior using straight-forward perturbation, without the need for more specialized treatments such as averaging or multiple scales, since secular terms do not arise in the formulation.The second order WBVA, PWE, and finite element (FE-HB) results match closely, up to and including the third harmonic, with small differences around the peak amplitude occurring in the fifth harmonic.Fig. 7 presents the first harmonic amplitude and phase results of the system for the second and third nonlinear resonances, where the harmonic excitation amplitude is set to 15 MN.Here, the peak nonlinear-to-linear force amplitude ratios are approximately 0.061 and 0.0017, respectively.While the nonlinear response near mode 2 is comparable to that near mode 1 for 7.5 MN excitation, there is very little nonlinear response near mode 3 for this excitation level.Consequently, we observe that unlike the first two modes, the third mode behaves linearly, as evidenced by the fact that the O(ε 0 ) solution already provides a very good estimate of the converged response.Otherwise, we can make similar observations as made for mode 1.   while the new right end adds a clamped boundary condition, and a single set of propagation conditions are added, All other relationships from Sec. 2.3 can be re-used.The finite element model for this case is constructed using 120 equal-length C 0 linear elements with the two discrete nonlinear joints connecting appropriate nodal DOFs as in the previous case.
As before, we start with a linearized analysis to determine the eigenmodes before investigating the nonlinear forced response.Fig. 9 depicts the characteristic polynomial and the first six resonant frequencies and mode shapes.The two joints in the system are highlighted with blue and red respectively.Unlike the single-jointed case, some modes of this structure oscillate in a manner that contains no exertion (i.e., relative displacement) of the joints.Consequently, the system behaves fully linearly around the associated resonances.Since the linear bars do not include dissipation, these same modes are effectively undamped.Thus we can expect unbounded response amplitudes near their eigenfrequencies.Among the plotted modes, the second and sixth modes behave in this manner.For the nonlinear system, Fig. 10 presents the response near the first mode in a format similar to Fig. 6.The response at a point immediately to the right of the first joint is used for all frequency response plots.As before, the unstable branches are plotted using dashed lines; however, for both WBVA and PWE solutions, we use the PER approach to assess their stability due to its ease of application.We note that it is still possible to use the strained parameter approach, expanding only one of the two linear stiffnesses, but the implementation complexity increases with the size of the problem.We again assess stability of the finite element HB results directly in the frequency-domain [45].The amplitude of harmonic excitation is set to 15 MN, where the peak nonlinear-to-linear restoring force ratio is 0.046.For the WBVA solutions, a clear convergence trend is observed with increasing orders of ε in the first harmonic amplitude plot (Fig. 10a).We note that the stability predictions for both wave-based approaches are again in very good agreement with the FE-HB predictions.
Fig. 11 plots the first harmonic responses of the system around the second and third resonances, forced with harmonic excitation amplitudes of 15 MN and 37.5 MN, respectively.As mentioned above, the second mode does not exert the joints appreciably.Consequently the nonlinear frequency responses are nearly unbounded at resonance.However, since a gradient-based implementation of numerical continuation is used to trace out the solution branches, the solver has a tendency to "jump" across the resonance when the peak becomes too steep.Since the local Jacobians as well as basins of attractions of the different approaches differ, in general, the exact point at which this happens differs, resulting in mismatches between the wave approaches in Fig. 11a.The exact nature of the resonance is more apparent in the corresponding phase diagram (Fig. 11c), which displays the discrete jump as a sharp phase change from 0 • to -180 • .The finite element approach predicts the peak at a slightly higher frequency (approximately 0.03%) than the wave-based approaches, similar to the trend in Fig. 7d.This mismatch is likely due to discretization error in which the finite element method is known to converge on frequency from above.Lastly, we note that the peak nonlinear-to-linear force amplitude ratios for the second mode are on the order of 1 × 10 −5 at both joints due to the absence of appreciable joint exertion.
Figs. 11b and 11d plot the third resonance peak where the peak nonlinear-to-linear restoring force ratios at both joints is approximately 0.075, comparable to that observed near the first mode.We also note similar trends to the first mode for convergence and stability.

Computational Performance
Apart from the theoretical interest in wave-based modeling for nonlinear structures, the reduction in computation cost is a remarkable feature of the wave-based techniques.In order to quantify this, we measure the total time taken to calculate forced responses at all the frequency ranges of interest, at multiple forcing amplitudes, using the wave-based approaches (WBVA and PWE) and the Finite Element Harmonic Balance (FE-HB) approach that we have employed for validation.
In order to avoid biasing the results due to poor performance of the specific implementation of   together with the same for the FE-HB approach, for both the single-and two-jointed examples studied in Sec. 4. All computations are conducted on an Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz (16GB RAM), with all computation run on a single core.For both examples, results in Table 2 document that the wave-based approaches provide an approximately 30× speedup in comparison to the FE-HB reference.Comparing WBVA to PWE, PWE can be seen to have marginally longer computation times.This is due to the fact that the WBVA uses efficient analytical expressions for scattering, leading to a minimal set of nonlinear algebraic equations, while the PWE must solve a full set of algebraic expressions.

Concluding Remarks
The paper develops multiple wave-based approaches for efficiently predicting periodic response in nonlinear, jointed structures.In doing so, we develop a new wave-based vibration approach for studying nonlinear systems informed by re-usable, perturbation-derived scattering functions.We also present a numerical formulation referred to herein as the plane wave expansion approach.We demonstrate that both wave-based approaches are highly-accurate and highly-efficient through comparisons with a more traditional approach (finite element harmonic balance).The WBVA has the most compact formulation and highest efficiency, but requires considerable up-front effort to derive the requisite scattering functions, while the PWE can be easily applied as new joint models are encountered and maintains solution accuracy to higher levels of nonlinearity than the WBVA.We remark that wave-based analysis of periodic and multi-physics vibration problems has received scant attention in the structural dynamics community, presenting further opportunities for follow-on research.One motivating example is that of a periodic bar or beam incorporating shunted piezoelectric elements, which can be configured to exhibit bandgaps, non-reciprocity, and/or topological states.
Figure 12: Ordered vibration sub-problems appearing following straight-forward perturbation applied to the forced vibration problem in Fig. 3.
missing on the right-hand side of Eq. ( 8) and (ii) inclusion of fixed boundary conditions; however, for the argument made herein, these equations are not strictly necessary and instead we only require inspection of Fig. 12.At O ε 0 appears a linear sub-problem amenable to the linear wave-based vibration approach.It is clear that the steady-state solution of such a problem yields unique displacement fields u 0 L (x, t) and u 0 R (x, t) responding at the excitation frequency ω.Following solution of the zeroth-order problem, substitution of u 0 L (x, t) and u 0 R (x, t) into Eqs.( 14)-( 15) yields the O ε 1 sub-problem appearing in Fig. 12.The boundary conditions at the joint for u 1 L (x, t) and u 1 R (x, t) now appear forced by terms with frequency content ω and 3ω resulting from the boundary term Γ (u R0 −u L0 ) 3 , which we indicate by F 1 (ω, 3ω).Once again it is clear that this linear subproblem yields unique displacement fields u 1 L (x, t) and u 1 R (x, t), now responding at the excitation frequency ω and its third harmonic.
This procedure can be carried-out to higher orders, but stopping at O ε 1 is sufficient to observe that, following solution reconstitution, the procedure yields a unique problem solution for each frequency ω, without the possibility of recovering Duffing-like frequency response curves and multiple solutions.We note that such a solution approach isn't without utility, particularly away from resonance, as it provides higher-order corrections to the fundamental frequency response and estimates of the higher harmonics otherwise missing from a linear analysis, but it is clearly inaccurate near resonance.We note further that this inaccuracy would not be encountered in a linear problem where solution uniqueness is guaranteed, and in fact, a similar approach is carried-out in the companion paper [41] where we consider stability of linear, parametrically-forced vibration problems.

Appendix C
Fig. 13 compares the harmonic balance (HB) results with time-transient simulations conducted on the Finite Element (FE) model of the single-jointed structure (see Sec. 4.1).An implicit Newmark scheme is used for the transient analysis.We chose the excitation frequency for the simulation as the frequency closest to the peak of the forced response in Fig. 7 ( 80.7079 krad/s), and the excitation amplitude is fixed at 15 M N , as before.We conduct the simulation for forty cycles of the fundamental period (40 × 2π/80.7079× 10 −3 s) in order to bring out any sub-harmonic features (which are expected at 1/3, 1/5, . . . of the excitation frequency) up to the 1/40 th sub-harmonic.The acceleration of the output node (directly after the joint, as before) is plotted since the higher harmonic effects are more pronounced in the acceleration than in the displacement (i.e., each frequency component gets scaled by the square of the frequency).It is readily observed that in the time-and state-space domains, the harmonic balance scheme that was employed (truncating to five harmonics) provides a reasonably accurate representation of the transient results, as expected.Another major aspect that is confirmed from this plot is that the true response is also periodic, thereby validating the choice of frequency-domain simulation techniques.
In the frequency-domain, the match is reasonable up to the fifth harmonic.The transient results show exhibit no prominent sub-harmonic behavior for this case (the 1/3 rd line is highlighted above).Further, intermediate (between H3 and H4, for instance) as well as higher harmonic components (H7 in the figure, for example) can be seen to exist, although these are not explicitly present in the HB formalism.The accuracy of the HB solution in spite of these components is due to the fact that their presence is very small in comparison to the other components.One way of quantifying this is through the ratio of the harmonic component with the first harmonic component.This turns out to be 3.89 × 10 −4 for the acceleration from the simulation above (2.51 × 10 −5 and 5.66 × 10 −5 , respectively, for displacement and velocity).It may also be observed that even harmonics do not appear on the plot since their magnitudes are near machine-precision.
The same argument is extended to justify the sufficiency of considering perturbation solutions only up to second order for the examples in Sec. 4. It can be seen from the wave-based scattering relationships developed in Sec.2.3 that the zeroth-order expansions are fully expressed in terms of a single harmonic alone, first order expansions include the third harmonic, and second-order expansions include the fifth harmonic.Since the above justifies the truncation of the solution up to five harmonics, it implies that comparable accuracy can be achieved through a second-order perturbation approach.In general, a convergence analysis has to be conducted to ascertain the appropriate order of perturbation by comparing the difference of the responses between predictions from the chosen order and one higher.

Figure 1 :
Figure 1: Two semi-infinite linearly-elastic bars joined by a nonlinear spring and a linear dashpot.Single-frequency plane waves are considered incident from the left-hand side.The incident wave is reflected and transmitted due to the presence of the nonlinear joint.
reflection and transmission of superharmonic content.Substitution of Eqs. (

Figure 2 :
Figure 2: Two semi-infinite linearly-elastic bars joined by a cubic spring and a linear dashpot.Multiple ordered plane waves are considered incident from both the left and right sides.

Figure 3 :
Figure 3: Two linearly elastic bars joined by a nonlinear spring and a linear dashpot.Each bar has one fixed end.Single-frequency 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.

Figure 4 :
Figure 4: A schematic representation of the evaluation of the joint nonlinearity in the Plane Wave Expansion (PWE) approach.

Figure 5 :
Figure 5: (a) Characteristic polynomial and (b-g) the first six mode shapes of the single-jointed system.The resonant frequencies are emphasised with markers in (a).

Figure 6 :Figure 7 :
Figure 6: Nonlinear forced response of the single-jointed system around the first resonance, undergoing a 7.5 MN amplitude periodic excitation: (a) first harmonic amplitude; (b) third harmonic amplitude; (c) first harmonic phase; and (d) fifth harmonic amplitude.

Figure 8 :
Figure 8: Three linearly elastic bars connected by two damped, nonlinear joints.

Figure 9 :
Figure 9: (a) Characteristic polynomial and (b-g) the first six mode shapes of the two-jointed system.The resonant frequencies are emphasised with markers in (a).

Figure 10 :
Figure 10: Nonlinear forced response of the two-jointed system around the first resonance, undergoing a 15 MN amplitude periodic excitation: (a) first harmonic amplitude; (b) third harmonic amplitude; (c) first harmonic phase; and (d) fifth harmonic amplitude.

Figure 11 :
Figure 11: Nonlinear forced response of the two-jointed system around the second (a,c) and third (b,d) resonances: (a,b) first harmonic amplitude; (c,d) first harmonic phase.The harmonic excitation for the second and third mode regimes are 15 MN and 37.5 MN, respectively.

Figure 13 :
Figure 13: Comparison between time-transient and frequency-domain simulations for the singlejointed structure (Sec.4.1) undergoing a harmonic excitation of amplitude 15 M N and frequency 80.7089 krad/s.Depicted clockwise are the time-domain acceleration response, the displacementvelocity state-space response, and the acceleration frequency-domain response.Extra labels are inserted in the x axis of the frequency domain showing relevant harmonics of the forcing frequency.

Table 2 :
Average computation times per point on the frequency response curve for (a) the singlejointed model and (b) the two-jointed model.Values scaled by the average response computation time for the FE-HB model in each case.numericalcontinuationemployed,wedivide the total time taken by the total number of steps in the frequency response curve.This provides an estimate of the average time taken to obtain the solution for a single point on the frequency response curve.In Table2we present the average computation times for the presented wave-based approaches