Efficient energy balancing across multiple harmonics of nonlinear normal modes

Predicting the forced responses of nonlinear systems is a topic that attracts extensive studies. The energy balancing method considers the net energy transfer in and out of the system over one period and establishes connections between forced responses and nonlinear normal modes (NNMs). In this paper, we consider the energy balancing across multiple harmonics of NNMs for predicting forced resonances. This technique is constructed by combining the energy balancing mechanism with restrictions (established via excitation scenarios) on external forcing and harmonic phase-shifts; a semi-analytical framework is derived to achieve both accurate/robust results and efficient computations. With known inputs from NNM solutions, the required forcing amplitudes to reach NNMs at resonances, along with their discrepancy, i.e. the harmonic phase-shifts, are computed via a one-step scheme. Several examples are presented for different excitation scenarios to demonstrate the applicability of this method and to show its capability in accurately predicting the existence of an isola where multiple harmonics play a significant part in the response.


Introduction
Close to the boundaries of their performance envelope, many engineering systems can experience nonlinear phenomena [7,23,42]. For example, an aircraft may be designed to be highly flexible to gain aerodynamic efficiency which, in turn, induces significant nonlinear behaviours during operation [36]; micro-electromechanical systems can be employed as innovative bandpass filters with improved performance by exploiting the modal coupling due to nonlinearity [8].
To understand these nonlinear behaviours, the extension of linear modal analysis to account for nonlinearity has been extensively considered. An early approach was proposed by Rosenberg [31,32], where a nonlinear normal mode (NNM) is defined as a synchronous periodic response of a undamped and unforced, or conservative, system. This definition was later relaxed to a periodic response (not necessarily synchronous) of a conservative system [19]. Alternative definitions of a nonlinear mode include a damped NNM [20,35,39], a spectral submanifold (SSM) [9], and a phase resonance nonlinear mode [41]. Using the concept of nonlinear modes, many nonlinear behaviours have been investigated, e.g. modal interactions, bifurcations and instability [1,[14][15][16]21,40], to facilitate the understanding of their mechanisms and aid practical design. In addition, owing to their invariance property, analysis of large-scale nonlinear systems also becomes feasible via model order reduction techniques [9,18,24,40].
In practice, engineering systems are usually operating under forced conditions, where the forced response curves (FRCs) are widely considered for analysis and design of nonlinear systems [19,42]. Employing an energy-based framework, the forced responses can be used to extract NNMs for modal analysis [25] or vice versa-using NNMs to interpret the forced responses [13].
In [25], a nonlinear extension of the force appropriation technique was developed in order to identify NNMs in experimental setups. To this end, a set of multi-point multi-harmonic excitations are tuned to exhibit phase quadrature with all harmonics of displacements, in order for the excitations to compensate the damping effect. Successful applications of the force appropriation technique were demonstrated in a large number of studies in combinations with, for example, resonant decay [26], control-based continuation [28,29] and phase-locked loop [5,27]. These studies show that, for most cases, even a single excitation is sufficient to accurately locate the NNM branches. However, numerical and experimental studies also show limitations and significant errors, arising from the insufficient compensation between the single-location forcing and distributed damping of the structure [30]. Alternative approaches, e.g. responsecontrolled stepped-sine testing [17] and velocity feedback [33], can as well be utilised for nonlinear modal testing.
For nonlinear systems operating under forced conditions, the resonances are of particular interest as they usually represent the most significant responses in a system. Near resonances, the forced responses may be seen as perturbations from NNMs [2,13]. Using these observations, an energy-balancing method was used to establish the relationships between forced responses and NNM branches [13]. This is achieved by assuming resonances as phase-shifts of NNMs and considering the energy balancing during periodic responses-the net energy transfer in and out of any mode must be zero. The energy balancing method has also been used to predict the presence of isolas [6,10,21]-i.e. regions of forced responses that are separated from the primary branches of the response; to quantify the relative importance of NNMs by predicting the forcing amplitude required to reach them at resonance [22]; and to identify an appropriate number of excitations and their locations in conducting force appropriation tests [30]. This energy balancing method was origi-nally derived based on systems with conservative nonlinearity and, recently, it has been extended to account for non-conservative nonlinearity with applications in full-scale structures [37] using the concept of damped NNMs [20]. This extended energy balancing method has been compared with nonlinear modal synthesis in predicting resonances [38].
To complicate matters, not all NNMs are pertinent to forced responses-only those NNMs whose fundamental components are coupled with specific phase relationships correspond to forced responses when forcing and damping are applied [11,12]. In [2], when a conservative NNM may be perturbed to forced responses are quantified via Melnikov analysis. Another limitation of the energy balancing technique, proposed in [13], lies in the assumption that a single harmonic is representative of the modal response. This single-harmonic assumption may fail to capture resonances where multiple harmonics show significant contributions. Such a limitation was demonstrated in [6], where significant errors, partly arising from the single harmonic assumption, can be observed during the emergence of an isola to the primary forced response curve.
Being a versatile technique with applications in both theoretical and experimental studies, one of the main advantages of the energy balancing method is its analytical framework, which in turn can bring inaccuracy in complex application scenarios as it is typically assumed that only the fundamental harmonic is significant. However, the analytical framework does not allow for a direct extension to account for multiple harmonics; instead, computationally expensive numerical schemes are required. In this study, we present an extended technique that overcomes the single-harmonic limitations, whilst preserving the computational efficiency by employing a semi-analytical scheme. To this end, the rest of the paper is organised as follows.
Firstly in Sect. 2, the mechanism of energy balancing during periodic responses, proposed previously in the literature, is reviewed highlighting the relationships between forced responses and NNMs. A nonlinear beam system is considered to show the applicability and limitations of single-harmonic energy balancing method. The single-harmonic energy balancing method is then extended to account for multiple harmonics in Sect. 3, when quadrature excitations are applied to the modal harmonics. A semi-analytical framework is obtained by combining energy balancing across har-monics and forcing/phase-shift constraints. In Sect. 4, this technique is further extended to more practical scenarios where the excitations are in quadrature with physical displacements. Throughout the paper, the application of the proposed method is compared with the force appropriation technique [25] and singleharmonic energy balancing method [13] via a series of examples.
2 Energy-transfer balancing in nonlinear systems: from system to harmonic levels This section first reviews the mechanism of energy balancing of nonlinear systems during periodic responses. This is achieved by demonstrating the energy-transfer balancing at the level of the system and of the mode 1 , and then to that of the harmonic. From such an energybased perspective, the relationships between nonlinear normal modes (NNMs) and forced responses are highlighted. Discussions in this section provide a basic concept which we build on in this paper, where an efficient technique for applying harmonic-level energy-transfer analysis (HETA) is proposed-as introduced in Sect. 3 and Sect. 4.
2.1 Energy-transfer balancing at the system level For a nonlinear system, its dynamics can be expressed, using its modal components, in the form whereq,q and q are vectors of modal accelerations, velocities and displacements, respectively. The diagonal matrix D contains the linear damping coefficients. This reflects a significant number of engineering systems with geometric nonlinearity where the damping may be assumed small when compared to external forcing and nonlinear stiffness. As such, the simple case of linear modal damping is often assumed [5,6,25,26,29]. A discussion on nonlinear damping is given in [37]. is a diagonal matrix containing the squares of the linear natural frequencies. Vector N q (q) contains the nonlinear stiffness terms (assumed to be conservative, and a function of the displacements, q). The external periodic forcing is captured by the vector p.
At the system level, energy-transfer balancing considers the net energy transfer in and out of a system over one period of motion. As the nonlinear terms considered here are conservative, this energy transfer can only be achieved via the damping and external forcing terms. The energy-transfer terms arising from damping and external forcing over a period are termed damping and forcing energy-transfer terms, denoted E D and E P , respectively. For periodic motions, the total energy loss by the system due to these nonconservative terms must sum to zero, such that In the modal domain, the nonlinear system, described by Eq. (1), can be seen as a collection of modes, where the equation of motion of the ith mode is written and where d i , ω ni and p i are the ith modal damping coefficient, linear natural frequency and external modal forcing, respectively. As such, the energytransfer terms in Eq. (2) can be translated, using modal components, into summation form where E Di and E Pi are the ith modal damping and forcing energy transfer, respectively, and N is the total number of modes in the system. Note that the nonlinearities, described by N qi (q), may lead to coupling and energy transfer between modes. Likewise, the nonlinearityinduced net energy transfer to the i th mode from other modes in an NNM over a period is termed nonlinear energy transfer and denoted E Ni . The effect of E Ni can be seen as internal energy rearrangement among modes; hence, it exhibits a conservative effect at the system level and has zero contribution to the energy balancing in Eqs. (2) and (4).

