## Abstract

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.

### 1. 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, which linear systems cannot exhibit and hence cannot be studied using methods from the linear vibration theory.

Examples of such behaviours include jump phenomena [1], limit-cycle oscillations [2,3], chaotic vibrations [4] and internal resonances [5,6]. During an internally resonant motion, the energy transfer between modes is significantly enhanced. This effect can be exploited for practical applications such as energy harvesting using electromagnetic and piezoelectric devices [7,8] and for enhanced stable micromechanical oscillators [9,10]. As such, the ability to reliably identify such behaviours is of great importance.

Several analytical techniques have been developed for analysing nonlinear dynamic systems, which allow approximate closed-form solutions to be obtained directly; three well-established techniques are the harmonic balance, multiple scales and normal form methods [1,11,12]. Alternatively, the response of nonlinear systems may be computed efficiently and accurately using numerical continuation packages, without the need for the analytical treatment [13,14]. However, complex engineering structures are often modelled through a finite element (FE) discretization, for which the application of such methods is infeasible. This is both due to the fact that the discretized equations of motion are not readily accessible when using commercial FE software and, more importantly, because the number of degrees-of-freedom (DOFs) of FE models is often far too large for such methods to be applied directly. Reduced-order modelling enables nonlinear FE models to be analysed using the aforementioned analytical and numerical methods by capturing the salient dynamic behaviour of the system in a highly efficient manner.

A vast body of literature is devoted to so-called *indirect* reduced-order modelling methods [15], which are the focus of this study. These methods can be used in conjunction with any commercial FE software package, as they do not rely on access to the FE code or knowledge of the full-order equations of motion. Instead, these methods make use of some standard outputs from the FE model to construct parametric reduced-order models (ROMs). Due to their non-intrusive nature, they are well suited for industrial applications and complex structures whose equations of motions cannot be derived analytically.

There exist two main families of indirect reduced-order modelling methods: displacement based and force based. In the former one, the dataset used to calibrate the ROM is obtained by constraining the structure into a series of prescribed displacements and extracting the corresponding internal restoring forces required to maintain each static configuration. This technique is often referred to as the enforced displacement procedure or stiffness evaluation procedure [16–18]. Using the static solution dataset, the coefficients of the nonlinear stiffness terms are evaluated for a subset of the modal coordinates of the full-order model, which form the reduction basis. The main drawback of this approach is its slow convergence, requiring that a large number of modes is included in the reduction basis and leading to large and computationally expensive ROMs [15,19–22]. Specifically, in addition to the low-frequency, dynamically important modes, a set of high-frequency in-plane modes must be retained to capture the effect of membrane stretching that occurs in flat structures such as beams and plates. Alternatively, the reduction basis may be augmented with the so-called dual or companion modes, but their identification is cumbersome [23–26]. Similarly, methods that augment the reduction basis with the so-called modal derivatives have recently been proposed [27,28].

In the second family of indirect methods, the dataset used for computing the nonlinear stiffness coefficients is obtained by applying a set of static forces to the structure, which are proportional to the reduced modeshapes, and extracting the corresponding displacement—this technique is referred to as the applied loads procedure or implicit condensation and expansion (ICE) [29–32]. This approach relies on static condensation to account for the effect of membrane stretching and thus removes the need for including any high-frequency in-plane modes in the reduction basis as independent DOFs. Instead, these are assumed to be quasi-statically coupled to the reduced modes, and their effect can be captured implicitly. The ICE method has been successfully applied for the reduction of a range of flat structures, see, for example, refs. [20,33]; however, it has traditionally suffered from issues related to the lack of invariance of the reduced subspace. A main 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–41].

Nevertheless, a major challenge associated with the ICE method, as well as reduced-order 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.

### 2. 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}\hspace{0.17em}\mathrm{\forall}i\in \{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. [34,43,44].

Its equations of motion can be expressed in the form

**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

*a*

*b*

*c*

**ϕ**

_{n}) and natural frequency (

*ω*

_{n}) of the

*n*th mode of the underlying linear system are computed by solving the eigenproblem $(\mathbf{K}\mathbf{-}{\omega}_{\mathbf{n}}^{\mathbf{2}}\mathbf{M}\mathbf{)}{\mathit{\varphi}}_{\mathbf{n}}\mathbf{=}\mathbf{0}$. After applying the transform $\mathbf{x}\mathbf{=}\mathit{\Phi}\mathbf{q}$ and pre-multiplying by ${\mathit{\Phi}}^{\text{T}}$, equation (2.1) can be rewritten as follows:

