An extended energy balance method for resonance prediction in forced response of systems with non-conservative nonlinearities using damped nonlinear normal mode

The dynamic analysis of systems with nonlinearities has become an important topic in many engineering fields. Apart from the forced response analyses, nonlinear modal analysis has been successfully extended to such non-conservative systems thanks to the definition of damped nonlinear normal modes. The energy balance method is a tool that permits to directly predict resonances for a conservative system with nonlinearities from its nonlinear modes. In this work, the energy balance method is extended to systems with non-conservative nonlinearities using the concept of the damped nonlinear normal mode and its application in a full-scale engineering structure. This extended method consists of a balance between the energy loss from the internal damping, the energy transferred from the external excitation and the energy exchanged with the non-conservative nonlinear force. The method assumes that the solution of the forced response at resonance bears resemblance to that of the damped nonlinear normal mode. A simplistic model and full-scale structure with dissipative nonlinearities and a simplistic model showing self-excited vibration are tested using the method. In each test case, resonances are predicted efficiently and the computed force–amplitude curves show a great agreement with the forced responses. In addition, the self-excited solutions and isolas in forced responses can be effectively detected and identified. The accuracy and limitations of the method have been critically discussed in this work.


Introduction
The dynamic responses of structures with various nonlinear characteristics have become an important research field that is able to improve the design of future engineering structures. Most of the nonlinearities involved in a complicated engineering system can be generally classified by two main types: conservative nonlinearities and non-conservative nonlinearities. Geometric nonlinearity can be seen as an example of a conservative nonlinearity if viscoelasticity of material is not considered and it has been widely explored by many researchers [2]. Frictional contact is instead an example of non-conservative nonlinearities, and it is important for many engineering applications. In some systems, the frictional contact, considered as a dissipative nonlinearity, is able to provide the external damping to reduce the vibrational stress under resonance. The frictional damping has been widely used and studied in the turbomachinery industry [15,19,20]. On the other hand, non-conservative frictional forces can also provide energy to a system, as in the case of brake-squeal models where self-excited vibrations can be induced by full sliding frictional belts [10]. The energy dissipated by or gained from the nonlinear force makes the system non-conservative and classical nonlinear modal analysis not applicable.
In the context of nonlinear modal analysis, the nonlinear modes are computed and used to study a dynamic system. The nonlinear modes are considered as dynamic responses of a conservative autonomous system with nonlinear characteristics. In a linear system, the responses of free vibration are defined as normal modes. Unlike the linear normal modes, superposition and orthogonality properties are not valid for nonlinear modes. The first definition of nonlinear modes given by Rosenberg [23] defined the nonlinear normal mode (NNM) as a vibration unison, where all points in a structure reach their equilibrium position and their extreme position simultaneously. However, the definition proposed by Rosenberg cannot explain internal resonances. After that, many further developments of the NNM have been proposed in [21,27]. Then, different methods are used to derive the NNM (e.g. normal form [12,28] or multiple scales [17]). Analytical analysis of NNM unmasks the dynamics of nonlinear autonomous system including energy dependence of NNM, interaction and bifurcation of NNM. To compute the NNM of large scale structures, various numerical methods were attempted [22], i.e. a shooting method with continuation procedure [18]. However, most of works mentioned above considered undamped NNM with a periodic solution, i.e. system with conservative nonlinearities. When non-conservative nonlinear characteristics are taken into account, the NNMs are no longer periodic and defined as damped nonlinear normal mode (dNNM) [24]. Shaw and Pierre defined the dNNM using the concept of invariant manifold [25]. Then, Laxalde and Thouverez described a numerical method to compute the dNNM using the extension of complex nonlinear modes by solving complex eigenproblems [16]. Thanks to the methods proposed in [16], the dNNM of full-scale structures with frictional contact can be efficiently computed. To obtain a periodic based dNNM, Krack [14] developed a numerical method, namely extended periodic motion concept (EPMC), where an artificial modal damping is introduced into the system to compensate the energy dissipated by the nonlinear forces. Compared to the numerical method described in [16], this EPMC method is able to provide a periodic solution and energy dissipated can be effectively quantified by the artificial modal damping ratio. Therefore, EPMC is considered as an efficient method to investigate frictional systems and used in this work [11].
In an engineering perspective, the resonances in forced responses are always considered as an important topic for the issue of high cycle fatigue. At resonance, the dynamic response has a similar solution than that of the nonlinear mode at a specific energy level. In fact, backbone curves are sometimes seen as the collection of resonance points for a varying value of the excitation force. Based on these considerations, an analytical method named energy balance method (EBM) has been attempted by Hill [8] to predict where the resonant solution of the forced response crosses the underlying nonlinear mode. Before this EBM, similar energy transfer analysis has been described in [5,13,29]. Grenat used similar energy analysis to calculate the non-conservative NNM based on NNM [6]. The EBM is originally developed for systems with conservative nonlinearities using the theory of NNM and successfully applied to a two-degree-of-freedom (DoF) nonlinear oscillator [8]. EBM is also able to capture isolas in forced responses and isolated branches of nonlinear modes [9]. The method relies on the fact that, for any steady-state response, the net energy transferred into a system over one period of vibration must be zero [8]. Moreover, the work of Cenedese and Haller [4] uses Melnikov function analysis that defines the mathematical conditions under which the EBM is valid.
A novel extended energy balance method (E-EBM) is described in this work and used to predict the resonance in systems with non-conservative nonlinearities. In such systems with non-conservative nonlinearities, the energy dissipated or gained by nonlinear forces needs to be taken into consideration. To do so, the dNNMs are used instead of the NNMs. E-EBM is tested and validated for three different test cases: a one-DoF system with frictional contact, a two-DoF system interacting with a moving belt and a full-scale joint beam with contact interface.
The paper is organised as follows: the general nonlinear dynamic equations for a forced system and autonomous system are given in Sect. 2; the method to compute dNNMs is explained in Sect. 3; the E-EBM is introduced and explained in Sect. 4; the results from three test cases are investigated in Sects. 5-6; finally, the benefits and the limitations of this E-EBM are critically discussed in Sect. 7 followed by the conclusions.