Energy-transfer balancing at the mode level
The concept of energy balancing analysis can be extended to the mode level by considering the energy balancing for each modal component. At this level, in addition to damping and forcing which both allow energy transfer in and out of modes, nonlinearity can also lead to energy transfer between modes, which must (a) -On an NNM Fig. 1 A schematic of the energy transfer at the mode level, for a two-mode system. Panel a shows the case where the forcing and damping match precisely, so that the response may equal an NNM. In panel b, the forcing and damping do not match and the response is not precisely on the NNM also be accounted for. To compute these energy transfer terms, the ith modal equation of motion (3) is firstly multiplied by the velocity of the ith mode,q i , and then the time integral over one period is considered. The obtained expression is given by where the terms arising from damping, external forcing and nonlinearity are defined as and where the terms related to acceleration and linear stiffness are zero due to orthogonality. Expression (5) can be interpreted as an energy-based mapping of periodic motions governed by Eq. (3), and it describes the energy balancing at the mode level-the net energy transfer in and out of any mode over one period must be zero, similar to Eq. (2) for that at the system level. This energy balancing relationship can also be intuitively obtained from the two-mode schematic shown in Fig. 1b, where the energy transfer arising from damping, external forcing and nonlinearity must be balanced for a steady-state response.
At the mode level, for an NNM (i.e. undamped and unforced periodic response), the net energy transfer between modes must be zero, i.e. E Ni = 0, as no mode may lose or gain energy over one period and still remain periodic. Therefore, if an NNM solution is precisely equal to a forced response, Eq. (5) implies that the forcing energy gain must be precisely equal to the damping energy loss for each mode, i.e.
A schematic of this energy transfer is shown in Fig. 1a, for a two-mode system. In most applications, however, the applied forcing cannot precisely satisfy Eq. (7), i.e. the forcing energy gain for each mode does not equal the damping energy loss (for example, a mode may be damped but unforced). Such a forcing case is referred to here as an imperfect forcing. Therefore, Eq. (5) reveals that, for an imperfect forcing, E Ni = 0, which violates a condition of an NNM (i.e. that E Ni = 0). A schematic of this energy transfer is shown in Fig. 1b. This implies that an NNM can never be precisely reached by a system with imperfect forcing. To quantify the deviation in the energy transfer from Fig. 1a to Fig. 1b, caused by imperfect forcing, the effect of nonzero nonlinear energy-transfer terms, E Ni , may be considered.
In Ref. [2,13], the forced response in the neighbourhood of an NNM, where E Ni = 0, is considered as a perturbation from this NNM. In estimating the nonzero E Ni , the only parameters that can quantify a change in the periodic response are response frequency, phase and amplitude. For a given periodic forcing, the response frequency must be fixed to remain periodic; the change of phase can directly lead to an energy exchange between modal components, whilst the change of amplitude will not contribute to an energy exchange unless accompanied by a phase change. As such, the amplitude change may be seen as a second-order effect compared to the phase change. Consequently, we consider the phase change, termed the phase-shift, of an NNM that is necessary to accommodate the internal energy transfer, E Ni , as in [13]. If this phase-shift is sufficiently small, it may be assumed that the forced response is close to that of an NNM; conversely, a large phase-shift corresponds to a large change in the response, indicating that the forced response is significantly different to the NNM. Note that, in [13], it assumes the dynamics of a mode may be approximated by a single harmonic and is established Fig. 2 A schematic of a cantilever beam with a nonlinear spring at the free end using Eqs. (5) and (7); hence, it is termed the modelevel energy-transfer analysis (META) in this paper.
To demonstrate the application of META, a motivating example is considered, shown in Fig. 2-a cantilever beam with a cubic nonlinear spring at the free end, and excited at a point on the part-span position, L 1 , by a single-harmonic force The deflection of this beam is modelled using the first two modes, and the modal equations of motion are given bÿ where θ i is the modeshape of the ith mode at the excitation point, and where d i = 2ω ni ζ where ζ is the modal damping. The nonlinear terms are given by where α i denote nonlinear modal coefficients. The linear stiffness of the spring is tuned such that the two modes exhibit a 1:3 response, i.e. the fundamental (largest) harmonic of the second mode responds at three times the frequency of the fundamental harmonic of the first mode. As this system (which is analogous to that considered in [34]) is used as an illustrative example, a derivation of the expressions is not provided here; however, interested readers are directed to Refs. [34,42] for an overview of the approach taken for this derivation.
Here, an example system with parameters given in Table 1 is considered. The first NNM branch (i.e. the locus of NNMs emerging from the first linear natural frequency) of the beam is found using the numerical continuation software COCO [4]. It is shown as a blue line in Fig. 3a(ii) in the projection of the response frequency, ω, against the maximum response amplitude of the second mode, Q 2 .
A periodic force, f (t), is applied at a part-span position, L 1 = 0.7L, to consider forced responses of the beam. The relationship between the NNM branch and the forced responses is identified via META using Eq. (7) and shown in Fig. 3a(i)-the forcing amplitude, F 1 , required for the forced resonant responses to share solutions with the NNM branch at response frequency, ω. Here, two forcing cases are considered and they are denoted by dashed lines in panel a(i). The intersections, labelled by a square and three dots, respectively, indicate predicted resonances on forced response curves. To verify these predictions, the corresponding forced responses are found using numerical continuation and shown as green and red lines, respectively, in panel a(ii). Excellent agreements are achieved for both cases-the predictions well capture the existence of the isolated forced responses and the development of resonances.
The excitation location is then moved from L 1 = 0.7L to L 1 = 0.2L, and the energy balancing relationship is again computed via META, shown in Fig. 3b(i). Here, we again consider two forcing cases, denoted using dashed lines. As in Fig. 3a(i), the lower-amplitude forcing case predicts one intersection, expected to be related to one primary forced response curve, whilst the higher-amplitude forcing case predicts three intersections, corresponding to one primary curve and one isolated curve. To verify these predictions, the forced responses are computed numerically and shown in panel b(ii). Results obtained via META capture the resonances on the primary forced responses-the predictions are in line with the primary resonant peaks for both cases. In contrast to the cases shown in Fig. 3a, the META is unable to predict the existence of the isolated curve for the lower-amplitude forcing case; similarly, the resonances are not well captured for the higheramplitude forcing case-see the discrepancies between dots and resonances on the isola.
To explain why META succeeds for the cases where the forcing location is at L 1 = 0.7L, but fails when L 1 = 0.2L, the energy-based framework is retained to aid the discussion, but the energy balancing at the mode level, described by Eq. (5), is further extended to consider the periodic energy transfer between harmonics. From such an energy-based perspective, we seek to explain the limitations of META and show the need for further extension, i.e. the contributions of this work.
As with the interpretation that a system may be seen as a collection of modes, a mode may be seen as a collection of harmonics. This allows the modal coordi- Table 1 The parameters of the two-mode beam model, described by Eqs. (9) and (10) Parameter  nates to be expressed as a sum of harmonics where u i, j is the jth harmonic of the ith mode and H i is a set denoting the harmonics in the ith mode. Combining Eqs. (11) and (5), the energy transfer at the mode level can be expressed in summation form where the terms within the summation represent the net energy transfer to the jth harmonic of the ith mode arising from damping, forcing and nonlinear terms, respectively.
Here, to analyse the limitations of META, the periodic energy balancing at phase resonances, i.e. forceddamped motions in quadrature with the forcing, on the isolated branch are considered. Two example cases, labelled '×' in the embedded plots of Fig. 3a, b, are con-sidered, and they are related to instances where META succeeds and fails in predicting the existence of the isolated branch. Using Eqs. (6), (11) and (12), the periodic energy-transfer terms may be computed for each harmonic, and they are shown in Figs. 4a and 4b with respect to labelled cases in embedded plots of Fig. 3a, b.
Applying an excitation with forcing frequency, , near the first natural frequency gives energy to the fundamental component of the first mode, u 1,1 , and the first harmonic of the second mode, u 2,1 . When the forcing is at L 1 = 0.7L, Fig. 4a shows that the major energy input is directed to u 1,1 , with only one percent of energy imparted to u 2,1 . Whilst the major internal energy transfer, triggered by nonlinear coupling terms, occurs between fundamental components of the two modes, u 1,1 and u 2,3 , where the transferred energy is lost via damping. Note that the response being considered exhibits a 1 : 3 modal interaction, so u 2,3 is the 0.987 fundamental component of the 2nd mode response. As such, ignoring the small contributions from harmonics, the energy balancing, Eq. (12), can be approximately captured by the fundamental components. This also gives rise to the accurate predictions via META (see Fig. 3a), where the responses of the fundamental components are used to represent the two modes. For the case in Fig. 4b where the excitation is at L 1 = 0.2L, the major energy input from external forcing is still imparted to u 1,1 , as with the case in Fig. 4a; however, with an increased amount given to u 2,1 (about nine percent of the total energy input). As for the internal energy transfer between harmonics, E Ni, j , besides the fundamental components, a significant involvement of u 2,1 can also be observed. As such, without considering the energy transfer arising from harmonics (u 2,1 for this case), a significant error in energy balancing analysis via META is shown in Fig. 3b, which is unable to capture the existence of isolated forced responses for the lower-amplitude forcing case.
By comparing cases in Figs. 3, 4, the limitations of META can be explained from an energy-based perspective by considering the energy transfer between harmonics. It also demonstrates the importance of considering multiple harmonics in energy-transfer analysisthe topic of this paper. To establish energy balancing analysis considering multiple harmonics, the modelevel energy balancing, i.e. Eqs. (5) and (12), is extended to the harmonic level.