^{1}in its columns, such that ${\mathit{\Phi}}^{\text{T}}\mathbf{M}\mathit{\Phi}\mathbf{=}\mathbf{I}$, $\mathit{\Lambda}={\mathit{\Phi}}^{\text{T}}\mathbf{K}\mathit{\Phi}$ is the diagonal matrix containing the corresponding squares of the natural frequencies along its leading diagonal, and $\mathbf{F}\mathbf{=}{\mathit{\Phi}}^{\text{T}}{\mathbf{F}}_{\mathbf{x}}$ and $\mathbf{f}\mathbf{(}\mathbf{q}\mathbf{)}\mathbf{=}{\mathit{\Phi}}^{\text{T}}{\mathbf{f}}_{\mathbf{x}}\mathbf{(}\mathit{\Phi}\mathbf{q}\mathbf{)}$ are, respectively, the vectors of external and nonlinear restoring forces in the modal space.

For this motivating example, the physical parameters of the system are set to the following values: *m* = 1 kg, ℓ_{0} = 1 m, *k*_{1} = 1000 N m^{−1}, *k*_{2} = 10 N m^{−1}, *k*_{3} = 1000 N m^{−1}, *k*_{4} = 1 N m^{−1} and *k*_{5} = 21.2 N m^{−1}. 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 $\mathit{\Lambda}\mathbf{q}\mathbf{+}\mathbf{f}\mathbf{(}\mathbf{q}\mathbf{)}\mathbf{=}\mathbf{F}$ for **q** for a series of load cases where *F*_{1} ∈ [ − 3, + 3] and ${F}_{n}=0\hspace{0.17em}\mathrm{\forall}n\in \{2,3,4\}$, where *F*_{n} denotes the *n*th element in **F**. Figures 2*a*,*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 *n*th mode and the reduced mode using the function *g*_{n} and approximate it as the *K*th-order polynomial, i.e.

*b*.

Figure 3*a* 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 $\Omega =1.55\hspace{0.17em}{\text{rad s}}^{-1}$ corresponds to a 1:3 internal resonance between the first and second modes. Figure 3*b* 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 3*a* 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 3*b*. 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}\approx {g}_{n}({q}_{1}),\hspace{0.17em}\mathrm{\forall}n\in \{3,4\},\hspace{0.17em}\mathrm{\forall}t,\hspace{0.17em}\mathrm{\forall}\Omega $. 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\Omega $. 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.

### 3. 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

*a*

*b*

