Detecting internal resonances during model reduction

Model order reduction of geometrically nonlinear dynamic structures is often achieved via a static condensation procedure, whereby high-frequency modes are assumed to be quasi-statically coupled to a small set of lower frequency modes, which form the reduction basis. This approach is mathematically justifiable for structures characterized by slow/fast dynamics, such as thin plates and slender beams, and has been shown to provide highly accurate results. Nevertheless, selecting the reduction basis without a priori knowledge of the full-order dynamics is a challenging task; retaining redundant modes will lead to computationally suboptimal reduced-order models (ROMs), while omitting dynamically significant modes will lead to inaccurate results, and important features such as internal resonances may not be captured. In this study, we demonstrate how the error associated with static condensation can be efficiently approximated during model reduction. This approximate error can then be used as the basis of a method for predicting when dynamic modal interactions will occur, which will guide the reduction basis selection process. Equivalently, this may serve as a tool for verifying the accuracy of ROMs without the need for full-order simulations. The proposed method is demonstrated using a simple oscillator and a finite element model of a clamped–clamped beam.

EN, 0000-0003-2066-0869; TLH, 0000-0002-4125-7895 Model order reduction of geometrically nonlinear dynamic structures is often achieved via a static condensation procedure, whereby high-frequency modes are assumed to be quasi-statically coupled to a small set of lower frequency modes, which form the reduction basis. This approach is mathematically justifiable for structures characterized by slow/fast dynamics, such as thin plates and slender beams, and has been shown to provide highly accurate results. Nevertheless, selecting the reduction basis without a priori knowledge of the full-order dynamics is a challenging task; retaining redundant modes will lead to computationally suboptimal reducedorder models (ROMs), while omitting dynamically significant modes will lead to inaccurate results, and important features such as internal resonances may not be captured. In this study, we demonstrate how the error associated with static condensation can be efficiently approximated during model reduction. This approximate error can then be used as the basis of a method for predicting when dynamic modal interactions will occur, which will guide the reduction basis selection process. Equivalently, this may serve as a tool for verifying the accuracy of ROMs without the need for full-order simulations. The proposed method is demonstrated using a simple oscillator and a finite element model of a clamped-clamped beam.

Introduction
Engineering structures vibrating at large amplitudes experience geometric nonlinearity, which couples the underlying linear normal modes of the model via nonlinear stiffness functions in the equations of motion. These couplings foster energy exchange between modes and can give rise to complex dynamic behaviours, drawback of the ICE method is the fact that the computed ROM parameters vary with respect to the scaling of the forces used to obtain the static solution dataset, resulting in ROMs that lack robustness [29]. In addition, the ICE method assumes that the inertia of the condensed modes is negligible; this assumption breaks down when considering structures that can undergo large in-plane displacements, such as cantilever beams, and can lead to significantly inaccurate results. Recent works have addressed these limitations. Specifically, it has been shown that the quasi-static coupling between the reduced and condensed modes generates higher-order nonlinear terms in the reduced dynamics. When these are taken into account, the ROM parameter estimation procedure becomes robust with respect to the scaling of the static solution dataset [33,34]. In addition, it has been shown how the effect of the inertia of the condensed modes can be accounted for in the reduced dynamics; this additional treatment is referred to as inertial compensation (ICE(-IC)) [35]. This introduces some additional acceleration-and velocity-dependent terms in the reduced equations of motion, whose coefficients can be easily computed using the existing static solution dataset. Similar expressions in the reduced dynamics are seen in the direct model order reduction techniques such as the quadratic manifold with modal derivatives [36,37], as well as the more general concept of an invariant manifold based on the theory of normal forms [38][39][40][41].
Nevertheless, a major challenge associated with the ICE method, as well as reducedorder modelling techniques in general, remains: selecting the reduction basis without a priori knowledge of the full-order dynamics. Retaining redundant modes will lead to computationally suboptimal ROMs, while omitting dynamically significant modes will lead to inaccurate results, and important features such as internal resonances may not be captured. In other words, the ICE method relies on a slow/fast decomposition and is unable to capture any internal resonances between the reduced and condensed modes [42]. In this study, we aim to address this limitation and propose a method that can be used to predict the existence of internal resonances in conservative systems and thus guide the reduction basis selection process, without the need for full-order simulations. Specifically, we represent each condensed coordinate as the superposition of two components: one that is statically coupled to the reduced coordinates and another that is dynamically independent of them: the latter may be considered as the error associated with the static condensation of the mode in question. By using this framework, we show how these errors may be approximated during model reduction in a computationally efficient manner. This may serve as a tool for predicting internal resonances between the reduced and condensed coordinates, or, equivalently, for verifying the accuracy of ROMs by ensuring that static condensation approximation is sufficiently accurate for all operating conditions of interest.
This article is structured as follows. In §2, we explore the nature of quasi-static and dynamic modal coupling in geometrically nonlinear systems using a simple oscillator as a motivating example. In §3, we show how the error associated with static condensation may be approximated during model order reduction, which can be used to predict the existence of internal resonances. In § §4 and 5, we demonstrate the proposed method using the simple oscillator and an FE model of a clamped-clamped beam, respectively. Finally, conclusions are presented in §6.