Energy-transfer balancing at the harmonic level
As with the extension of energy balancing analysis from the system level to the mode level, it can be further extended to the harmonic level. Multiplying the ith modal equation of motion (3) with the velocity of the jth harmonic of the ith mode,u i, j , and integrating over a period of time, the energy balancing at the harmonic level can be obtained where the energy-transfer terms are given by Here, we note that the integrals of acceleration and linear stiffness are zero. Expression (13) describes the balancing of net energy transfer in and out of any harmonic, and this may be seen as the fundamental level as it is established using the fundamental elementsharmonics, and describes the fundamental mechanism of energy transfer of nonlinear systems during periodic responses. At the harmonic level, for an NNM, the net energy transfer between harmonics must be zero, i.e. E Ni, j = 0. Therefore, the relationship between NNMs and

Fig. 5
A schematic of the energy transfer between the harmonics of a two-mode system forced responses can be obtained from Eq. (13), which reads This indicates that the harmonic forcing energy gain is balanced by the harmonic damping energy lossan extension of Eq. (7) to the harmonic level. Using this relationship, the force appropriation technique has been proposed to extract NNM branches experimentally [25]. In order to satisfy Eq. (15), a quadrature criterion, where the excitation has to compensate for the damping terms, is achieved by applying 90 • phaselagged forces that contain all harmonic components of all modes in an NNM. However, in practice, such a perfect forcing set can never be precisely achieved, which causes the response to deviate from the NNM-a nonzero energy transfer to harmonics, i.e. E Ni, j = 0. Such a nonzero E Ni, j may couple all harmonics of all modes within an NNM, and this will be derived analytically in the following discussions. A general energy-transfer schematic for a two-mode multi-harmonic system is shown in Fig. 5, which represents an extension of that shown in Fig. 1b; and is a generalised version of that used for the example we have just discussed (Fig. 4). This illustrates a more complex energy-transfer network among all harmonics at the harmonic level, when compared to taking a higher-level view at the mode level.
Similar to the mode level, to evaluate the deviation of forced responses (E Ni, j = 0) from NNM solutions (E Ni, j = 0), the effect of nonzero E Ni, j may be considered. As with the META, where a modal phaseshift is used to evaluate a nonzero E Ni , the harmonic phase-shift is used to evaluate the nonzero E Ni, j for the energy balancing at the harmonic level. This allows the harmonic-level energy-transfer analysis (HETA) to be analytically formulated, detailed in the following.