System with non-conservative nonlinearities
The objective of this section is to describe the equation of motion for a general nonlinear system. The general dynamic equation to compute the nonlinear forced response is expressed as Eq. (1), where M is the mass matrix; K is the stiffness matrix; C is the viscous damping matrix. Q f (t) is the displacement of the system, the underline of the symbol is used to represent the vector, and variable t does not appear explicitly in all equations in this work. F e is the external excitation force. F nl is the non-conservative nonlinear force. is the excitation frequency; γ is the excitation forcing level; ϕ is the absolute phase of the excitation force. The forced responses of a nonlinear system are computed by solving Eq. (1).
The nonlinear modes are defined as the solutions of the autonomous nonlinear systems. In an autonomous system, the external excitation force and viscous damping term are not considered. Therefore, the general equation of motion to calculate the nonlinear modes is represented in Eq. (2), where Q is the displacement along the nonlinear mode of the system.

Nonlinear modal analysis
The general dynamic equations for forced responses and nonlinear modes have been explained in the previous section. The objective of this section is to describe the methodology used to compute the dNNMs for systems with non-conservative nonlinear characteristics. For such systems, the dynamic solutions are typically not periodic and the periodic-based definition of nonlinear modes is not applicable in the classic sense. Therefore, the EPMC approach [14] is used to compute the dNNMs. Based on the concept of dNNM, the nonlinear modes and other modal properties are modal energy-dependent and vary with the level of energy within the system. The energy dependency is achieved by introducing a new parameter, modal amplitude α . The modal amplitude α is used to quantify the modal energy within the system. The modal displacements of the system are expressed as Q = α · Q 0 , where Q 0 is the mass normalised modal displacement. Furthermore, since the total energy of the autonomous system is not conservative, an artificial damping term is introduced into the system to balance the energy dissipated or gained from the nonlinear forces F nl . The artificial damping term is presented in Eq. (3), where ω 0 is the energy-dependent resonant frequency of studied dNNM and ζ is the modal damping ratio.
The modal damping ratio ζ is either positive or negative and depends on the nonlinear force. In a system with dissipative nonlinear forces, the modal damping ratio is positive and the total energy of the autonomous system is dissipated by the nonlinear force. However, in a nonlinear system showing self-excited vibration, the modal damping ratio is negative. In this case, the energy is injected into the system due to the nonlinear force. After introducing the modal amplitude and modal damping, the dynamic equation of the nonlinear autonomous system is rewritten as in Eq. (4). The mass normalised modal displacements Q 0 , the resonant frequency ω 0 and the modal damping ratio ζ vary with the modal amplitude α. In an autonomous system, the absolute phase of the mass normalised solution Q 0 is arbitrary. Therefore, to solve the equation above, a phase normalisation constraint and a mass normalisation constraint are required. Both constraints are applied in a similar way with [26].
The harmonic balance method (HBM) with alternating frequency/time method [3] is used to solve Eqs. (1) and (4). In this method, the continuous properties in time domain Q are discretised based on Fourier series into frequency domain (Q c p + Q s p i) as shown in Eq. (5), where N h is total number of the harmonics used and p is order of harmonic. The dynamic equations are solved in frequency domain by using HBM. The ω can be excitation frequency in a forced responses and resonant frequency ω 0 in the computation of dNNMs. The detailed description of the method can be found in [26].
In forced responses analysis, the evolution of the nonlinear forced responses with respect to the excitation frequency is expected, whereas in nonlinear modal analysis, the dNNMs and other modal properties are required for a range of modal amplitude. The continuation method is the numerical technique to track the evolution of the system for a specific chasing parameter. Readers are invited to refer [1] for a detailed description of the formulations. In this work, a secant predictor and arc-length corrector are used in the continuation method.