Motivation
In this section, we explore the nature of modal coupling in general conservative, geometrically nonlinear systems. Specifically, we investigate the quasi-static coupling approximation, which is often used in reduced-order modelling frameworks, and its applicability in different scenarios. To this end, we consider a discrete 4-DOF system composed of two point masses, m, as a motivating example. The masses are free to move in the x-y plane and are connected to a fixed frame and to each other through a set of linearly elastic springs, with stiffness k i ∀i ∈ {1, 2, 3, 4, 5} and unstretched length 0 , as shown in figure 1. At the equilibrium position, all springs are undeformed and oriented either horizontally or vertically. This system may be considered an extension to the single-mass, 2-DOF oscillator previously studied in refs. [  Its equations of motion can be expressed in the form where x is the vector of physical displacements, M and K are the linear mass and stiffness matrices, respectively, and F x and f x are the vectors of external and nonlinear restoring forces, respectively. These are given by where d i (x) ∀i ∈ {1, 2, 3, 4, 5} are the lengths of the springs in the deformed configuration. The derivation of these can be found in appendix A. Here, we consider the dynamics of the system in its linear modal space, where the equations of motion are linearly uncoupled. The mode shape (φ n ) and natural frequency (ω n ) of the nth mode of the underlying linear system are computed by solving the eigenproblem (K − ω 2 n M)φ n = 0. After applying the transform x = Φq and pre-multiplying by Φ T , equation (2.1) can be rewritten as follows: where Φ is the matrix containing the mass-normalized mode shapes 1 in its columns, such that Φ T MΦ = I, Λ = Φ T KΦ is the diagonal matrix containing the corresponding squares of the natural frequencies along its leading diagonal, and F = Φ T F x and f(q) = Φ T f x (Φq) are, respectively, the vectors of external and nonlinear restoring forces in the modal space.  Note that the horizontal grounding springs (k 1 and k 3 ) are significantly stiffer than the vertical grounding ones, which creates a dichotomy between the natural frequencies of the first two modes (ω 1 = 1.0 rad s −1 and ω 2 = 4.6 rad s −1 ), and those of the other two modes (ω 3 = 31.6 rad s −1 and ω 4 = 31.9 rad s −1 ). This system aims to emulate, in a highly simplified manner, the slow/fast dynamics that characterize the low-frequency transverse modes and high-frequency in-plane modes in plate-or beam-like structures.
First, we consider the quasi-static behaviour of the system when only the first mode is forced directly. We numerically solve the static equations Λq + f(q) = F for q for a series of load cases where F 1 ∈ [−3, +3] and F n = 0 ∀n ∈ {2, 3, 4}, where F n denotes the nth element in F. Figures 2a,b show the quasi-static modal response of the system against the static force applied to the first mode and against the corresponding response of the first mode, respectively. In reduced-order modelling methods such as the ICE(-IC) [30,35], this dataset, or more commonly a subset thereof, is used to approximate the functions describing the nonlinear stiffness of the reduced mode(s) and sometimes the quasi-static relationship between the condensed modes and the reduced mode(s). Specifically, the latter task is carried out only for a small set of high-frequency in-plane (or membrane) modes, for which the inertial forces are assumed to be small relative to the internal restoring forces. Here, we generalize this approach by treating all condensed modes equally, irrespective of their natural frequency or characteristics of their mode shapes. We denote the quasi-static relationship between the nth mode and the reduced mode using the function g n and approximate it as the Kth-order polynomial, i.e.
where the coefficients B (n) k are identified via least-squares regression based on the static solution dataset, i.e. fitting to the curves shown in figure 2b. Figure 3a shows the first backbone curve of the 4-DOF oscillator, computed according to equations (2.1) and (2.2) (with F x = 0) using the MATLAB-based numerical continuation toolbox Continuation Core (COCO) [13]. The branch emerging near Ω = 1.55 rad s −1 corresponds to a 1:3 internal resonance between the first and second modes. Figure 3b shows the time history of the modal response of the system (solid black lines), for three different nonlinear normal modes (NNMs) associated with the first backbone curve, which are represented by red dots in figure 3a and correspond to fundamental response frequencies of 1.01, 1.39 and 1.53 rad s −1 , respectively. The quasi-static response of modes 2-4 is computed by evaluating the functions g n (q 1 ) during the NNM motion: this is represented by dash-dotted lines in figure 3b. The difference between the dynamic and the quasi-static response of each mode, i.e. q n − g n (q 1 ), is represented by dashed blue   lines. This may be considered as the error arising from the quasi-static approximation/implicit condensation. It can be seen that, as expected, the quasi-static approximation is sufficiently accurate when applied to the high-frequency modes, as q n ≈ g n (q 1 ), ∀n ∈ {3, 4}, ∀t, ∀Ω. However, in the case of the second, low-frequency mode, the quasi-static approximation is initially moderately accurate near the first linear natural frequency, but becomes increasingly inaccurate as the system approaches internal resonance. Interestingly, for all NNM solutions, the error arising from this approximation appears to be a single-harmonic signal of frequency 3Ω. This suggests that the response of each condensed mode may be naturally decomposed into two parts: a component that is quasi-statically coupled to the reduced mode(s) and a component that is dynamically independent of the reduced mode(s). In the next section, we exploit this idea to show how the existence of dynamic interactions can be predicted and discuss how this is fundamentally relevant to reduced-order modelling.

Predicting dynamic interactions (a) Derivation
We start by considering the semi-discretized equations of motion of a conservative, linearly elastic, geometrically nonlinear structure, which are typically obtained through a finite element procedure. As mentioned earlier, these can be written as a series of N second-order ordinary differential equations in the form andq in the physical and mass-normalized linear modal spaces, respectively, where N is the number of DOFs of the FE model. Following ref. [35], but included here in the condensed form for completeness, we now separate the modal coordinates into three distinct classes. The first class, denoted by the subscript • r (reduced), consists of a small set of modes, which, for a given set of operating conditions, contain the majority of the total energy in the system and are dynamically independent: these modes must form the basis of the ROM. The number of modes in this class, R N, dictates the lower limit to the number of DOFs that an accurate ROM must have. For a single backbone curve in the absence of any internal resonance, R = 1. The second class, denoted by the subscript • s (static), is composed of S N modes, which may contain a substantial fraction of the energy of the full-order system, yet their response can be approximated as being quasi-statically coupled to the reduced modes. These modes need not be included as independent DOFs in the reduction basis, as their effects can be incorporated implicitly in the reduced dynamics [30, 31,35]. Finally, the third class, denoted by the subscript • u (unmodelled), contains the remaining U = (N − R − S) modes, which are very weakly coupled to the reduced modes, such that they always contain a negligible amount of energy under the operating conditions of interest. As such, it is assumed that these modes can be ignored during the reduction process with negligible loss of accuracy. By using this framework, the modal equations of motion of the FE model, equation (3.1b), can be rewritten as follows: 2 Ignoring the third group of weakly coupled modes, these can be approximated as follows: wheref r (r, s) := f r (r, s, 0) andf s (r, s) := f s (r, s, 0), such that q r ≈ r, q s ≈ s, q u ≈ u = 0 and x ≈ Φ r r + Φ s s. Equivalently, the kinetic energy of the full-order system, T , can be approximated as follows: As discussed in ref. [35], when the modes in the second group, s, can be expressed as functions of the reduced modes, r, i.e. s = g(r), then equation (3.3) can be exactly reduced tö where g(r) is an S × 1 vector of quasi-static coupling functions, ∂g/∂r is its S × R Jacobian matrix, ∂ 2 g/∂r 2 is its S × R × R second derivative tensor andf r (r) :=f r (r, g). Reduced-order models based on equation (3.6) were found to produce remarkably accurate results for systems where a clear slow/fast dynamic behaviour can be observed, e.g. between the low-frequency transverse modes and the highly stiff in-plane modes in thin plates and slender beams.
Here, we aim to broaden the scope of the implicit condensation approach and seek to quantify the error introduced by the static condensation. We define this error using the S × 1 time-dependent vector h(t), i.e. s = g(r) + h. (3.7) Using equation (3.7), and noting thatṡ = (∂g/∂r)ṙ +ḣ, the Lagrangian of the system, L, can be expressed as follows: as shown in appendix B. By using a Taylor series expansion about s = g, the stiffness expressions in equation (3.9) can be approximated as follows: and As expected, it can be seen that, when h = 0, equation (3.12a) is equivalent to the reduced dynamics obtained for the perfectly statically coupled case (equation (3.6)). The additional h-andḧ-dependent terms that appear on the right-hand side of equation (3.12a) may be considered as a forcing arising due to the dynamic coupling between r and s. In addition, when the kinetic energy of the statically coupled modes is negligibly small (i.e.ṡ Tṡ ≈ 0), the g-dependent terms in equation (3.6) are negligible, and the reduced dynamics can be simply expressed as follows:r + Λ r r +f r (r) = F r . (3.14) As discussed in ref. [35], this more traditional form of the ROM equations of motion is suitable for structures in which in-plane displacements are limited. The remainder of this article demonstrates how the equations of motion for h (equation (3.12b)) may be used to efficiently predict the presence of dynamic coupling between the reduced modes (r) and the condensed modes (s). This allows features such as internal resonances to be predicted and ROMs to be validated without the need for full-order FE simulations.

(b) Use in reduced-order modelling frameworks
The linear properties of the ROM, Λ r , Φ r and Φ s , can be computed directly using the linear mass and stiffness matrices of the FE model. The reduced nonlinear stiffness functions,f r (r), and the quasi-static coupling functions, g(r), are approximated indirectly in a least-squares manner using a force-displacement dataset for a series of nonlinear static solutions extracted from the FE model, as detailed in refs. [29,35]. These are approximated as the Kth-order polynomials of the reduced coordinates. Once the ROM parameters are identified, the reduced backbone curves can be computed, e.g. using numerical continuation, based on either equation (3.6) or equation (3.14).
Using equation (3.12b), additional insight can now be gained by simulating the error dynamics for each NNM of the ROM. Given that the static modal coupling is well captured through the functions g(r), then, in the absence or near the onset of a dynamic interaction, the error h is expected to be relatively small such that any nonlinear monomials of h become negligible, i.e. O(h 2 ) ≈ 0. As such, the linearized version of equation (3.12b) may be considered, which can be solved efficiently using, for example, the harmonic balance method [1]. To this end, the S × S matrix containing the linear coefficients of h must be approximated as functions of r. As with f r (r) and g(r), these can be computed based on least-squares polynomial regression. In this case, the tangent stiffness matrix, K tan , must be extracted for each nonlinear static load case, in addition to the vectors of applied forces and resulting displacement, which are needed for the standard ICE(-IC) method. The coefficients of each function in C(r) are then computed using the {r, Φ T s K tan Φ s } dataset evaluated at each load case. The practical implications of this are considered later in §5, where the proposed techniques are demonstrated using the commercial FE software Abaqus.
Once these functions are approximated, then the matrix C, as well as the vector on the righthand side of equation (3.12b), p, can be evaluated during a periodic solution of the ROM. These can then be expressed as summations of sinusoidal components, i.e. and where N h is the number of harmonics, and a (ij) n , α  where c h and c p are vectors containing the harmonic coefficients of h and p, respectively, and Υ is a matrix that can be algorithmically populated with the harmonic coefficients of C. Further details are given in appendix D. Finally, the computed amplitudes in c h are used to estimate the time history of the error corresponding to each condensed mode during a periodic motion of the ROM (equation (3.16a)). As discussed in §1, the main limitation of the ICE(-IC) method remains the fact that it relies on a slow/fast decomposition between the reduced and condensed modes. With the method proposed herein, this limitation is overcome, as the condensed basis can now be formed using any modes, irrespective of their natural frequency. This method can be used to monitor the error associated with the static condensation and thus identify if/when a condensed mode becomes resonant. This would suggest the quasi-static coupling assumption is no longer appropriate for the mode in question, and instead the mode must be included as an independent DOF in the reduction basis (Φ r ). Equivalently, this method can be used as a tool for validating ROMs, by ensuring that the component of the condensed modes that is dynamically independent of the reduced modes (i.e. h) remains sufficiently small for all operating conditions of interest.

Application to the 4-DOF oscillator (a) Single-mode ROM results
We now revisit the 4-DOF oscillator considered in §2 to demonstrate our proposed method. We first compute a seventh-order (K = 7), single-DOF ROM of the first mode using the ICE(-IC) method (equation (3.6)), while the remaining three modes are included in the statically condensed basis, i.e. Φ r = [φ 1 ], Φ s = [φ 2 φ 3 φ 4 ], and Φ u is unpopulated as no modes are unmodelled. The dataset used to calibrate the ROM consists of a series of static solutions where the the static force applied to the first mode is equally spaced between −3 and +3. An example of the fitting procedure for the reduced nonlinear stiffness function (f r (r)) and the quasi-static coupling functions (g(r)) is shown in figures 4a,b, respectively.   The backbone curve of the ROM when the static condensation error is taken into account, using our proposed method, is represented by red lines. The first backbone curve of the full-order system is also shown for reference (dashed grey lines). (Online version in colour.) Figure 5 shows the backbone curve of the computed single-DOF ROM, in the projection of modal amplitudes against fundamental response frequency (blue lines). Note that modes 2-4 are not modelled directly, but their response is approximated using the quasi-static coupling functions. It can be seen that the primary response of the reduced mode, and that of the highfrequency condensed modes, can be accurately predicted by the ROM. As expected, however, the internal resonance near Ω = 1.6 rad s −1 cannot be captured. This is due to the dynamic energy transfer between modes that takes place during an internal resonance, an effect that the static condensation approach is unable to capture.

(b) Internal resonance prediction
By using equation (3.12b) and the harmonic balance method, as described in §3b and appendix D, the error associated with the static condensation of each mode can now be estimated. To achieve this, the elements of the tangent stiffness matrix corresponding to the condensed modes, i.e. Λ tan,s = Φ T s K tan Φ s , must be approximated as functions of the reduced coordinates. 3 Similar to the approximation of the reduced nonlinear stiffness functions and the quasi-static coupling functions, these are computed in a least-squares manner, based on the same set of full-order nonlinear static solutions. Examples of this are shown in figures 4c,d.
The improved prediction of the reduced backbone curves, which takes into account the approximated error arising from the static condensation (with N h = 7), is now represented by red lines in figure 5. From this, it can be seen that, as the response frequency increases, the magnitude of h 2 relative to g 2 becomes increasingly large until Ω ≈ 1.6 rad s −1 , before it rapidly decreases again. 4 This singularity-type behaviour is caused by the third harmonic component of h 2 , and it indicates that a dynamic interaction between the first and second modes exists, without the need for simulating the dynamics of both modes simultaneously. This suggests that, in this region, a single-mode ROM is no longer appropriate, and the second mode must be included in the reduction basis. Note that, when considering FE models using commercial packages, the backbone curves of the full-order model cannot be readily computed. As such, significant features of the system, such as internal resonances, can be often overlooked during model reduction. With our proposed method, the existence of such dynamic interactions can be predicted, without the need for expensive full-order simulations.

(c) Two-mode ROM results
We now compute a two-DOF, seventh-order ROM (Φ r = [φ 1 φ 2 ], Φ s = [φ 3 φ 4 ]) using the same procedure. The dataset used to calibrate the ROM consists of a series of static load cases where one or both of the reduced modes are forced directly: the amplitude of the maximum force applied to either mode is |F 1 | = 3, as mentioned earlier, and |F 2 | = 8. The computed backbone curves are shown in figure 6. The ROM is now able to accurately capture the dynamic behaviour of the fullorder system for the whole range of frequencies considered. In this case, the additional treatment proposed herein acts as a tool for validating the ROM, as it indicates that the static condensation approximation is sufficiently accurate: this is determined by observing that the contribution from h is negligible for both condensed modes. In §5, we demonstrate the proposed method using a finite element model of a clamped-clamped beam.

Application to a finite element model
We now consider a geometrically nonlinear model of a clamped-clamped beam, constructed using the commercial FE software Abaqus, with the Abaqus2Matlab [45] toolbox used for preand post-processing. The beam model is identical to the one studied in ref. [35]. It has a length of 300 mm and a rectangular cross-section of 25 mm × 1 mm and is made of steel with a Young's modulus of 205 GPa, mass density of 7800 kg m −3 and Poisson's ratio 0.3. The mesh consists of 120 shear deformable, three-node quadratic beam elements (B32), resulting in 1434 DOFs.  We compute a quintic 5 single-DOF ROM of the first (bending) mode using the ICE method and include modes 3, 6, 72 and 129 in the condensed basis: these correspond to the second and third symmetric bending modes, and the first and second symmetric axial modes, respectively. 6 The static solution dataset used to calibrate the ROM consists of four load cases where the static force applied to the first mode is F 1 = {−45, −22.5, +22.5, +45}, resulting in a maximum static transverse deflection of 1.13 mm at the beam midspan. For this model, the computation time for each load case is about 25 seconds on a standard computer. As mentioned earlier, the fitting procedure for the reduced nonlinear stiffness functions, the quasi-static coupling functions and the tangent stiffness functions for the condensed modes 7 is carried out in a least-squares manner. Figure 7a shows the backbone curve of the single-DOF ROM (top), as well as the corresponding normalized amplitude of the error in the condensed modes (bottom). As the backbone curves of the full-order FE model cannot readily be computed, and thus cannot be directly compared with those of the ROM, we assess the accuracy of the ROM by comparing a set of reduced NNMs with the dynamic response of the FE model when the corresponding initial conditions are applied for a wide range of response frequencies. The initial conditions are enforced in the form of initial applied modal forces, as proposed in ref. [47]. The periodic orbits of the ROM are compared with the FE response in the time domain, as shown in figure 7b. From this, it can be seen that the ICE method can provide accurate results for the reduction of the clamped-clamped beam model, even with a single-mode ROM: this agrees with observations seen in the literature [20,29,30].
Interestingly, while the standard ROM results, as shown in figure 7b, indicate good agreement with the full-order model, the results obtained from the error-monitoring treatment suggest that there is a strong dynamic interaction between the first and third modes of the beam near Ω = 410 rad s −1 , which the single-mode ROM is unable to capture. Note that, while the h  components of modes 6, 72 and 129 are also large, the component of mode 3 is the largest and begins its rapid growth at a lower frequency than other modes. As such, mode 3 is treated as the candidate for a dynamic interaction. Specifically, this is associated with the amplitude of the fifth harmonic of the third mode. This result points to the existence of a 1:5 internal resonance between the first and third modes, a feature of flat clamped-clamped beams that has been widely observed in the literature, e.g. [19,20].
To verify this observation, we compute a quintic two-DOF ICE ROM of the beam, with maximum force applied in each reduced mode is F 1 = 45 and F 3 = 360. The resulting backbone curve is shown in figure 8a. It can be seen that the ROM now exhibits a 1:5 internal resonance, while the error predicted using our proposed method remains relatively small. As mentioned earlier, the accuracy of the reduced internally resonant NNMs is verified by comparing them with the corresponding set of responses of the full-order model (figure 8b).
It should be noted that, in addition to enabling the internal resonance to be modelled, the two-DOF ROM also leads to more accurate response predictions on the primary backbone curve. We quantify the accuracy of each reduced NNM using the periodicity metric , as defined in ref. [48], i.e.
where x 0 is the displacement vector of the FE model at t = 0, which is imposed as an initial condition based on the reduced NNM, and x T is the displacement vector of the FE model after one period. A smaller value indicates a response, which is closer to being periodic and thus a more accurate ROM. The computed periodicity values for different NNMs on the primary backbone curves, both for the single-and two-DOF ROMs, are shown in figure 9. The results suggest that the third mode of the FE model becomes increasingly important for response frequencies higher than ∼ 420 rad s −1 . This is in qualitative agreement with the results shown in figure 7a and further confirms the validity of the error-approximation procedure.

Conclusion
In this article, we have shown how the existence of internal resonances may be predicted in a computationally efficient manner during model order reduction of geometrically nonlinear structures. Specifically, we have used a simple oscillator as a motivating example to show how each modal coordinate in a nonlinear system may be represented as the sum of a component that is quasi-statically coupled to a small number of coordinates, which must form the reduction basis in a ROM, and a component that is dynamically independent of them. The latter part may be considered as the error arising from static condensation, which is the concept on which methods such as the implicit condensation and expansion rely, and is typically applied to structures characterized by slow/fast dynamics. We have employed the harmonic balance method to approximate the error dynamics independently of the reduced dynamics, thus enabling any dynamic interaction between the reduced and condensed modes to be identified. This can be achieved in a very computationally efficient manner, as a linear approximation of the error dynamics is considered. We have demonstrated the proposed method using the simple oscillator, as well as a finite element model of a clamped-clamped beam, and shown how the existence of a 1:3 and a 1:5 internal resonance, respectively, could be predicted based on single-mode ROMs.