Harmonic-level energy-transfer analysis for systems under harmonic forcing
This section considers the extension of the modelevel energy-transfer analysis, proposed in Ref. [13], to consider energy transfer between harmonics. This is achieved by firstly interpreting the forced-damped near-resonant responses as perturbations from the NNM solutions due to harmonic phase-shifts. The energy-transfer analysis over all harmonics of all modes is then formulated using Eq. (13), which leads to an underdetermined problem, i.e. the unknowns outnumber the equations. To solve this, the force reduction, along with extra phase constraints imposed by the quadrature conditions, are introduced to formulate a determined, or solvable, equation set.

Harmonic phase-shifts of NNMs under quadrature forcing
To account for multiple harmonics, the displacement of the ith mode, for an NNM response, is approximated as a sum of a finite number of harmonics, such that where the overbar,•, indicates that this is an NNM response, and whereū i, j is the jth harmonic of the ith modal displacement,q i ; and H i is the set of harmonics used to approximateq i . The amplitude and phase ofū i, j are represented byŪ i, j andφ i, j , respectively, and the frequency ω is defined as ω = 2π T −1 , where T is the period of the response of the system (i.e. considering all modes). As with the modal displacements, the external forcing applied to the ith mode may be separated into harmonic components, written (17) Note that this assumes that the harmonics of the excitation force are all in quadrature (i.e. at π/2 out-ofphase, or π/2 phase-lagged) with the corresponding modal harmonic. This also assumes that all harmonics are forced; however, the unforced harmonics may be specified as those where the excitation amplitude, P i, j , is zero-this is revisited later in Sect. 3.3.
Following the approach used in [13], as described in Sect. 2, the forced-damped near-resonant response may be seen as a perturbation from an NNM response with a shift in the phase of the response. This can also be justified using the Melnikov analysis in [2], where the forced responses are proven to be perturbations from NNMs when forcing and damping are small in comparison to conservative terms. Here, this phase-shift is applied to all harmonics of the response, rather than just the fundamental component (as considered in [13]). The phase-shift of the jth harmonic of the ith mode is writtenφ i, j , and it is assumed that the phase-shift is small in viewing forced responses as perturbations from NNMs. Introducing these phase-shifts to Eq. (16), the forced response of the ith mode is written Note that, as in [13] and discussed in Sect. 2.2, the change of amplitude may be seen as a second-order effect when compared to the phase change in energytransfer contributions, it is assumed that all amplitudes in the forced response are equal to those in the NNM (U i, j =Ū i, j ), i.e. the amplitudes are unaffected by the application of forcing and damping. Equation (18) also indicates that, provided the phase-shift is a small term, an approximation of the forced response can be achieved by accommodating the NNM solution with this phase-shift.
Based on this smallness assumption of phase-shift, we now formulate the energy-transfer analysis at the harmonic level via Eq. (13), where the energy-transfer terms, E Di, j , E Pi, j and E Ni, j , may be approximated as linear functions of the excitation amplitude, P i, j , and the phase-shifts,φ i, j .
3.2 Problem formulation: energy-transfer balancing at the harmonic level Equation (13) shows the balancing of the damping, forcing and nonlinear energy-transfer terms for the jth harmonic of the ith mode. These terms may be computed, respectively, using Eq. (14). In addition, in Appendix A, it is shown that these terms may be written as where the dagger, • † , denotes a known term-i.e. a term that may be computed using the NNM response (which is assumed to be known), whilst the forcing amplitudes, P i, j , and phase-shifts,φ n,k , are assumed to be unknown. The known terms are found using Equation (19b) shows that the energy transferred to the jth harmonic of the ith mode from external forcing is only due to the force directly applied to that harmonic (i.e. it is only a function of P i, j ), of which examples are shown in Fig. 4. The nonlinear counterpart, Eq. (19c), represents the nonlinear energy transfer from all harmonics of all modes in an NNM to the jth harmonic of the ith mode; and it reveals that a phase-shift in any harmonic of any mode, i.e. the kth harmonic of the nth mode, may lead to an energy transfer to the jth harmonic of the ith mode. This demonstrates the feature of nonlinearity-induced energy transfer in coupling all harmonics of all modes in an NNM, as shown in Fig. 5.
Substituting expressions (19) into the energy balancing expression, Eq. (13), gives By collecting the energy-transfer terms due to damping, forcing and nonlinearity, the energy balancing can be formulated from the jth harmonic, to all harmonics of the ith mode, and to all modes of the system, summarised in Table 2, where the total number of harmonics used to represent all N modes of the system is given by Energy balancing for q i Expressions for coefficients, {matrix dimension} Energy balancing for the system Expressions for coefficients, {matrix dimension} Up to this point, the energy-transfer analysis is formulated as a set of R H equations, given by Eq. (27), with 2R H unknowns of phase-shifts,φ i, j , and forcing amplitudes, P i, j . To solve this underdetermined equation set, where the unknowns outnumber the equations, either reducing the number of unknowns, introducing extra constraints, or both, are needed.