*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 ${\bullet}_{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 ${\bullet}_{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 ${\bullet}_{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}

**q**

_{r}≈

**r**,

**q**

_{s}≈

**s**,

**q**

_{u}≈

**u**=

**0**and $\mathbf{x}\approx {\mathit{\Phi}}_{\mathbf{r}}\mathbf{r}\mathbf{+}{\mathit{\Phi}}_{\mathbf{s}}\mathbf{s}$. Equivalently, the kinetic energy of the full-order system, $\mathcal{T}$, can be approximated as follows:

*a*

*b*

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 to

**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 and ${\stackrel{~}{\mathbf{f}}}_{r}(\mathbf{r}\mathbf{)}\mathbf{:=}{\hat{\mathbf{f}}}_{\mathbf{r}}\mathbf{(}\mathbf{r}\mathbf{,}\mathbf{g}\mathbf{)}$. 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.

**r**and

**h**can be derived using the Euler–Lagrange equation. This leads to

*a*

*b*

By using a Taylor series expansion about **s** = **g**, the stiffness expressions in equation (3.9) can be approximated as follows:

*a*

*b*

*a*

*b*

**g**(

**r**) is computed based on the

*static*response of the system, with a static force only applied to the reduced modes (i.e.

**F**

_{s}=

**0**). Substituting equation (3.10) into equation (3.9), and noting that

**B**+ (∂

**g**/∂

**r**)

^{T}

**C**=

**0**as shown in appendix C, leads to

*a*

*b*

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 $\ddot{\mathbf{h}}$-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. ${\dot{\mathbf{s}}}^{\text{T}}\dot{\mathbf{s}}\approx 0$), the **g**-dependent terms in equation (3.6) are negligible, and the reduced dynamics can be simply expressed as follows:

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, ${\mathit{\Lambda}}_{r}$, ${\mathit{\Phi}}_{r}$ and ${\mathit{\Phi}}_{s}$, can be computed directly using the linear mass and stiffness matrices of the FE model. The reduced nonlinear stiffness functions, ${\stackrel{~}{\mathbf{f}}}_{r}(\mathbf{r}\mathbf{)}$, 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 *K*th-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. $\mathcal{O}({\mathbf{h}}^{\mathbf{2}}\mathbf{)}\approx \mathbf{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 ${\stackrel{~}{\mathbf{f}}}_{r}(\mathbf{r}\mathbf{)}$ 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 $\mathbf{C}\mathbf{(}\mathbf{r}\mathbf{)}$ are then computed using the $\{\mathbf{r}\mathbf{,}{\mathit{\Phi}}_{\mathbf{s}}^{\text{T}}{\mathbf{K}}_{\text{tan}}{\mathit{\Phi}}_{\mathbf{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 right-hand 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.

*a*

*b*

*N*

_{h}is the number of harmonics, and ${{a}^{(ij)}}_{n}$, ${\alpha}_{n}^{(i)}\hspace{0.17em}\mathrm{\forall}n\in [0,{N}_{h}]\cap \mathbb{Z}$, and ${{b}^{(ij)}}_{n}$, ${\beta}_{n}^{(i)}$ $\hspace{0.17em}\mathrm{\forall}n\in [1,{N}_{h}]\cap \mathbb{Z}$, are coefficients, which can be identified via a discrete Fourier transform for each element in

**C**and

**p**, during an NNM motion with fundamental response frequency

*ω*. Similarly, each element in

**h**can be expressed as a sum of its Fourier components, i.e.

*a*

*b*

**h**computed by solving the resulting set of simultaneous linear equations. This can be expressed as follows:

**c**

_{h}and

**c**are vectors containing the harmonic coefficients of

_{p}**h**and

**p**, respectively, and $\mathbf{{\rm Y}}$ 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 (${\mathit{\Phi}}_{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.

### 4. 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. ${\mathit{\Phi}}_{r}=[{\mathit{\varphi}}_{1}]$, ${\mathit{\Phi}}_{s}=[{\mathit{\varphi}}_{2}\hspace{0.17em}{\mathit{\varphi}}_{3}\hspace{0.17em}{\mathit{\varphi}}_{4}]$, and ${\mathit{\Phi}}_{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 (${\stackrel{~}{\mathbf{f}}}_{r}(\mathbf{r}$)) and the quasi-static coupling functions (**g**(**r**)) is shown in figures 4*a*,*b*, respectively.

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 high-frequency condensed modes, can be accurately predicted by the ROM. As expected, however, the internal resonance near $\Omega =1.6\hspace{0.17em}{\text{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. ${\mathit{\Lambda}}_{\text{tan},s}={\mathit{\Phi}}_{s}^{\text{T}}{\mathbf{K}}_{\text{tan}}{\mathit{\Phi}}_{\mathbf{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 4*c*,*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 $\Omega \approx 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 (${\mathit{\Phi}}_{r}=[{\mathit{\varphi}}_{1}\hspace{0.17em}{\mathit{\varphi}}_{2}]$, ${\mathit{\Phi}}_{s}=[{\mathit{\varphi}}_{3}\hspace{0.17em}{\mathit{\varphi}}_{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 full-order 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.

### 5. 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 pre- and 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 7*a* 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 7*b*. 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 7*b*, 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 $\Omega =410\hspace{0.17em}{\text{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 ${\mathit{\Phi}}_{r}=[{\mathit{\varphi}}_{1}\hspace{0.17em}{\mathit{\varphi}}_{3}]$ and ${\mathit{\Phi}}_{s}=[{\mathit{\varphi}}_{6}\hspace{0.17em}{\mathit{\varphi}}_{72}\hspace{0.17em}{\mathit{\varphi}}_{129}]$. The calibration dataset consists of 24 load cases, where the maximum force applied in each reduced mode is *F*_{1} = 45 and *F*_{3} = 360. The resulting backbone curve is shown in figure 8*a*. 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 8*b*).

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.

**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 7

*a*and further confirms the validity of the error-approximation procedure.

### 6. 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.

The significance of this development is twofold. First, the method presented herein enables the identification of the modes, which must form the reduction basis, without relying on knowledge of the full-order dynamics. Second, this method can serve as a computationally cheap validation of ROMs without the need for full-order simulations, which can sometimes be infeasible to obtain. This removes the need for cumbersome trial-and-error processes and case-by-case treatment guided by empirical rules and is a key step towards developing reduced-order modelling methods which can be applied systematically to a broad range of structures.

### Data accessibility

FE model submitted as electronic supplementary material.

### Authors' contributions

E.N.: the development of the work, with supervisory support from T.L.H. and S.A.N. on the development of the idea. All authors contributed to the preparation of the manuscript.

### Competing interests

We declare we have no competing interests.

### Funding

E.N. is supported by an EPSRC DTP studentship and S.A.N. is supported by an EPSRC Programme (grant no. EP/R006768/1).

## Appendix A. Derivation of the equations of motion of the 4-DOF oscillator

The potential energy of spring *i* is given by

*d*

_{i}is its length in the deformed configuration, as defined in equation (2.2c). The Lagrangian of the system can be expressed as follows:

## Appendix B. Derivation of the reduced and error dynamics

Using equation (3.8), the partial derivatives of $\mathcal{L}$ with respect to $\dot{\mathbf{r}}$, $\dot{\mathbf{h}}$, **r** and **h** can be written, respectively, as follows:

*a*

*b*

*c*

*d*

*a*) and (B 1

*b*) with respect to time gives

*a*

*b*

**r**and

**h**can be derived using the Euler–Lagrange equation, i.e.

*a*

*b*

## Appendix C. Relationship between B and C

Starting from equation (3.3), it can be seen that the quasi-static coupling functions, **g**(**r**), are defined such that the following equation is satisfied:

**r**leads to

**B**+ (∂

**g**/∂

**r**)

^{T}

**C**=

**0**, since

**C**=

**C**

^{T}.

## Appendix D. Estimation of the harmonic coefficients of h

Using equations (3.15a) and (3.16), and neglecting harmonics higher than *N*_{h}*ω*, the left-hand side of equation (3.12b) can be expressed as follows:

*a*

*b*

*c*

*S*(2

*N*

_{h}+ 1) × 1 vectors

**c**

_{h}and

**c**

_{p}contain the harmonic amplitudes of

**h**and

**p**, respectively, i.e.

*a*

*b*

*c*

*d*

*S*(2

*N*

_{h}+ 1) ×

*S*(2

*N*

_{h}+ 1) matrix $\mathbf{{\rm Y}}({{a}^{(ij)}}_{n},{{b}^{(ij)}}_{n})$ is populated according to equation (12). Finally, the harmonic coefficients of

**h**are computed by inverting equation (13), i.e. ${\mathbf{c}}_{\mathbf{h}}\mathbf{=}{\mathbf{{\rm Y}}}^{\mathbf{-}\mathbf{1}}{\mathbf{c}}_{\mathbf{p}}$.

## Footnotes

1 Note that mode shapes have the same units as the physical displacement vector, such that all modal quantities are dimensionless.

2 Note that, by definition, only the reduced modes have static forces applied during the model reduction procedure, i.e. **F**_{s} = **0** and **F**_{u} = **0**.

3 Note that $\mathit{\Lambda}:={\mathit{\Lambda}}_{\text{tan}}(\mathbf{q}\mathbf{=}\mathbf{0}\mathbf{)}$.

4 Note that, in the region where |*h*_{2}| is relatively large, the smallness assumption is violated, resulting in mispredictions of the static condensation errors in modes 3 and 4. Nevertheless, the early growth of *h*_{2} indicates the onset of a dynamic interaction, whilst the exact modal vibration amplitudes are of little importance.

5 Note that a quintic ROM is adopted, as it was found to be robust with respect to the scaling of the static forces used to calibrate it. This suggests that, for the response range considered here, a higher truncation order is not necessary.

6 Here, the condensed basis was chosen based on the relative amplitude of each mode evaluated at the static solutions, as discussed in ref. [35]. For more complex structures, more intricate methods of selecting good candidate modes for the condensed basis may be necessary—this remains a topic for future investigation. The reader is referred to ref. [46] for a recent work on mode-selection criteria based on the theory of spectral submanifolds.

### References

- 1.
- 2.
Patil MJ, Hodges DH, Cesnik CE . 2001 Limit-cycle oscillations in high-aspect-ratio wings.**J. Fluids Struct.**, 107-132. (doi:10.1006/jfls.2000.0329) Crossref, Web of Science, Google Scholar**15** - 3.
Thothadri M, Moon F . 2005 Nonlinear system identification of systems with periodic limit-cycle response.**Nonlinear Dyn.**, 63-77. (doi:10.1007/s11071-005-1914-0) Crossref, Web of Science, Google Scholar**39** - 4.
Wiggins S . 2003**Introduction to applied nonlinear dynamical systems and chaos**,. New York, NY: Springer Science & Business Media. Google Scholar**vol. 2** - 5.
Amabili M, Pellicano F, Vakakis A . 2000 Nonlinear vibrations and multiple resonances of fluid-filled, circular shells, part 1: equations of motion and numerical results.**J. Vib. Acoust.**, 346-354. (doi:10.1115/1.1288593) Crossref, Web of Science, Google Scholar**122** - 6.
Vakakis AF, Manevitch LI, Mikhlin YV, Pilipchuk VN, Zevin AA . 2001**Normal modes and localization in nonlinear systems**. New York, NY: Springer. Crossref, Google Scholar - 7.
Chen L-Q, Jiang W-A . 2015 Internal resonance energy harvesting.**J. Appl. Mech.**, 031004. (doi:10.1115/1.4029606) Crossref, Web of Science, Google Scholar**82** - 8.
Lan C, Qin W, Deng W . 2015 Energy harvesting by dynamic unstability and internal resonance for piezoelectric beam.**Appl. Phys. Lett.**, 093902. (doi:10.1063/1.4930073) Crossref, Web of Science, Google Scholar**107** - 9.
Antonio D, Zanette DH, López D . 2012 Frequency stabilization in nonlinear micromechanical oscillators.**Nat. Commun.**, 1-6. (doi:10.1038/ncomms1813) Crossref, Web of Science, Google Scholar**3** - 10.
Chen C, Zanette DH, Czaplewski DA, Shaw S, López D . 2017 Direct observation of coherent energy transfer in nonlinear micromechanical oscillators.**Nat. Commun.**, 1-7. (doi:10.1038/s41467-016-0009-6) PubMed, Web of Science, Google Scholar**8** - 11.
Jezequel L, Lamarque C-H . 1991 Analysis of non-linear dynamical systems by the normal form theory.**J. Sound Vib.**, 429-459. (doi:10.1016/0022-460X(91)90446-Q) Crossref, Web of Science, Google Scholar**149** - 12.
Nayfeh AH . 1981**Introduction to perturbation techniques**. New York, NY: Wiley-Interscience. Google Scholar - 13.
Dankowicz H, Schilder F . 2013**Recipes for continuation**,. Philadelphia, PA: SIAM. Crossref, Google Scholar**vol. 11** - 14.
Doedel EJ, Champneys AR, Fairgrieve TF, Kuznetsov YA, Oldeman BE, Paffenorth RC, Sandstede B, Wang X, Zhang C . 2007 Auto-07p: continuation and bifurcation software for ordinary differential equations. Montreal, Canada: Concordia University. Google Scholar - 15.
Mignolet MP, Przekop A, Rizzi SA, Spottswood SM . 2013 A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures.**J. Sound Vib.**, 2437-2460. (doi:10.1016/j.jsv.2012.10.017) Crossref, Web of Science, Google Scholar**332** - 16.
Muravyov AA, Rizzi SA . 2003 Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures.**Comput. Struct.**, 1513-1523. (doi:10.1016/S0045-7949(03)00145-7) Crossref, Web of Science, Google Scholar**81** - 17.
Perez R, Wang X, Mignolet MP . 2014 Nonintrusive structural dynamic reduced order modeling for large deformations: enhancements for complex structures.**J. Comput. Nonlinear Dyn.**, 031008. Crossref, Web of Science, Google Scholar**9** - 18.
Rizzi SA, Przekop A . 2005 The effect of basis selection on static and random acoustic response prediction using a nonlinear modal simulation. Technical report NASA/TP-2005-213943. Google Scholar - 19.
Givois A, Grolet A, Thomas O, Deü J-F . 2019 On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models.**Nonlinear Dyn.**, 1747-1781. (doi:10.1007/s11071-019-05021-6) Crossref, Web of Science, Google Scholar**97** - 20.
Kuether RJ, Deaner BJ, Hollkamp JJ, Allen MS . 2015 Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes.**AIAA J.**, 3273-3285. (doi:10.2514/1.J053838) Crossref, Google Scholar**53** - 21.
Rizzi SA, Przekop A . 2008 System identification-guided basis selection for reduced-order nonlinear response analysis.**J. Sound Vib.**, 467-485. (doi:10.1016/j.jsv.2007.12.031) Crossref, Web of Science, Google Scholar**315** - 22.
Vizzaccaro A, Givois A, Longobardi P, Shen Y, Deü J-F, Salles L, Touzé C, Thomas O . 2020 Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements.**Comput. Mech.**, 1293-1319. (doi:10.1007/s00466-020-01902-5) Crossref, Web of Science, Google Scholar**66** - 23.
Hollkamp J, Gordon R, Spottswood S . 2003 Nonlinear sonic fatigue response prediction from finite element modal models: a comparison with experiments,*44th AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics, and Materials Conference*, p. 1709. Google Scholar - 24.
Hollkamp JJ, Gordon RW, Spottswood SM . 2005 Nonlinear modal models for sonic fatigue response prediction: a comparison of methods.**Journal of Sound and Vibration**, 1145-1163. (doi:10.1016/j.jsv.2004.08.036) Crossref, Web of Science, Google Scholar**284** - 25.
Kim K, Radu AG, Wang X, Mignolet MP . 2013 Nonlinear reduced order modeling of isotropic and functionally graded plates.**Int. J. Non-Linear Mech.**, 100-110. (doi:10.1016/j.ijnonlinmec.2012.07.008) Crossref, Web of Science, Google Scholar**49** - 26.
Wang X, Mignolet M, Eason T, Spottswood S . 2009 Nonlinear reduced order modeling of curved beams: a comparison of methods,*50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference*, p. 2433. Google Scholar - 27.
Mahdiabadi MK, Tiso P, Brandt A, Rixen DJ . 2021 A non-intrusive model-order reduction of geometrically nonlinear structural dynamics using modal derivatives.**Mech. Syst. Signal Process.**, 107126. (doi:10.1016/j.ymssp.2020.107126) Crossref, Web of Science, Google Scholar**147** - 28.
Sombroek CS, Tiso P, Renson L, Kerschen G . 2018 Numerical computation of nonlinear normal modes in a modal derivative subspace.**Comput. Struct.**, 34-46. (doi:10.1016/j.compstruc.2017.08.016) Crossref, Web of Science, Google Scholar**195** - 29.
Gordon R, Hollkamp J . 2011 Reduced-order models for acoustic response prediction,*Technical Report AFRL-RB-WP-TR-2011-3040*, Air Force Research Laboratory, Dayton, OH. Google Scholar - 30.
Hollkamp JJ, Gordon RW . 2008 Reduced-order models for nonlinear response prediction: Implicit condensation and expansion.**J. Sound Vib.**, 1139-1153. (doi:10.1016/j.jsv.2008.04.035) Crossref, Web of Science, Google Scholar**318** - 31.
McEwan M, Wright JR, Cooper JE, Leung AYT . 2001 A combined modal/finite element analysis technique for the dynamic response of a non-linear beam to harmonic excitation.**J. Sound Vib.**, 601-624. (doi:10.1006/jsvi.2000.3434) Crossref, Web of Science, Google Scholar**243** - 32.
Segalman D, Dohrmann C . 1996 A method for calculating the dynamics of rotating flexible structures, part 1: derivation.**J. Vib. Acoust.**, 313-317. (doi:10.1115/1.2888183) Crossref, Web of Science, Google Scholar**118** - 33.
Frangi A, Gobat G . 2019 Reduced order modelling of the non-linear stiffness in mems resonators.**Int. J. Non-Linear Mech.**, 211-218. (doi:10.1016/j.ijnonlinmec.2019.07.002) Crossref, Web of Science, Google Scholar**116** - 34.
Nicolaidou E, Melanthuru VR, Hill TL, Neild SA . 2020 Accounting for quasi-static coupling in nonlinear dynamic reduced-order models.**J. Comput. Nonlinear Dyn.**, 071002. Crossref, Web of Science, Google Scholar**15** - 35.
Nicolaidou E, Hill TL, Neild SA . 2020 Indirect reduced-order modelling: using nonlinear manifolds to conserve kinetic energy.**Proc. R. Soc. A**, 20200589. (doi:10.1098/rspa.2020.0589) Link, Google Scholar**476** - 36.
Jain S, Tiso P, Rutzmoser JB, Rixen DJ . 2017 A quadratic manifold for model order reduction of nonlinear structural dynamics.**Comput. Struct.**, 80-94. (doi:10.1016/j.compstruc.2017.04.005) Crossref, Web of Science, Google Scholar**188** - 37.
Rutzmoser JB, Rixen DJ, Tiso P, Jain S . 2017 Generalization of quadratic manifolds for reduced order modeling of nonlinear structural dynamics.**Comput. Struct.**, 196-209. (doi:10.1016/j.compstruc.2017.06.003) Crossref, Web of Science, Google Scholar**192** - 38.
Pesheck E, Boivin N, Pierre C, Shaw SW . 2001 Nonlinear modal analysis of structural systems using multi-mode invariant manifolds.**Nonlinear Dyn.**, 183-205. (doi:10.1023/A:1012910918498) Crossref, Web of Science, Google Scholar**25** - 39.
Shaw SW, Pierre C . 1993 Normal modes for non-linear vibratory systems.**J. Sound Vib.**, 85-124. (doi:10.1006/jsvi.1993.1198) Crossref, Web of Science, Google Scholar**164** - 40.
Touzé C, Thomas O, Huberdeau A . 2004 Asymptotic non-linear normal modes for large-amplitude vibrations of continuous structures.**Comput. Struct.**, 2671-2682. (doi:10.1016/j.compstruc.2004.09.003) Crossref, Web of Science, Google Scholar**82** - 41.
Vizzaccaro A, Shen Y, Salles L, Blahoš J, Touzé C . 2020 Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. (http://arxiv.org/abs/2009.12145). Google Scholar - 42.
Shen Y, Béreux N, Frangi A, Touzé C . 2021 Reduced order models for geometrically nonlinear structures: assessment of implicit condensation in comparison with invariant manifold approach.**Eur. J. Mech.-A/Solids**, 104165. (doi:10.1016/j.euromechsol.2020.104165) Crossref, Web of Science, Google Scholar**86** - 43.
Touzé C, Amabili M . 2006 Nonlinear normal modes for damped geometrically nonlinear systems: application to reduced-order modelling of harmonically forced structures.**J. Sound Vib.**, 958-981. (doi:10.1016/j.jsv.2006.06.032) Crossref, Web of Science, Google Scholar**298** - 44.
Touzé C, Thomas O, Chaigne A . 2004 Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes.**J. Sound Vib.**, 77-101. (doi:10.1016/j.jsv.2003.04.005) Crossref, Web of Science, Google Scholar**273** - 45.
Papazafeiropoulos G, Muñiz-Calvente M, Martínez-Pañeda E . 2017 Abaqus2matlab: a suitable tool for finite element post-processing.**Adv. Eng. Softw.**, 9-16. (doi:10.1016/j.advengsoft.2017.01.006) Crossref, Web of Science, Google Scholar**105** - 46.
Buza G, Jain S, Haller G . 2021 Using spectral submanifolds for optimal mode selection in nonlinear model reduction.**Proc. R. Soc. A**, 20200725. (doi:10.1098/rspa.2020.0725) Link, Google Scholar**477** - 47.
Kuether RJ, Allen MS . 2014 A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models.**Mech. Syst. Signal Process.**, 1-15. (doi:10.1016/j.ymssp.2013.12.010) Crossref, Web of Science, Google Scholar**46** - 48.
VanDamme CI, Allen MS . 2017 Nonlinear normal modes of a curved beam and its response to random loading. In*Nonlinear Dynamics*, vol. 1, pp. 115–126. New York, NY: Springer. Google Scholar