Extended energy balance method
The methodology to compute the dNNMs for a system with non-conservative nonlinearities is explained in Sect. 3. The objective of this section is to describe the E-EBM to predict the resonance in forced responses. The idea behind EBM is that the solution of the nonlinear forced response and that of the nonlinear mode share similarities and, for small enough levels of damping and forcing, these solutions can be assumed to be identical and have been named 'resonant shared solutions' [8]. This concept also applies to dNNM in that the solution of the autonomous system bears resemblance with that of the nonlinear forced response at resonance. In this work, we extend the application of EBM from NNMs to dNNMs.
In systems with non-conservative nonlinear characteristics, the dNNMs of the system are computed firstly and the E-EBM is applied to predict the resonances in forced responses. The dynamic equation to calculate the nonlinear forced responses is shown in Eq. (1). The equation of motion for the computation of dNNMs is represented by Eq. (4). The solutions at resonance Q f in nonlinear forced response are assumed to be same as the dNNMs with certain level of modal amplitude α, when the excitation frequency equals the resonant frequency ω 0 . The resonant shared solutions are here called ψ. An example of resonant shared solution is shown in Fig. 1, where the nonlinear forced response and the dNNM are shown. In other words, the dNNM with particular modal amplitude α is regarded as the solution of nonlinear forced response at resonance with a particular excitation forcing level, which is expected to be predicted based on the E-EBM.
Therefore, under resonance, the dynamic equation of the forced system and autonomous system are rep-  7) and (8).
For each level of the modal amplitude α, the resonant shared solution ψ, the resonant frequency ω 0 and the modal damping ratio ζ are computed using the numerical methods described in Sect. 3. For a given system with known viscous damping term C and the vector of the excitation force, the excitation forcing level γ and the phase of the excitation force ϕ are expected to be determined for each level of modal amplitude α using such E-EBM. Recalling Eqs. (7,8) and assuming that the solution of the nonlinear forced response equals the resonant shared solutions ψ, the subtraction of Eqs. (8) to (7) must hold in its weak form.
The energy transfers caused by damping terms E d and forcing term E f (γ , ϕ) are calculated by integrating the damping force and excitation force times the velocity over one vibrational period 2π/ω 0 , as shown in Eqs. (9) and (10).
For each solution point along the dNNM curve (see an example in Fig. 1), i.e. for each value of the modal amplitude α, the excitation force that generates a resonance at that point can be found by the E-EBM. In Fig. 2, a typical scenario is provided; the energy transfer caused by damping terms E d is constant for a certain level of modal amplitude α, whereas the energy Fig. 2 The energy transferred by damping (E d ) for a fixed modal amplitude α and the energy transferred by forcing terms (E f (γ , ϕ)) for different excitation levels γ with respect to the forcing phase ϕ exchanged by the external excitation E f (γ , ϕ) varies with the forcing phase ϕ and the excitation forcing level γ . To determine the values of ϕ and γ is to find the single intersection between the curve of E f (γ , ϕ) and the constant value of E d that will occur at the maximum (if If this operation is applied to each value of the modal amplitude α along the dNNM, the classic force-amplitude curve can be built.
When the nonlinear force is non-conservative, the modal damping ratio becomes nonzero. Depending on the value of modal damping ratio, four different situations are possible: (i) ζ > 0: When the modal damping ratio is positive, the non-conservative nonlinear force in the system provides external damping and the total energy of the autonomous system is dissipated. In this case, the total energy injected by the external excitation force is required to balance the summation of the energy dissipated by viscous damping and external damping generated by the nonlinear force. (ii) ζ < 0 and E d > 0: The modal damping ratio is negative meaning that the nonlinear force is injecting energy into the system. Although, the energy dissipated by damping terms is higher than the energy injected by nonlinear force leading to a positive E d . In this case, the external excitation force has to provide energy to the system to balance the difference between the viscous damping and the negative modal damping generated by the nonlinear force. (iii) ζ < 0 and E d = 0: The modal damping ratio is negative, and it perfectly balances the energy dissipated by the viscous damping term, so E d equals zero. In this case, self-excited solutions can be found and no external force is required to observe a periodic solution, which coincides with the limit cycle oscillation. (iv) ζ < 0 and E d < 0: The modal damping ratio is negative, and it overcomes the energy dissipated by the viscous damping term, so E d is negative. In this case, for periodic solution to exist, the external force has to take energy out of the system.
In the following sections, three different test cases are investigated: the first two cases refer to the (i) situation and the third one to the (ii-iv) situations.

Systems with dissipative nonlinearities
In this section, a one-DoF system and a full-scale joint beam structure with frictional contact are studied. The frictional contact involved is able to provide the frictional damping, and the total energy of the system is dissipated. In the system with dissipative nonlinearities, a positive modal damping ratio is expected. The dNNMs are calculated, and E-EBM is applied. The results are numerically validated using the nonlinear forced responses.

One-DoF system with frictional contact
A one-DoF system with frictional contact is tested as shown in Fig. 3. In this test case, a simple mass-spring model with tangential friction is considered. A Jenkins element with constant normal load is used to model the frictional force [31]. In this simplistic contact model, there are two contact conditions, which are sticking and sliding. The parameters of the system are given in Table 1. A simple harmonic excitation force is considered in horizontal direction.
The resonant frequency ω 0 , the modal damping ratio ζ and the predicted excitation forcing level γ are plotted with respect to the amplitude of the tangential displacement as shown in Fig. 4. The resonant frequency for sticking contact condition is 213.3 rad/s. When the amplitude is smaller than 5 mm, the mass is sticking to the ground and the system is linear. The resonant frequency ω 0 remains as a constant; the modal damping ratio ζ is zero. In force-amplitude plot (Fig. 4c), the predicted excitation forcing level γ is linearly proportional to the amplitude. An increase in amplitude leads to a sliding contact condition. A softening effect can be observed and the resonant frequency ω 0 starts to decrease when the sliding contact condition occurs. The total energy of this autonomous system is dissipated caused by the rubbing motion between the mass and ground. A positive modal damping ratio ζ is achieved, and maximum frictional damping is around 4.1%. The predicted excitation forcing level γ increases nonlinearly with the amplitude. To validate the predicted excitation forcing level γ (force-amplitude plot), nonlinear forced responses of this system under different excitation forcing levels (1 N, 20 N, 40 N, 50 N and 60 N) are computed. The numerical validation is completed and shown in Fig. 5. On the left, the evolution of amplitude with respect to the predicted excitation forcing level (forceamplitude plot computed from the E-EBM) is represented by black curve. For given excitation forcing level, the solutions of the system are found and marked in the figure by different colours. On the right, the evolution of amplitude with respect to the resonant frequency (resonant frequency-amplitude plot) as well as the nonlinear forced responses is shown. The 3D view (frequency-force-amplitude) of the dNNM and several nonlinear forced responses is represented in Fig. 5b. Generally, the E-EBM is able to provide accurate prediction of the amplitude of the system under resonance for a given excitation force as shown in Fig. 5a. To achieve a better evaluation of the accuracy, the hysteresis loops under resonance from dNNM and nonlinear forced responses are shown in Fig. 6. From the figure, the hysteresis loops from dNNM are almost identical to that from nonlinear forced responses under resonance. To further demonstrate the accuracy of the E-EBM, for each excitation forcing level, the amplitude |Q| calculated based on the E-EBM is compared with the amplitude obtained from the resonance |Q f | in nonlinear forced responses. The resonant amplitude and percent-  Table 2. From the table, the resonant amplitude calculated from the E-EBM shows a great agreement with the resonant amplitude obtained from nonlinear forced responses. The maximum percentage difference achieved is around 1.36%, when the excitation force is 60 N.
The error obtained in this case is mainly from the numerical method. There is only mode in this one-DoF system. The classic HBM is used to compute the periodic solutions, which are decomposed into Fourier series truncated to certain order of harmonics. In nonlinear forced response, the first-order harmonic excitation is applied to this one-DoF system, whereas, in non- Fig. 7 Percentage differences for different harmonics and different excitation forcing level for one-DoF system with friction linear modal analysis, the harmonics up to fifth order are used to compute the dNNM. The physical theory behind the E-EBM is that the energy exchanged by the damping terms is balanced by the energy exchanged by the excitation force. When the damping terms are completely equal to the excitation force, the E-EBM is expected to theoretically predict the actual resonance with zero error for this one-DoF system. However, this cannot be achieved in nonlinear systems. In a nonlinear system, harmonics with higher order can be excited, even the excitation force is in first order. In Fig. 7, the percentage differences for different harmonics (used in computation of dNNM) and different excitation forcing level are shown. The best prediction or smallest percentage difference can be observed when only the first order is used in the computation of dNNM. Apart from that, according to Table 2 and Fig. 7, the percentage difference increases with the excitation forcing level. In this one-DoF system, the higher excitation force leads to more sliding and the nonlinear force has higher impact on the periodic solutions. When the excitation force is 1 N, as shown in Fig. 6a, the system is almost linear (sticking contact condition) and harmonics with higher order are negligible. Therefore, percentage error is almost zero, which can be considered as an idea situation. In summary, the differences in amplitude between the predicted resonance and actual resonance are mainly caused by the higher order of harmonics through the numerical computation.

Full-scale joint beam with contact interface
A full-scale joint beam with contact interface is tested by the E-EBM, and the structure is shown in Fig. 8. The contact interface between two beams is coloured in yellow. This joint beam contains two 3D beams, which are modelled by linear brick element. This full-scale joint beam model has been well explained in [33], and parameters of the model are given in Table 3. There are 3003 nodes in total within the model and including 121 contact pairs (242 contact nodes) within the contact interface. For nonlinear dynamic analysis, nonlinear static analysis is completed in the first place by using finite element software. The nonlinear static solution and contact pressure within the contact interface are obtained and used for further dynamic analysis. A 3D node-to-node contact element which is made of two 1D Jenkins elements in tangential direction and a coupling spring in normal direction is used to simulated the frictional force within the contact interface, as shown in Fig. 9. In this contact model, k t1 and k t2 are contact stiffness in tangential direction and k n is normal contact stiffness. The gap or normal force N 0 for each pair of contact nodes is computed through nonlinear static analysis. The detailed description for this 3D node-to-node contact element can be found in [32]. There are three contact conditions, which are sticking, sliding and separation. In the first place, the frictional contact within the contact interface is modelled by linearised contact stiffness. Then, the first four modes of the linearised joint beam are computed and shown in Fig. 10. Figure 10a, c shows the linearised mode shapes of the first and second in-plane bending modes. The first  Poisson's Ratio 0.32

Fig. 9
A 3D node-to-node contact model: two 1D Jenkins elements in tangential direction and a coupling spring in normal direction edge-wise flapping mode is shown in Fig. 10b, and the first torsion mode is represented in Fig. 10d. The first edge-wise flapping mode is studied in this test case. The resonant frequency ω 0 , the modal damping ratio ζ and the predicted excitation forcing level γ for the first edge-wise flapping mode are represented with respect to the displacement of the tip of the beam as shown in Fig. 11. At the low level of vibrational amplitude, the contact nodes within contact interface stay in sticking contact condition and system is linearly. When the system vibrates at lager amplitude, some contact nodes start to slide and the total energy of the system is dissipated. The results shown in Fig. 11 are typical to the structures with frictional contact as explained in Sect. 5.1.
To the validate the predicted excitation forcing level in Fig. 11c (force-amplitude plot), the nonlinear forced responses with excitation force 0.1 N, 0.2 N, 0.4 N, 0.8 N and 2 N are computed and shown in Fig. 12a. The 3D view of the dNNM and nonlinear forced responses is shown in Fig. 12b. From the figures, the amplitude calculated from the E-EBM for given excitation forcing levels is very close to the resonant amplitude obtained from the nonlinear forced responses. The comparison of resonant amplitude calculated from the E-EBM and nonlinear forced responses is listed in Table 4. The maximum percentage difference is around 0.76%. To further demonstrate the accuracy of the method, the maps of contact status and energy dissipated within contact interface are compared between the dNNM and nonlinear forced responses under resonance. There are 121 contact nodes on each beam. In Fig. 13a, b, the contact status is plotted within the contact interface, where black squares represent the sticking condition and red diamonds represent sliding condition. The map of contact status from the dNNM looks similar with the one from nonlinear dynamic response. As for the map of energy dissipated, it is computed based on each contact node by calculating the area within hysteresis loop and shown in Fig. 13c, d. The total energy dissipated from dNNM and nonlinear forced responses is 0.228 J and 0.241 J. From the distribution, the nodes near the corners of the map shows the maximum damping ability (more energy dissipated). By comparing the contact status and energy dissipation, the results show a great agreement between the predicted resonance using dNNM and resonant solution in nonlinear forced response.
As discussed in Sect. 5.1, the main error shown in the E-EBM comes from the numerical method used in the computation the periodic solutions either in dNNM or nonlinear forced responses. In these test cases, harmonics of the same order are used in HBM. In Sect. 5.1, a one-DoF system is tested and there is only one mode within that system. However, in this test case, a fullscale joint beam is analysed. For this 3D structure, the forced responses are normally coupled by more than one modes. The amplitude raised by other modes cannot be taken into consideration. Therefore, the differences in resonant amplitude are mainly caused by the other modes within the system.

System with self-excited solutions
The third test case is a two-DoF system interacting with a moving belt, and the description of such system is shown in Fig. 14. A similar model has been introduced to study the self-excited vibration in brakes [7,10]. The variation presented here has been exploited to simulate the tip rubbing problem between the blades and cas-  ing within aero-engines [30]. In [30], a contact model with only two contact conditions, which are separation and sliding, is used. The parameters of the system are given in Table 5. A simple harmonic excitation force is considered in normal direction.
In [30], the undamped NNM of the system has been studied and self-excited solutions are found. The modal interaction between the first and second modes occurs through a 1:2 internal resonance, when the resonant frequency of the first mode coincides with half of the resonant frequency of the second mode. Underneath this internal resonance, two self-excited solutions are found.  In this section, dNNM is used to find the branch of all possible self-excited solutions; this branch lies underneath the internal resonance and connects two sides of the undamped NNM. In addition, the isolated branches of forced periodic solutions as well as the res-onant shared solutions are detected and predicted using the E-EBM.
The amplitude of tangential displacement with respect to the resonant frequency ω 0 from nonlinear modal analysis of the system is shown in Fig. 15. As shown in the figure, the undamped NNM (black dash curve) shows the modal interaction between the first mode and second-order super-harmonic of the second mode (a peak when resonant frequency ω 0 is around 1.13). The solutions represented by undamped NNM are conservative, and the modal damping ratio is zero. Therefore, the self-excited solution cannot be observed in undamped NNM. Applying the numerical method in Sect. 3, the dNNM of the system is computed and it is shown by the black solid curve in Fig. 15. A branch of dNNM (black solid curve) is observed under the modal interaction, and the solutions represented by this branch are non-conservative and a negative modal damping ratio is found. When the viscous damping applied to the system can completely balance the negative modal damping (E d = 0), the self-excited solutions can be observed. In other words, the branch of the dNNM can be considered as the branch of the self-excited solutions as shown in Fig. 16b, i.e. when the viscous damping terms c 1 and c 2 are 0.06 Ns/m, there are two self-excited solutions of the system (grey diamonds in Fig. 16b). When a higher viscous damping is applied into the system, the two self-excited solutions will move inside and eventually merge to one solution when viscous damping terms c 1 and c 2 reach 0.07 Ns/m. If the viscous damping added into the system is greater than the critical viscous damping (0.07 Ns/m), the self-excited solutions cannot be observable in the system.
The mode shapes of the dNNM with different amplitude are shown in Fig. 16. In Fig. 16a, there are five points highlighted by grey markers: A and E are the bifurcation points between the undamped NNM and dNNM; B and D are two self-excited solutions when viscous damping terms c 1 and c 2 reach 0.06 Ns/m. The mode shapes or periodic solutions are demonstrated by Q n − Q t plot, where Q n and Q t are displacement in normal and tangential directions, respectively. The first mode has a strong normal component, whereas the second mode has strong tangential component. In all periodic solutions, the 1:2 internal resonance is observable in that the tangential component has double the period than the normal one. The modal damping ratios at point A and E are zero and the area of their trajectories is null because they are the points where the dNNM bifurcates The resonant frequency ω 0 , modal damping ratio ζ and the predicted excitation forcing level γ of the dNNM are plotted with respect to amplitude of the tangential displacement as shown in Fig. 17, when the viscous damping terms c 1 and c 2 are 0.06 Ns/m. Two self-excited solutions of the system are marked by grey diamonds in Fig. 17. A negative modal damping ratio ζ is observed as shown in Fig. 17b. In this case, the contact force makes extra work on the system and total energy of the autonomous system is increased due to the non-conservative contact force. The force-amplitude plot is shown in Fig. 17c. When the excitation forcing level is zero, there are two resonant shared solutions between the nonlinear forced responses and dNNM, which are the self-excited solutions described above. For a larger excitation forcing level (0 < γ < 0.0025), it is expected to find four resonant shared solutions. In this case, two isolated branches of the nonlinear forced responses are expected. When the excitation forcing level γ is great than 0.0025 N, there are only two resonant shared solutions observed leading to one isolated  Fig. 18. For each subfigure in Fig. 18, on the left side, the force-amplitude plot computed with E-EBM is represented by a black curve; on the right side, resonant frequency-amplitude plot and the nonlinear forced responses are shown. The self-excited solutions are marked by grey diamonds. The predicted resonant shared solutions are highlighted by black markers. The 3D view (frequency-force-amplitude) of the dNNM and several nonlinear forced responses is represented in Fig. 19. In Fig. 18a, b, two isolated branches of the nonlinear forced responses and four resonant shared solutions are observed when excitation forcing level is 0.001 N and 0.002 N, respectively, which corresponds to the prediction in Fig. 17b. Recalling the four possible situations introduced in Sect. 4, two different behaviours can be observed in this these cases (0.001 N and 0.002 N). The two inner most resonant shared solutions have a negative E d , which means the excitation force has to take the energy out of the system (see situation (iv) in Sect. 4). Instead, the two outer most resonant shared solutions fall into the category of positive E d (see situation (ii) in Sect. 4). The accuracy of the prediction for the case with excitation force 0.002 N is slightly lower than the case with 0.001 N. Figure 18c, d shows one isolated branch for each nonlinear forced response, when excitation forcing level is 0.0025 N and 0.003 N, respectively. For the case with 0.0025 N, there are three predicted resonant solutions: two outer resonant shared solutions and one in the middle. The two outer most belong to situation (ii) in Sect. 4, whereas the one in the middle belongs to situation (iv). As for Fig. 18d, there two resonant shared solutions predicted (situation (ii) in Sect. 4). All self-excited solutions (marked by grey diamonds in Fig. 18) are examples of the situation (iii) in Sect. 4. The accuracy of the E-EBM is lower than that of the first two test cases, especially in the cases with 0.0025 N and 0.003 N. To achieve a better understanding, the case with 0.0025 N is selected and periodic solutions are shown in Fig. 20.
The dNNM (black curve) and nonlinear forced response (orange curve) for the case with 0.0025 N are shown in Fig. 20. All three predicted resonant shared solutions are highlighted by the black squares and two self-excited solutions in grey diamonds. The periodic solutions of the two outer predicted resonant solutions are shown by the black dash dot curve in the figure. The orange dash dot curves represent the solutions obtained The predicted resonance shared solution in the middle is underneath an internal resonance. In this case, the first mode and second mode are strongly coupled, leading to a lower accuracy. In other words, the accuracy of the E-EBM cannot be guaranteed under internal resonance.
As described in Sect. 5.1, the resonance can be predicted based on the fact that net energy exchanged by the damping terms and excitation force is zero. When the damping terms mathematically equal the excitation force, the method can provide perfect predictions. In this test case, the excitation force is a simple harmonic excitation applied in the normal direction, whereas, from the periodic solution, the displacement in both normal and tangential directions is found due to the internal resonance. This similarity between excitation force and damping term cannot be achieved in case of internal resonance, thus explaining the inaccurate predictions observed in this test case. When the excitation force has a different shape than the resonant shared solution, some part of the energy can be injected into other modes instead of the mode studied. This discrepancy between the shapes of the excitation force and resonant shared solutions is considered as the major source of the inaccuracy. In addition, the classic HBM is used as the numerical method to solve the nonlinear dynamic equation. In this test case, the frictional contact shows a great nonlinearity and non-smoothness. The harmon- Fig. 19 A 3D view of damped nonlinear normal mode (black curve) and several nonlinear forced responses for two-DoF system interacting with a moving belt ics with higher order also have a great impact on the accuracy of the prediction.

Discussion
The EBM is extended to systems with non-conservative nonlinearities by using their dNNMs instead of their NNMs. The E-EBM is used to predict the resonant condition (resonant amplitude, resonant frequency and modal damping under resonance) for a given excitation force. This method has been applied to three different test cases, and results are numerically validated by the nonlinear forced responses. From three test cases studied, the advantages and limitations of the method will be critically discussed in this section.
Firstly, the E-EBM is able to predict the resonant condition for systems with non-conservative nonlinearities, i.e. system with frictional contact. According to the results from the first two test cases in Sect. 5, the periodic solution of the nonlinear forced responses under resonance can be effectively predicted with higher accuracy. The maximum percentage difference in the resonant amplitude for first two cases is 1.36% and 0.76%, respectively. A negligible difference is obtained between the predicted resonant solutions and solution from nonlinear forced response. Secondly, the classic force-amplitude plot can be generated during the computation of the dNNMs. It is able to identify the resonant amplitude according to the excitation forcing level in an efficient way. Then, for systems with self-excited mechanism, the E-EBM is able to easily find the self-excited solutions. The method makes the detection and track of the isolated forced responses applicable. Finally, the EBM was only used in simple systems (i.e. system with few DoFs) in the literature.  This work has proven that the E-EBM can be applied for real engineering problems by using their dNNM, i.e. joint structures with contact interface.
The limitations of the E-EBM have been identified in this paper from the investigation of three test cases. The inaccuracy of the method can be mainly considered from two different aspects. The effect caused by the numerical method is considered in the first place. All three test cases are solved by using the classic HBM method. In HBM, the continuous solutions are discretised by Fourier series truncated to certain order. When the solution of an system has a strong nonlinear characteristic, higher order of harmonics shows a great impact on the accuracy of the method (see Fig. 7). Secondly, the effect from other modes can lead to a lower accuracy in the prediction of the resonance. In a complicated structures, other mode can also contribute to the resonant amplitude in a forced response. When the resonant solution from the forced responses shows a strong coupling between two modes, a lower accuracy can be found for the prediction of the resonant amplitude.
Actually, there is no proof that the method can work when an internal resonance is present. One can imagine that when there is an internal resonance, a small forcing on one mode will produce a large amplitude response on the coupled mode. In that case, the method can fail at predicting this response due to the orthogonality between such modes. However, in the case tested in this work, the E-EBM was applied to the branch of self-excited solutions underlying the internal resonance and not at the peak of the internal resonance where only the coupled mode would have been responded. In the computation of dNNM, the solutions calculated is already the mixture of the two modes involved in the internal resonance; thus, the method is able to provide quite accurate predictions from the dNNM.
Frictional contact is a common type of nonconservative nonlinear force and has been widely explored in the literature. All three test cases in this work can be seen as system with frictional contact. However, other types of non-conservative nonlinearities are expected to be tested using dNNM coupled with E-EBM.

Conclusion
The major contribution of this paper is the extension of the energy balance method to systems with non-conservative nonlinearities. A one-DoF system with frictional contact, a two-DoF system interacting with a moving belt and a full-scale joint beam model are tested. The method has been validated in a numerical way by comparing the predicted resonant solutions using the extended energy balance method to the actual resonant solutions obtained from forced responses. Generally, the method is able to predict the resonance for a given excitation force without the computation of the nonlinear forced response. The classic forceamplitude plot can be constructed during the computation of the damped nonlinear normal mode. Besides the simplistic model, this work also attempts to use the extended energy balance method to solve a real engineering problem. In addition, the main features of the isolated branches for a forced response can be directly detected and tracked using the method. However, the accuracy of the method is limited by the contribution of the other modes in the forced responses, especially in the case of internal resonance. The accuracy of the method is also affected by the numerical method used to solve the nonlinear dynamic problem. The extended energy balance method coupled with damped nonlinear normal mode can be used for engineering design problem to dramatically reduce the computational cost while ensuring an accurate prediction of the resonance in the future.
Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.