Problem solving: quadrature harmonic forcing
As previously discussed, the current formulation allows all harmonics of all modes to be forced. However, in practical applications, the forcing can only be applied to a limited number of harmonics. Here, a special case is considered where a limited number of quadrature harmonic forces 2 are applied, whilst the more complex case, where quadrature physical forces are applied, will be discussed later in Sect. 4. To specify this here the harmonics of the nonzero forcing (applied to the ith mode) are defined as belonging to the setP i , i.e. the nonzero subset of P i . As such, the total number of nonzero forcing applied to the system is given by With this, P may be simplified by discarding the zerovalued elements to reduce the number of unknowns that need to be estimated. This reduction is achieved by relating the vector of nonzero forcing amplitudes P (measuring {R F ×1}) to P via the force reduction matrix, C P , where P = C PP .
The force reduction matrix may be constructed using a {R H × R H } identity matrix and removing the columns associated with unforced harmonics. Hence, whereẼ † P = E † P C P is also used, and henceẼ † P measures {R H × R F }. Through force reduction, the total number of unknowns is now reduced to (R H + R F ), i.e. R H phase-shifts and R F forcing amplitudes. However, it is still underdetermined and additional R F constraints are required to solve the equations.
As it is assumed that all forces are in quadrature with the harmonic they are forcing, it may also be assumed that the forced harmonics do not exhibit a phase-shift 3 . This condition may be enforced using a phase constraint matrix, C φ , where Here, C φ is a {R F × R H } matrix that constrains the phase-shifts of all forced harmonics to be zero. In the case where all harmonic forces are independent, the phase constraint matrix is the transpose of the force reduction matrix, i.e. C φ = C P ; hence, C φ may be constructed by removing the rows associated with unforced harmonic of a  (31) and (32) are now combined and written into matrix form as where and where the unknowns, i.e.φ andP, are collected in the vector v, and the known coefficients are collected in matrices A and B. As the number of constraints now matches the number of unknowns, the unknown terms are found using This allows the phase-shifts of any unforced harmonics and the amplitudes of any forces, which are represented in v, to be computed. As discussed above, the NNM solutions are viewed as known inputs that can be computed via numerical simulations. With these known parameters, the proposed method can be formulated by combining energy balancing expressions with force reductions and phase constraints, and closed-form solutions can be obtained analytically via Eq. (35). Note that this analysis usually results in large matrices, so it is expected that whilst this can all be defined analytically, it is more practical to compute digitally. It should also be noted that the derivations above are based on the smallness assumption of phase-shifts; relaxing this assumption leads to computationally expensive numerical simulations where such closed-form solutions (35) are no longer available.
For the case where the forces applied to the harmonics are independent, the use of the forcing reduction matrix, C P , and the phase constraint matrix, C φ , could be viewed as unnecessary:Ẽ † P may be constructed directly by removing the columns of E † P that correspond to the unforced harmonics and a lowerdimension problem could be formulated by simply removing the forced phase-shifts and their corresponding columns in E † N . However, this approach allows additional constraints to be introduced, as shown later in Sect. 4. The implementation of this method is now demonstrated using a simple example system.

Example 1: a quadrature harmonic forcing case
To demonstrate the HETA, the cantilever beam, schematically shown in Fig. 2, is again considered. To simplify this case, a single-harmonic force is applied to each mode independently, such that the equations of motion are written where the nonlinear forces are given by Eq. (10). Note that the case where forcing is applied to the part-span will be revisited later in Sect. 4. The parameters of the system are the same as considered in Sect. 2.3, given in Table 1. It should be noted that the forcing amplitudes, F 1 and F 2 , are assumed to be unknown (and to be found using HETA). Also note that the linear natural frequencies have a ratio of approximately 1 : 3, i.e. ω n2 /ω n1 ≈ 3. This leads to a 1 : 3 modal interaction, as discussed in detail in [34].
Before applying HETA, the NNM branches of the system need to be obtained to provide known parameters to compute the coefficients in Eq. (35). The first NNM branch is presented in Fig. 6 by a solid-blue line. The top and bottom panels of Fig. 6 show the maximum displacement amplitudes of the first and second modes, respectively. In both of these panels, a distinctive loop region is seen, indicating a strong internal resonance between the two linear modes. Further discussion of this NNM branch can be found in [34]. To aid comparison, the forced responses, obtained via quadrature forcing, or force appropriation technique [25], are also presented as thin red lines, with dots. Note that the NNMs to be identified are synchronous periodic motionsa required condition in applying force appropriation [25]. For the system described by Eq. (36), quadrature is achieved when the force applied to the first mode is at 90 • to the displacement of the first harmonic, u 1,1 , and when the second modal forcing is at 90 • to the displacement of the third harmonic, u 2,3 -see [25,30] for further details. As the NNMs and the forced responses appear to be close, it is expected that the phase-shifts will be small.
To apply HETA, a finite number of harmonics must be used to approximate the modal displacements, as previously discussed in Eq. (16). For this example, it is assumed that the first and third harmonics of each mode are sufficient 4 , i.e.
such that Whilst additional harmonics would provide greater accuracy, just two harmonics are considered here for simplicity. For this example, the construction of energy balancing analysis via Eq. (33) is not presented here, but instead, results and analysis are given directly to show the power of the proposed method. For details of the construction, readers are directed to Appendix B. Figure 7 shows the forcing amplitudes, |F 1 | and |F 2 |, and the phase-shifts,φ 1,3 andφ 2,1 , found via HETA outlined above and force appropriation. Note that, in applying force appropriation, the phase-shifts denote the phase differences between forced responses and NNM solutions, obtained from COCO solutions [4]. In all panels of Fig. 7, the computationally cheap analytically predicted values, obtained via HETA, are represented by blue lines, and the computationally expensive numerically simulated values, obtained via the force appropriation technique, are shown by red lines with dots. The difference between the analytically predicted and numerically simulated forcing amplitudes, |F 1 | and |F 2 |, is indistinguishable, showing that these have been predicted with a very high level of accuracy. Although there is some discrepancy, the phase-shift values,φ 1,3 andφ 2,1 , also show a very good agreement, despite the low number of harmonics used to approximate these responses. The greatest inaccuracy can be seen near the linear natural frequency, ω n1 , (where NNMs are close to the linear response) in theφ 1,3 andφ 2,1 phaseshifts. This is likely due to numerical error in both the simulated and estimated values, as the harmonics are very small for low-amplitude NNMs. Additionally, as NNMs approach the linear case, the mixedmode NNM response approaches a linear single-mode response whose phase is no longer defined.
Besides the strong agreement between results obtained via these two techniques, tracing the evolution of phase-shifts, shown in panels (b) and (c), can monitor and estimate the applicability and accuracy of energy balancing analysis in 1. using NNMs to interpret forced responses or vice versa, 2. using forced responses to identify NNMs, i.e. the force appropriation technique, as the phase-shifts capture the perturbations from NNMs to forced responses. In this case, the smallness of phase-shifts justifies the applicability of the proposed method. In addition, by accommodating the To summarise, the results of HETA, shown in blue in Fig. 7, have been constructed in Appendix B following the procedure outlined in this section using only the NNM responses. For the sake of simplicity, this approach has assumed that the response consists of just two harmonics, which limits the accuracy of the predictions. However, owing to the semianalytical framework, it is convenient to increase the number of harmonics without adding much computational effort. As such, this provides a more efficient technique, whilst preserving the necessary high resolution, for resonance predictions when compared with computationally expensive numerical analysis. The following section will consider additional, more complex, examples and demonstrate the application of this approach to cases with a greater number of parameters.

Accounting for quadrature physical forcing
The formulation of HETA presented in Sect. 3 is restricted to cases where independent harmonic forces are applied to each mode of the system; however, in many practical applications, a force may be applied to multiple modes simultaneously. For example, a force applied to the mid-span of a beam is able to excite all modes of the system simultaneously (aside from those that have a node at the forcing location). If one excitation force is applied to multiple modes, the forcing cannot be in quadrature with all forced modal responses, and hence the phase constraints used in Sect. 3 (i.e. the forced harmonics exhibit zero phase-shifts) become invalid; instead, the modes must be free to exhibit different phase-shifts. In this section, it is shown that HETA may be extended to account for additional forcing conditions-namely where a force is in quadrature with a displacement at a physical location, rather than a harmonic of a mode. This is reflective of most practical scenarios where a force is tuned to reach quadrature with the excitation point. Here, such a forcing is termed a quadrature physical forcing, in comparison to the quadrature harmonic forcing considered in Sect. 3.

Formulation of HETA considering quadrature physical forcing
To formulate HETA with quadrature physical forces, the excitation points are defined as ∈ L, whose number of elements is R L ; at the excitation point, , the harmonic components of the forcing are defined as j ∈ F , whose number of elements is R F . Note that it is assumed that the forced harmonics are included in all modes, i.e.
Consequently, the number of forcing harmonics, across all excitation points, is given by At the excitation point , the ith modal component of the external forcing,p ,i , can be obtained by introducing the linear modal transform to the physical forcing, and is given bỹ where θ ,i denotes the modeshape of the ith mode at the excitation point, . Separating into harmonic components,f can be expressed as whereF , j andψ , j are the amplitude and phase of the harmonic forcing,f , j . As such,p ,i may be approximated as a sum of harmonics using Eqs. (41) and (42), of which the jth harmonic component is given bỹ The energy transfer from external forcing,p ,i , to u i, j can be computed via Eq. (14b), given bỹ where The details of this derivation are given in Appendix A. The phase relationship between the external forcing and NNMs is captured by φ di, j . As such, the sign of φ di, j accounts for the effect of external forcing on energy transfer-either an energy gain or an energy loss, reflective of the fact that one force may affect modes differently. In contrast, for the case where independent harmonic forces are applied in quadrature with the harmonic displacements, as discussed in Sect. 3, the applied force must lead to a forcing energy-gain term for the enforced harmonic-see Eqs. (19b) and (20b). As shown in Appendix A, φ di, j may be related to NNM solutions via quadrature conditions, and hence, it is seen as a known parameter. As discussed in Sect. 3.2, the forcing energy-transfer terms can be separated into forcing energy coefficients, E † P ,i, j , and harmonic forcing amplitudes,P ,i, j , i.e. where (48) Combining the forcing amplitudes for all harmonics of the ith mode gives Next, collecting modal forcing amplitudes for all enforced modes of the system gives wherẽ P = P ,1P ,2 . . .P ,i . . .P ,N , (52a) As with discussions in Sect. 3, the reduced forcing energy coefficient matrix,Ẽ † P, , can be formulated viã where E † P is a diagonal matrix with leading diagonal elements, E † P ,i, j , defined in Eq. (48); and the force reduction matrix, C P , may be constructed using a {R H × R H } identity matrix and removing columns associated with unforced harmonics. The forcing energy-transfer terms are then collected with respect to known terms-parameters that may be computed via NNM solutions and unknown termsunknown forcing amplitudes to be solved, i.e.
whereẼ † F = E † P C P θ . For each excitation point,Ẽ † F andF may be constructed using the procedure outlined above. The direct collection of them across all excitation points gives Other energy-transfer terms that arise from damping and nonlinearity, E Di, j and E Ni, j , are the same as those defined in Sect. 3, and they can be constructed in the same form as E † D and E † N , respectively, via Eqs. (28a) and (28d). As such, the energy balancing expression with quadrature physical forcing is formulated as To construct phase-shift constraints, the displacement at the excitation point is written y and may be expressed as a sum of modal displacements, i.e.
Considering the modal displacement as a sum of harmonic components, as in Eq. (16), gives and hence the jth harmonic of the physical displacement y may be written Here, as y , j represents the harmonic displacement which is in quadrature with the forcing, it may be written as a sinusoid with amplitude Y , j and phaseψ , j , i.e.
where Eq. (18) has been used to express u i, j as a sinusoid and where φ i, j =φ i, j +φ i, j . If the jth harmonic of the force is in quadrature with y , j , then the phase,ψ , j , of y , j must be equal for both the NNM and the forced case. To find how this leads to a constraint between the harmonic phase-shifts, the timeand phase-dependent components of Eq. (60) may be separated by first writing From this, the cos ( jωt) and sin ( jωt) components may be balanced to give the relationships The displacement amplitude Y , j may now be removed by dividing Eq. (62b) by Eq. (62a) to give If the phase,ψ , j , is equal for both the NNM and forced case, it therefore follows that Eq. (63) may be satisfied for both φ i, j =φ i, j +φ i, j (i.e. forced responses) and φ i, j =φ i, j (i.e. NNMs). Assuming the phase-shifts are small, this may be approximated to This restriction leads to a constraint between the harmonic phase-shifts and may be simplified to In this case, each harmonic of each forcing location will be associated with an unknown forcing amplitude, but will also lead to an additional constraint given by Eq. (65). For the case where multi-point multi-harmonic forces are applied, such phase constraints introduce as many constraints as the number of unknown forcing amplitudes, R F , which may be expressed where C φ denotes the phase constraint matrix, similar to that in Eq. (32). Note that, in the previous case in Sect. 3, where quadrature harmonic forcing is considered, the phase constraint matrix restricts zero phaseshifts to the enforced harmonics; however, here, the phase constraint matrix enforces relationships between harmonic phase-shifts. Therefore, the number of unknowns (R H phaseshifts and R F forcing amplitudes) matches the number of equations (given by the energy balancing, i.e. Eq. (56), and the phase constraints, i.e. Eq. (66). Consequently, Eqs. (56) and (66) can be combined to give the energy balancing analysis for quadrature physical forcing cases, in the same form as Eq. (35) where This allows the phase-shifts of all harmonics,φ, and the physical forcing amplitudes,F, to be computed via known parameters obtained from NNM solutions.

Example 2: a quadrature physical forcing case
To formulate the HETA for a quadrature physical forcing case, the two-mode beam model, schematically shown in Fig. 2, is again considered. Here, a twoharmonic physical force is considered, given by which is applied to the part-span position at L 1 = 0.1L, such that the equations of motion are written 70b) where the nonlinear forces are given by Eq. (10). The parameters of the system are the same as those considered in Sect. 2.3 and Sect. 3.4, given in Table 1. Each modal displacement, q i , is approximated by three harmonics, namely the odd-numbered harmonics up to the 5th order, i.e.
such that As this example is of more practical relevance 5 , the details of the HETA construction are given here for reference. Firstly, the known energy-transfer coefficients are assembled into vector or matrix forms. Using Eqs. (20a), (26a) and (28a), the damping energytransfer vector is given by The diagonal matrix of forcing energy coefficients is defined by Eq. (53), with elements given by Eq. (48), i.e.
where the elements in this matrix may be computed using Eq. (20c). The phase-shift constraints, defined in Eq. (65), may be assembled into matrix form as Eq. (66), where the matrix of phase-shift constraint coefficients, C φ , and vector of phase-shifts,φ, are given where Combining equations from Eq. (73) to Eq. (78), the energy-transfer analysis may be formulated in the form as Eq. (67). Therefore, the unknowns, i.e. the physical forcing amplitudes,F , and the harmonic phaseshifts,φ, may be computed. Here the results obtained from HETA are compared with the forced responses obtained via the force appropriation technique, proposed in [25]. In implementing force appropriation, the same mono-point forcing, given by Eq. (69), is applied with the quadrature criterion satisfied. This is reflective of many practical modal tests where the number of forcing inputs is limited. In fact, it is often sufficient to use just a single-point excitation in identifying NNMs [5,[26][27][28][29]. Figure 8 shows the solved parameters, i.e. physical forcing amplitudes, |F 1 | and |F 2 |, and harmonic phase-shifts, i.e.φ 1,1 ,φ 1,3 ,φ 1,5 ,φ 2,1 ,φ 2,3 andφ 2,5 . In all panels of Fig. 8, the predicted results obtained via HETA are shown as solid blue lines whilst the simulated results obtained using force appropriation are shown as red lines with dots. The forcing amplitudes obtained from both techniques are indistinguishablesee panels (a) and (e) in Fig. 8. Excellent agreement is generally achieved between analytical predictions and numerical simulations for most phase-shifts. One discrepancy lies in the predictions ofφ 1,3 near the region where a strong internal resonance is shown between two linear modes (the loop region in panels (a) and (e)). This discrepancy may arise from the large phase-shift in this region, whilst this technique is derived based on the smallness assumption of the phase-shift. Nonetheless, the numerical results indicate a drastic change in the sign of phase-shift from positive to negative or vice versa, which is also captured by the results obtained from HETA. The other discrepancy lies in the prediction ofφ 2,1 near the natural frequency region; however, the trend over response frequency is again captured in good agreement.
In this case, some harmonics, namely u 1,3 , u 1,5 and u 2,5 , exhibit large harmonic phase-shifts near the internally resonant region. However, these harmonics have negligible amplitudes, as such the effect on the overall response is still minimal. In contrast, the fundamental components, i.e. u 1,1 and u 2,3 , which have large amplitudes and significant effect on the responses, exhibit small phase-shifts along the NNM branch-this can be used to estimate the applicability of energy balancing analysis, as with the case considered in Fig. 7. In turn, the accurately captured internally resonant region via force appropriation can also be interpreted by the small phase-shifts of large-amplitude fundamental components. To further improve the accuracy, increasing the number of forcing harmonics or excitation points can be considered [30] the harmonic phase-shifts,φ 1,1 ,φ 1,3 ,φ 1,5 ,φ 2,1 ,φ 2,3 andφ 2,5 , respectively. In all panels, the blue lines represent analytically predicted values (using the NNM data), the red lines (with dots) represent the numerical results, and the dotted-black line denotes the first linear natural frequency sis of nonlinear responses [6,10,21]. Using the energy balancing analysis to predict the existence of isolas has proven to be an efficient method [10]. However, as discussed in Fig. 3b, when responses contain multiple significant harmonics, predicting these isolas can be challenging-a key motivation for this work. In this section, an improved and more robust method is demonstrated in predicting the existence of isolas using the energy balancing analysis proposed in this paper. For comparison, the example system, considered in Sect. 2, is revisited. The nonlinear beam is modelled by a two-mode model, described by Eqs. (9) and (10), where the parameters are given in Table 1. A monopoint single-harmonic periodic forcing, expressed by Eq. (8), is applied to the part-span position at L 1 = 0.2L (the same excitation scenario as the case shown in Fig. 3b). Here, each mode is approximated by three harmonics, i.e. the odd-numbered harmonics up to the 5th order. As such, the formulation of META is equivalent to that given in Sect. 4.2 with F 2 = 0, and hence it is not re-formulated here.
The solved parameters, i.e. forcing amplitude, F 1 , and phase-shifts,φ 1,1 ,φ 1,3 ,φ 1,5 ,φ 2,1 ,φ 2,3 andφ 2,5 , are shown as blue lines in Fig. 9. In panel (a), the forcing amplitudes, solved via the single-harmonic META, are shown as a green line in comparison with that obtained via the multi-harmonic HETA. Near the natural frequency, the difference between results obtained via META and HETA are indistinguishable; as response frequency increases, the discrepancy between the two lines grows, indicating an increasing significance of the harmonics. Considering the forcing amplitude of F 1 = 0.18, HETA predicts one intersection, also shown in Fig. 3b, where the existence of an isola is not captured, whilst META predicts three intersections, labelled by solid dots and marked by (i), (ii) and (iii) in each panel of Fig. 9. This denotes that there are three resonant crossing points between the NNM branch and the forced responses. The perturbations, from NNM solutions to forced responses, for these three predictions can be evaluated by referring to the phase-shifts from panels (b) to (g). In each panel, it is shown that the phase-shift at point (i) is the largest of the three, whilst that at point (iii) is the smallest. Thus, it can be expected that the largest perturbation occurs at point (i), and the smallest one is at point (i).
To verify these results, the forced responses of the beam system are computed via numerical continuation and shown as red lines in the left panel of Fig. 10. Three HETA-predicted crossing points on the NNM branch are likewise denoted by solid dots, whilst the numerically obtained phase resonances, where the forceddamped motions are in quadrature with the forcing, are labelled by '×' signs. It can be seen that the resonances, predicted via HETA, show strong agreement with the numerically computed results, capturing the resonances on both the primary response curve and the isola. To aid comparison, the corresponding time-parameterised responses for the NNM solutions, HETA-predicted resonances-accommodating the NNM solutions with predicted phase-shifts via Eq. (18), and the numerically computed phase resonances are shown in the right-hand panels. In these panels, the perturbations from NNM solutions to resonances are accurately captured by the phase-shifts-HETA-predicted responses well match the numerically computed resonances. In addition, the levels of perturbations are well captured by the phase-shifts-the largest perturbation occurs at point (i) and the smallest one at (iii). Note that the current formulation of HETA does not allow for direct estimation of the stability of forced-damped responses. Nonetheless, the wellpredicted forced responses can be used as accurate initial solutions in stability studies. An alternative method, using NNMs to predict the stability of the perturbed forced responses, has been considered via a Melnikov analysis in [3].
In this example, predicting the existence of an isola is considered via HETA. Results show that the HETA can accurately predict all the resonant crossing points between the NNM branch and forced responses, showing an improved accuracy when compared with META. The perturbations from NNM solutions to resonances are shown to be accurately captured by the phase-shifts. This provides an approximation of the forced response by accommodating the NNM solutions with predicted phase-shifts.

Conclusions
In this paper, we have considered energy balancing across multiple harmonics of nonlinear normal modes (NNMs) for predicting resonances.
The energy balancing mechanism is established by accounting for the net energy transfer, arising from nonlinearity, external forcing, and damping, between harmonics of an NNM over a period of response. Combining the energy balancing principle with force reduction and phase-shift constraints, a semi-analytical framework has been developed for accurate and efficient predictions of resonances using NNM solutions. Two excitation scenarios, namely where the excitations are in quadrature with ( 1) harmonic displacements and ( 2) physical displacements, have been considered to cover practical applications.
We have then considered three specific application cases to demonstrate the application of this technique. For all these cases, the proposed computationally cheap technique showed excellent agreement with the computationally expensive numerical simulations. In the third example, the technique was used to predict the existence of an isola that contains multiple significant harmonics. The locations of resonances and the devia-tions from NNMs to resonances were both accurately predicted.
In implementing this technique, the applicability and accuracy of the proposed method can be directly evaluated via the computed harmonic phase-shift parameters. By accommodating the known NNMs with phaseshifts, an approximation of forced responses can be achieved, which could be used for further studies, e.g. stability and bifurcation analysis. Owing to the semi-analytical framework, this technique can be easily extended to account for more complex scenarios, e.g. multiple harmonics of multiple modes of NNMs, and multiple forcing harmonics at multiple locations, without adding much computational effort.
where P i, j and ψ i, j are the amplitude and phase of the forcing, respectively. Substituting the expressions for p i, j andu i, j into Eq. (86), gives where it has been assumed that the phase-shift,φ i, j , is small. The phase difference, ψ i, j −φ i, j , accounts for the effect of external forcing on harmonic u i, j -either an energy gain or an energy loss.
For the case where quadrature harmonic forcings are considered, i.e. Sect. 3, the phase of the harmonic forcing, ψ i, j , may be assumed to be equal the phase of the velocity,φ i, j . This further simplifies Eq. (88) to Assuming the displacement amplitude, U i, j , is known, but the forcing amplitude, P i, j , is unknown, the forcing energy term may be separated into the known, E † Pi,k , and unknown, P i,k , components, i.e.
where the dagger, • † , denotes a known term. Whilst, for the case where quadrature physical forcings are considered, as in Sect. 4, the phase relationship between the external forcing and the harmonic can be determined using the phase constraints Eq. (63) with In this context, the phase difference, ψ i, j −φ i, j , can be evaluated using the NNM solutions and thus, it can be seen as a known parameter.

A.3 Nonlinear energy-loss term
As with the forcing energy-loss term, the known and unknown components must be separated in the non-linear energy-loss term, Eq. (14c). Before considering the (unknown) phase-shifts,φ i, j , the responses are divided into the NNM responses (assumed to be known) and the perturbations (assumed to be unknown)i.e. q i =q i +q i where, as previously,q i represents the NNM response and whereq i is the perturbation to the ith modal displacement due to the forcing. Following this, the vector of modal displacements may be written q =q +q, and the harmonics are written u i, j =ū i, j +û i, j . Substituting these into the nonlinear energy-loss term, Eq. (14c), leads to Assuming the perturbation terms are small, the Taylor series expansion may be applied to the nonlinear force to give which may be expanded out to dN qi dq n q=qq n u i, j + N n=1 dN qi dq n q=qq n u i, j dt .
As the unperturbed response represents an NNM, the first term in the integral of Eq. (94) represents the net energy transfer due to the nonlinear coupling of an NNM response. As an NNM cannot exhibit a net energy transfer, it follows that Assuming that the perturbations,q n andû i, j , are small it also follows that the final term in the integral of Eq. (94) is negligible. Therefore, Eq. (94) may be further approximated to E Ni, j ≈ The perturbations,q n andû i, j , are now written in terms of the phase-shifts,φ i, j . To achieve this, Eq. (18) is used to write the jth harmonic of the ith mode as u i, j =ū i, j +û i, j = U i, j cos jωt +φ i, j +φ i, j , (97a) = U i, j cos jωt +φ i, j cos φ i, j −U i, j sin jωt +φ i, j sin φ i, j , (97b) where it has been assumed thatφ i, j is small. It therefore follows that As with the forcing energy-loss term, this may be separated into the known components (i.e. the NNM responses) and unknown components (the phase-shifts) as follows: Ni, jφ n,k , where the known component is written As in A.2, the dagger, • † , denotes a known term.

B Constructing energy balancing analysis for the example in Sect. 3.4
To construct the energy balancing equations using Eq. (31) for the beam example in Sect. 3.4, the coefficients and unknowns are first assembled into vector/matrix form. The damping energy vector is defined, using Eqs. (20a), (26a) and (28a), as Next, using (26b), (26c), (28b) and (28c), the matrix of forcing energy coefficients, and vector of unknown harmonic forcing amplitudes, are written as As the third harmonic of the first mode and the first harmonic of the second mode are unforced, the force reduction matrix and reduced vector of forcing amplitudes may be written where the forcing amplitudes P 1,1 = F 1 and P 2,3 = F 2 , from the equations of motion, i.e. Eq. (36), have been substituted. Using the force reduction matrix, the reduced forcing energy coefficient matrix is given bỹ where Eq. (20b) has been used to define the forcing energy coefficients.