Comparison of nonlinear mappings for reduced-order modelling of vibrating structures: normal form theory and quadratic manifold method with modal derivatives

Abstract

The objective of this contribution is to compare two methods proposed recently in order to build efficient reduced-order models for geometrically nonlinear structures. The first method relies on the normal form theory that allows one to obtain a nonlinear change of coordinates for expressing the reduced-order dynamics in an invariant-based span of the phase space. The second method is the modal derivative approach, and more specifically, the quadratic manifold defined in order to derive a second-order nonlinear change of coordinates. Both methods share a common point of view, willing to introduce a nonlinear mapping to better define a reduced-order model that could take more properly into account the nonlinear restoring forces. However, the calculation methods are different and the quadratic manifold approach has not the invariance property embedded in its definition. Modal derivatives and static modal derivatives are investigated, and their distinctive features in the treatment of the quadratic nonlinearity are underlined. Assuming a slow/fast decomposition allows understanding how the three methods tend to share equivalent properties. While they give proper estimations for flat symmetric structures having a specific shape of nonlinearities and a clear slow/fast decomposition between flexural and in-plane modes, the treatment of the quadratic nonlinearity makes the predictions different in the case of curved structures such as arches and shells. In the more general case, normal form approach appears preferable since it allows correct predictions of a number of important nonlinear features, including the hardening/softening behaviour, whatever the relationships between slave and master coordinates are.

1 Introduction

Reduced-order modelling of thin structures experiencing large-amplitude vibration is a topic that has attracted a large amount of researches in the last years. A number of methods have been proposed, with variants driven by either the structure under study and its peculiarity [63], the dynamical behaviour exhibited by the system [64], the model [54] or the discretisation method [33].

Roughly speaking, one can divide the techniques proposed in the literature into two different categories, the first one using linear change of coordinates and in the second family nonlinear mappings being defined. When referring to linear methods, one can also distinguish techniques where the best orthogonal basis is computed once and from all. Modal basis [3, 12, 34, 56], Ritz vectors [20], dual modes [25] and proper orthogonal decomposition (POD) [4, 23, 27] fall into that family. The proper generalised decomposition (PGD) [10, 14] under its progressive variant (pPGD) as defined in [32] also belongs to that case since additional vectors are added when the dynamics is becoming more complex. On the other hand, the linear change of coordinate can be adaptive, depending on the dynamics, the computation (single point or a whole branch of solution) or the location in phase space. Nonlinear principal component analysis (NLPCA) [22] as well as the optimised PGD (oPGD) [32] belongs to this family of improved linear methods, sometimes coined as nonlinear since the basis may change depending on some parameters.

In the third class of methods, a nonlinear change of coordinate is derived once and from all. Nonlinear normal modes [26, 37, 43, 46, 57], spectral submanifolds [16, 40] and the quadratic manifold derived from modal derivatives [19, 44] belong to this family. As shown in [5], when a linear method (e.g. POD) tries to find the best orthogonal axis fitting a learning set that has a complex shape, then the number of vectors will be larger than the number of curved subspaces one can use to describe the same data sets. In this particular example, it was shown that invariant manifolds pass exactly through the learning set, thus diminishing the number of coordinates needed to describe the dynamics.

Nonlinear normal modes (NNMs) and spectral submanifolds (SSM) offer a rigorously established conceptual framework for reducing geometrically nonlinear structures. In particular, the invariance property of reduction spaces is encapsulated in their definition, ensuring that the dynamical solutions computed from a reduced-order model (ROM) also exist for the full system [17, 45, 50, 52]. This key ingredient allows deriving accurate ROMs, which, for example, are able to predict the correct hardening/softening behaviour of nonlinear structure, which is not the case for their linear counterparts [57]. More specifically, recent contributions by Haller and collaborators have shown that SSMs are unique continuations of spectral subspaces of the linear system under the nonlinear terms [16] and are thus the best mathematical object to be used in the present context. For nonlinear conservative vibratory systems, SSMs simplify to the classic Lyapunov subcentre manifolds (LSM) that are filled with periodic orbits, thus unifying a number of definitions given for NNMs in the past decades, see, e.g. [21, 43, 45, 59].

On the other hand, modal derivatives (MDs) have been proposed independently [18, 64], and they share a number of common points with NNMs. In particular, MDs are defined by assuming that the mode shape (eigenvector) and its eigenfrequency have a dependence on amplitude, so that one can differentiate the classical Sturm–Liouville eigenvalue problem that defines linear normal modes, in order to make appear a quantity which is defined as the modal derivative. Symmetrically, NNMs also rely on the fact that modal quantities depend on amplitude, a key feature in nonlinear oscillations. The backbone curve and the dependence of the eigenmode shape with amplitude are then a result from the computation of NNMs, defined as invariant manifold in phase space. However, a complete comparison of both methods has not been drawn out yet. The only related paper uses the modal derivatives as a reduction method, from which the NNM, seen in this case as the family of periodic orbit in phase space—and thus reducing their information to the backbone curve only, without using the geometrical information in phase space—can be computed [49].

A recent development in the use of modal derivatives is to form a quadratic manifold for more accurate model order reduction. The properties of this nonlinear mapping are such that it is tangent to a subspace spanned by the most relevant vibration modes, and its curvature is provided by modal derivatives [19]. An idea also claimed in [44] is that such a quadratic manifold should be able to cancel the quadratic forces in the ROM. Incidentally, NNMs defined in the framework of normal form theory, as proposed in [53, 57], already present these features. Indeed, a third-order nonlinear change of coordinate is given, which has the property to be identity-tangent when the initial model is expressed in modal coordinates, thus conserving the linear modes as first approximation. Also, in case of no second-order internal resonance, the mapping exactly cancels all quadratic terms. Finally, the invariance property is directly inherited from the definition of an NNM as an invariant manifold in phase space, while the invariance of the quadratic manifold computed from MDs is not at hand.

The aim of this contribution is thus to investigate more properly the common points and differences in the two methods and explain their advantages and drawbacks in the context of building reduced-order models for geometrically nonlinear structures. The paper is organised as follows. Section 2 is concerned with the theoretical developments. The framework of geometrically nonlinear structures is briefly recalled; then, both methods of interest, normal form theory, modal derivatives and their extension to quadratic manifold (QM), are recalled and analysed in depth. The general derivation of the QM framework for both modal derivatives (MDs) and static modal derivatives (SMDs) is highlighted, whereas previous contributions generally use the simplifying assumption of SMDs in the developments. As a consequence of this development, the distinctive treatment of the quadratic nonlinearity between MDs and SMDs is specifically underlined. Of particular interest is the comparison of methods when a slow/fast decomposition of the system can be assumed. In the course of the paper, we will contrast the results given by MDs, SMDs and normal form and underline that the simplifying assumption of slow/fast approximation allows retrieving partly the correct results. By doing so, an illustration of the general theorem given in [17] is thus provided for a more restrictive framework. Indeed, theorems given in [17] encompass more generality and exact results, allowing to deal with the case of damping and forcing. We give, however, here more detailed comparisons and in particular analyse how the SMD can produce incorrect predictions for structures having a strong quadratic coupling such as arches and shells. Section 3 illustrates the findings of the previous section on two simple two-degree-of-freedom (dof) systems. Finally, Sect. 4 applies the previous results to continuous structures discretised with the finite-element (FE) procedure.

2 Models and methods

2.1 Framework

Geometric nonlinearity refers to the case of thin structures vibrating with large amplitudes, while the material behaviour remains linear elastic. In this framework, the semi-discretised version of the equations of motion, generally obtained from a finite-element procedure, reads :

$$\begin{aligned} \mathbf {M}\ddot{\mathbf {u}} + \mathbf {F}(\mathbf {u})=\mathbf {Q}, \end{aligned}$$

(1)

where $\mathbf {M}$ is the mass matrix, $\mathbf {u}$ the displacement vector at the nodes, $\mathbf {F}$ the nonlinear restoring force and $\mathbf {Q}$ the external force. The number of degrees of freedom (dofs) is N, being thus the dimension of vectors $\mathbf {u}$, $\mathbf {F}$ and $\mathbf {Q}$. Note that damping is presently not taken into account since most of the present work deals with efficient treatments of nonlinearities in the restoring force. While the concepts of NNMs and spectral submanifolds (SSM) can be straightforwardly extended to the cases with damping, as already shown, for example, in [53] for normal form or in [16] for SSM, a clear extension of MDs to damped systems does not seem to be present in the literature, to the best of our knowledge. Consequently, we restrict ourselves in this contribution to the treatment of the nonlinear stiffness without considering the effect of damping, but we acknowledge that damping has important effects in nonlinear vibrations that should thus need further investigations.

Geometric nonlinearity for slender structures is assumed so that $\mathbf {F}$, for the sake of simplicity, only depends on the displacement vector $\mathbf {u}$, but other cases can also be treated. More particularly, a number of models have been derived for thin structures such as plates and shells, relying on simplifying assumptions (e.g. von Kármán models for beams and plates [7, 28, 51], Donnell’s assumption for shallow shells [1, 2]), showing that the partial differential equations of motion only contain quadratic and cubic terms with respect to the displacement. On the other hand, general equations for three-dimensional elasticity with geometric nonlinearity (linear stress/strain relationship but nonlinear strain/displacement relationship) also show that the nonlinear terms in the restoring force should be of this type [12, 30, 33, 58]. Consequently, we consider in this contribution a nonlinear force that can be expressed as a function of the displacement up to cubic-order terms, reading:

$$\begin{aligned} \mathbf {F}(\mathbf {u})=\mathbf {K}\mathbf {u}+ \mathbf {G}\mathbf {u}\mathbf {u} + \mathbf {H}\mathbf {u}\mathbf {u}\mathbf {u}. \end{aligned}$$

(2)

In this equation, we use a simplified notation of the tensor product for the quadratic and cubic terms, already introduced in [19, 44]. The notation is fully explained in “Appendix A”, where the indicial expressions of the products are detailed for the sake of clarity. $\mathbf {G}$ is a third-order tensor of quadratic coefficients with current term $G^p_{ij}$, while $\mathbf {H}$ is the fourth-order tensor grouping the cubic coefficients $H^p_{ijk}$. For example, the vector $\mathbf {G}\mathbf {u}\mathbf {u}$ of the quadratic terms can be written as:

$$\begin{aligned} \mathbf {G}\mathbf {u}\mathbf {u}=\sum _{i=1}^N \sum _{j=1}^N \mathbf {G}_{ij}{u}_{i}{u}_{j}, \end{aligned}$$

(3)

with $\mathbf {G}_{ij}$ being the N-dimensional vector of coefficients $G^p_{ij}$, for $p=1,\ldots , N$. Note also that in this contribution, the representation of quadratic and cubic terms does not use the fact that the usual product is commutative ($u_{i}u_{j}=u_{j}u_{i}$), so that the second summation in (3) could be limited to the indices such as $j\ge i$, assuming also $\mathbf {G}_{ij}=0$ for $i\ge j$. In the representation selected throughout the paper, all summations will be full, as in (3) with a fully populated tensor of coefficient $\mathbf {G}$. The same rule applies for the cubic term also. This choice has been made since it allows shorter and simpler expressions for a number of equations given in the presentation, but of course it is not a limiting assumption and the other choice could have also be done.

The first (linear) term in Eq. (2) makes appear the usual tangent stiffness matrix $\mathbf {K}$ defined by :

$$\begin{aligned} \mathbf {K} = \left. \frac{\partial \mathbf {F}}{\partial \mathbf {u}}\right| _{\mathbf {u}=0}, \end{aligned}$$

(4)

from which one can define the eigenmodes, solution of the eigenvalue problem:

$$\begin{aligned} (\mathbf {K} - \omega _i^2 \mathbf {M})\mathbf {\phi }_i = 0, \end{aligned}$$

(5)

with $\mathbf {\phi }_i$ the ith being eigenvector and $\omega _i$ its associated eigenfrequency, for $i=1,\ldots , N$. Using $\mathbf {u} = \mathbf {\Phi }\mathbf {X}$, with $\mathbf {\Phi }$ being the matrix of all eigenvectors $\mathbf {\phi }_i$ and $\mathbf {X}$ the modal coordinates, the problem can be rewritten in the modal basis by premultiplying Eq. (1) by $\mathbf {\Phi }^T$, arriving at:

$$\begin{aligned} {\ddot{\mathbf {X}}} + \mathbf {\Omega }^2 \mathbf {X} + \mathbf {g}\mathbf {X}\mathbf {X}+\mathbf {h}\mathbf {X}\mathbf {X}\mathbf {X}=\mathbf {q}, \end{aligned}$$

(6)

where we have introduced $\mathbf {\Omega }$ the matrix of eigenfrequencies $\omega _i$, $\mathbf {g}$ and $\mathbf {h}$ the tensors of quadratic and cubic coefficients in the modal basis and $\mathbf {q} = \mathbf {\Phi }^T \mathbf {Q}$ the modal external force. The equation of motion in modal space can be written in explicit form with these coefficients as:

$$\begin{aligned} \forall \, p= & {} 1,\ldots , N:\quad {\ddot{X}}_p + \omega _p^2 X_p + \sum _{i=1}^{N} \sum _{j =1}^{N} g_{ij}^p X_i X_j \nonumber \\&+ \sum _{i=1}^{N} \sum _{j=1}^{N} \sum _{k =1}^{N} h_{ijk}^p X_iX_jX_k = q_p. \end{aligned}$$

(7)

The relations between the nonlinear tensors in physical coordinates $\mathbf {G}$ and $\mathbf {H}$ and those in modal coordinates $\mathbf {g}$ and $\mathbf {h}$ are derived from the linear change of coordinates and involve products with the matrix of eigenvectors $\mathbf {\Phi }$. They are provided in “Appendix C” for the sake of completeness.

2.2 NNMs and normal form

Nonlinear normal modes or NNMs have been used since the pioneering work by Rosenberg [43] in numerous vibratory problems. It offers a sound conceptual framework in order to understand the organisation of the dynamics in the phase space. Different definitions have been given in the past, e.g. family of periodic orbits [24, 43], invariant manifold in phase space and tangent at the linear eigenspaces near the origin [46]. More recently, a mathematically well-justified definition of NNM has been provided [16], allowing to settle down the different treatments in a unified way. For that purpose, Haller and Ponsioen proposed to refer to the smoothest member of an invariant manifold family tangent to a modal subbundle along an NNM as a spectral submanifold (SSM). In that sense, SSMs provide a rigorous framework allowing to define the corresponding concepts in all the situations encountered in mechanical vibrations: conservative or dissipative systems and autonomous or non-autonomous systems. Interestingly, the authors also provide in [40] automated formulations in order to derive SSMs up to large order, allowing them to draw out comparisons with numerous other methods proposed in the recent years, see, e.g. [9]. Enforcing the invariant property is key in a perspective of reduced-order modelling, since it is the only way to ensure that the trajectories of the ROM will also exist for the full system. Elaborating on this idea, NNMs have been used in the perspective of model-order reduction using either centre manifold theorem [39, 46], normal form theory [52, 54, 57], or spectral submanifolds [9, 16, 40, 60].

In this contribution, the normal form theory, as defined in [53, 57], is used. The main idea is to define a nonlinear change of coordinates, from the modal coordinates to new ones defined as the normal coordinates. The nonlinear mapping is inherited from Poincaré and Poincaré–Dulac theorems, based on the idea of finding out a nonlinear relationship capable of eliminating as much as possible of nonlinear terms. In this contribution, only the main results are recalled; the interested reader is referred to [52, 53, 57] for more details. The nonlinear change of coordinates is identity-tangent and formally reads:

$$\begin{aligned} X_p&= R_p + {{\mathcal {P}}}_p (R_i,S_j), \end{aligned}$$

(8a)

$$\begin{aligned} Y_p&= S_p + {{\mathcal {Q}}}_p (R_i,S_j), \end{aligned}$$

(8b)

where ${{\mathcal {P}}}_p$ and ${{\mathcal {Q}}}_p$ are the third-order polynomials, the analytical expressions of which are given in [57] for the undamped case and in [53] for the damped case. $X_p$ is the modal coordinate, $Y_p$ the modal velocity, and $(R_i,S_j)$ are the new coordinates related to the invariant manifolds and called normal coordinates.

The method used to derive the nonlinear mapping is based on the recognition of nonlinear resonances involving the eigenfrequencies of the system. In case where no internal resonance is present, one can show, for example, that all the quadratic terms can be cancelled from the normal form which is thus much simpler than the original system.

The dynamics, expressed with the newly introduced normal variables $(R_i,S_j)$, is written in an invariant-based span of the phase space and thus prone to open the doors to efficient reduced-order models, as already shown in [52]. The general equation for the dynamics expressed in the new coordinates reads:

$$\begin{aligned} \forall \; p&= 1,\ldots , N: {\ddot{R}}_p+\omega _p^2 R_p \nonumber \\&\quad +\left( A_{ppp}^p + h_{ppp}^p\right) R_p^3+ B_{ppp}^p R_p S_p^2 \nonumber \\&\quad + R_p \underset{j\ne p}{\sum _{j=1}^N} \left( \left( 3\,h_{pjj}^p+2\,A_{jjp}^p + A_{pjj}^p\right) R_j^2 + B_{pjj}^p S_j^2 \right) \nonumber \\&\quad + S_p \underset{j\ne p}{\sum _{j=1}^N} \left( 2\,B_{jjp}^p R_j S_j\right) =0. \end{aligned}$$

(9)

Note that the expression is slightly different from the one proposed in [57], a direct consequence of the choice of the representation of quadratic and cubic terms, with full summations. The coefficients $A_{ijk}^p$ and $B_{ijk}^p$ stem from the cancellation of the quadratic terms. Their expressions read:

$$\begin{aligned} A_{ijk}^p= & {} \sum _{s =1}^{N} 2\;{\bar{g}}_{is}^p a_{jk}^s, \end{aligned}$$

(10a)

$$\begin{aligned} B_{ijk}^p= & {} \sum _{s =1}^{N} 2\;{\bar{g}}_{is}^p b_{jk}^s, \end{aligned}$$

(10b)

where ${\bar{g}}_{is}^p = \left( g^p_{is} + g^p_{si}\right) /2$ is the mean value between two adjacent terms implying the same monomial term. The coefficients $a_{jk}^s$ and $b_{jk}^s$ appearing in the expression of $A_{ijk}^p$ and $B_{ijk}^p$ are related to the quadratic terms of the change of coordinate. For the sake of completeness, the interested readers can find their full expressions in “Appendix B”. As known from the theory, these second-order coefficients have a singular behaviour in the vicinity of internal resonances. In this case, a strong coupling is present between the nonlinear oscillators whose eigenfrequencies are commensurate, and the associated coefficient in the change of coordinate is set to zero, so that the corresponding monomial terms stay in the normal form.

From Eq. (9), one can observe that invariant-breaking terms are no longer present in the equations of motion. Invariant-breaking terms are defined as quadratic monomials of the form $g^k_{pp}X_p^2$ and cubic monomials $h^k_{ppp}X_p^3$ on kth oscillator equation. As soon as mode p has some energy, these invariant-breaking terms directly excite oscillator k, thus breaking the invariance of the linear mode subspace. As these terms are no longer present in Eq. (9), it shows that the dynamics is now expressed in an invariant-based span. One can also note that the only monomial terms present in Eq. (9) are those related to trivially resonant terms.

A ROM is simply selected by keeping in the truncation only the normal coordinates $(R_p,S_p)$ of interest, depending on the problem at hand. By doing so, one restricts the motion in the invariant manifold described by the master normal coordinates retained, giving rise to efficient reduced models, that simulate trajectories existing in the complete phase space, and allowing to recover the correct type of nonlinearity [55, 57] as well as nonlinear frequency response curves [54]. The simplest ROM is built by restricting the motion to a single NNM by keeping only one pair $(R_p,S_p)$ and cancelling all the other: $\forall \; k\ne p, \quad R_k = S_k = 0$. In this case, the nonlinear change of coordinates for the master coordinates reads:

$$\begin{aligned} X_p&= R_p+a_{pp}^p R_p^2 + b_{pp}^p S_p^2, \end{aligned}$$

(11a)

$$\begin{aligned} Y_p&= S_p + \gamma _{pp}^p R_p S_p, \end{aligned}$$

(11b)

whereas for the slave coordinates, one has:

$$\begin{aligned} \forall \, k&\ne p: X_k = a_{pp}^k R_p^2 + b_{pp}^k S_p^2\nonumber \\&\quad + r_{ppp}^k R_p^3 + u_{ppp}^k R_pS_p^2, \end{aligned}$$

(12a)

$$\begin{aligned} Y_k&= \gamma _{pp}^k R_pS_p+ \mu _{ppp}^k S_p^3 + \nu _{ppp}^k S_pR_p^2. \end{aligned}$$

(12b)

Again, all the introduced coefficients, $\gamma _{pp}^p$, $r_{ppp}^k$, $u_{ppp}^k$, $\mu _{ppp}^k$ and $\nu _{ppp}^k$, originate from the explicit expression of the polynomials ${{\mathcal {P}}}_p$ and ${{\mathcal {Q}}}_p$ of Eq. (8). They are all analytic, and their expressions are given in [57]. Interestingly, Eq. (12) describes the geometry of the invariant manifold in phase space, up to order three, but of course one can limit the development of this equation to second order only. The dynamics on the invariant manifold ($p^{\text{ th }}$ NNM) is found by cancelling all $(R_k,S_k)$ for $k\ne p$ in Eq. (9). In the case of a single-NNM motion, the equation is particularly simple and reads:

$$\begin{aligned} {\ddot{R}}_p + \omega _p^2 R_p + \left( A_{ppp}^p + h_{ppp}^p\right) R_p^3 + B_{ppp}^p R_p {\dot{R}}_p^2 = 0.\nonumber \\ \end{aligned}$$

(13)

Of particular interest here is the fact that the correcting coefficients $A_{ppp}^p$ and $B_{ppp}^p$ appearing in this equation are provided by the second-order terms in the nonlinear change of coordinates. Consequently, the third-order terms have no influence on this reduced dynamics, which is thus exactly the one given by the second-order truncation of the normal form nonlinear mapping.

All these formulas can be used to reconstruct the mode shape dependence on amplitude, assuming the motion is enslaved to a single NNM, i.e. lying in the invariant manifold associated with mode p. Assuming this single-NNM motion, the physical displacement is reconstructed from

$$\begin{aligned} \mathbf {u} = \sum _{k=1}^N X_k \mathbf {\phi }_k = X_p \mathbf {\phi }_p + \underset{k\ne p}{\sum _{k=1}^N}X_k \mathbf {\phi }_k, \end{aligned}$$

(14)

where $X_p$ is replaced using Eq. (11a) and $X_k$ using Eq. (12a), so that one finally obtains the amplitude-dependent mode shape as:

$$\begin{aligned} \mathbf {u}= & {} \left( R_p+a_{pp}^p R_p^2 + b_{pp}^p S_p^2 \right) \mathbf {\phi }_p \nonumber \\&+ \underset{k\ne p}{\sum _{k=1}^N} \left( a_{pp}^k R_p^2 + b_{pp}^k S_p^2+ r_{ppp}^k R_p^3 + u_{ppp}^k R_pS_p^2\right) \mathbf {\phi }_k.\nonumber \\ \end{aligned}$$

(15)

This formula has already been used in order to represent the amplitude dependence of mode shapes on amplitude, see, e.g. [47, 57], and will be further analysed and compared to the prediction given by the method of quadratic manifold from modal derivatives in Sect. 2.4.2.

Note that as a comparison to quadratic manifold is targeted, a detailed description of the effects of order truncation in the normal form approach is in order. In the present approach of the normal form, the change of coordinates is up to order three, but the reduced-order dynamics can be considered as up to the second order, since the effect of cancelling the cubic terms to the higher orders has not been taken into account due to the third-order truncation of all asymptotic developments. Also, most of the comparisons in the remainder of the paper will be drawn between single-mode reduced-order dynamics. In this simplified context, Eq. (13) clearly shows that the cancellation of the third-order non-resonant monomials has absolutely no effect on this equation which is left unchanged. Consequently, Eq. (13) is the reduced dynamics obtained with a second-order normal form nonlinear mapping. The only difference between second-order and third-order nonlinear mapping is in Eq. (12), which describes how the exact invariant manifold is approximated in phase space, and one can analyse the effect of either second-order or third-order nonlinear mapping in this respect. In the remainder of the paper, a clear attention will be devoted to these two specific truncations in order to draw out a fair comparison with the quadratic manifold approach.

We now turn to the definition of modal derivatives and the associated nonlinear mapping: the so-called quadratic manifold, before comparing the two methods in detail.

2.3 Modal derivatives

Modal derivatives have been first introduced by Idehlson and Cardona to solve structural vibrations problems with a nonlinear stiffness matrix [18]. They have been used in recent years in the context of reduced-order modelling [64], and the last developments propose to use them in order to create a nonlinear mapping with a quadratic manifold [19, 44]. In this section, we derive again the most important definitions, make the distinction between modal derivatives (MDs) and static modal derivatives (SMDs) and introduce the quadratic manifold approach.

2.3.1 Definition of modal derivatives and static modal derivatives

The modal derivatives have been first introduced with the aim of offering a framework taking into account the dependence of mode shapes and eigenfrequencies on amplitude for nonlinear system. This is a common point with nonlinear normal modes, which also recognises this fact as a major outcome that needs to be addressed correctly in the modelling. The introduction of the modal derivatives proposed in this section is mostly heuristic and based on previous works. Let us denote $\tilde{\mathbf {\phi }}_i(\mathbf {u})$ this amplitude-dependent eigenvector. The already introduced eigenvector $\mathbf {\phi }_i$, solution of the Sturm–Liouville problem, Eq. (5), represents the value of $\tilde{\mathbf {\phi }}_i(\mathbf {u})$ when $\mathbf {u}=\mathbf {0}$. The ijth modal derivative (MD) is defined as the derivative of $\tilde{\mathbf {\phi }}_i$ with respect to the jth coordinate used for the reduced basis, denoted here as $R_j$. For the sake of clarity, $X_i$ is the modal coordinates and $R_j$ the reduced coordinates, following the notations introduced for the normal form approach. At first order, one has $X_i=R_i$, but as we consider nonlinear change of coordinates, these relationships will be enriched by higher-order terms. For the quadratic manifold approach, this will be explained in the next subsections, so that for the present definitions, one can assume $R_i=X_i$. In that context, the ijth MD $\mathbf {\Theta }_{ij}$ is the derivative of $\tilde{\mathbf {\phi }}_i$ with respect to a displacement enforced along the direction of the jth eigenvector $\mathbf {\phi }_j$ as introduced in [18, 19, 44, 64] and writes:

$$\begin{aligned} \mathbf {\Theta }_{ij} \doteq \frac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {u})}{\partial R_j} \bigg |_{\mathbf {u}=\mathbf {0}}. \end{aligned}$$

(16)

In order to derive an equation from which the MD can be computed, one has to rewrite the eigenproblem given by Eq. (5) assuming the known dependencies on the amplitude, as:

$$\begin{aligned} \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} - {\tilde{\omega }}_i^2(\mathbf {u}) \mathbf {M}\right) \tilde{\mathbf {\phi }}_i(\mathbf {u}) = \mathbf {0}, \end{aligned}$$

(17)

where the linear stiffness matrix is replaced by the full nonlinear restoring force, and both eigenvalues and eigenvectors are amplitude dependent. Note that, in this contribution, the mass matrix is assumed to be independent of the amplitude, since this is the selected framework for this paper focusing on geometric nonlinearity. However, further development could include a dependence of the mass matrix on the amplitude in order to extend the use of MDs to other cases. The nonlinear eigenproblem of Eq. (17) must be complemented with the nonlinear mass normalisation equation:

$$\begin{aligned} \tilde{\mathbf {\phi }}_i(\mathbf {u})^T \mathbf {M} \tilde{\mathbf {\phi }}_i(\mathbf {u})=1. \end{aligned}$$

(18)

Equations (17)–(18) can then be Taylor-expanded as function of the amplitude, assuming moderate vibrations in the vicinity of the position at rest defined by $\mathbf {u}=\mathbf {0}$. Assuming that the displacement $\mathbf {u}$ depends on the coordinates introduced for the reduced basis, $R_1$ to $R_n$, each term can then be expanded along these new coordinates. The full derivation of this Taylor expansion is given in “Appendix E”.

The Taylor expansion of Eqs. (17) and (18) in the $R_j$ coordinates, up to first order, generates constant terms that coincide with the linear eigenproblem and mass normalisation. The next order terms, linear in $R_j$, allow deriving the following system, where the two unknowns are the MD vector $\mathbf {\Theta }_{ij}$ and the scalar describing the variation of the squared eigenfrequency with respect to amplitude, $\frac{\partial \omega _i^2}{\partial R_j}$:

$$\begin{aligned} \left[ \begin{array}{cc} \mathbf {K}-\omega _i^2 \mathbf {M} &{}\quad -\mathbf {M} \mathbf {\phi }_i \\ -\mathbf {\phi }_i^T \mathbf {M} &{}\quad 0 \end{array}\right] \left\{ \begin{array}{c} \mathbf {\Theta }_{ij} \\ \frac{\partial \omega _i^2}{\partial R_j} \end{array}\right\} = \left\{ \begin{array}{c} -2\mathbf {G} \mathbf {\phi }_j \mathbf {\phi }_i \\ 0 \end{array}\right\} ,\nonumber \\ \end{aligned}$$

(19)

where the quadratic tensor $\mathbf {G}$ of the restoring force introduced in Eq. (2) has been used. The detailed proof for the derivation of this system is given in “Appendix E”.

In most of the studies concerned with application of modal derivatives to model order reduction, the so-called static modal derivatives (SMDs) are used instead. Let us denote as $\mathbf {\Theta }^{(S)}_{ij}$ the SMD of $\mathbf {\Theta }_{ij}$, obtained by neglecting the terms related to the mass matrix in (19), which then simplifies to:

$$\begin{aligned} \mathbf {K}\mathbf {\Theta }^{(S)}_{ij}=-2\mathbf {G} \mathbf {\phi }_j \mathbf {\phi }_i. \end{aligned}$$

(20)

This simplification evidently highlights the fact that MDs and SMDs are able to retrieve the quadratic coupling generated by the nonlinear restoring force, since being directly proportional to the tensor of coefficients $\mathbf {G}$. Equation (20) also shows that the computation of SMDs is drastically reduced as compared to MDs, for two main reasons. The first one is that given the usual symmetry of the quadratic tensor $\mathbf {G}$ at hand in structural problems, one has $\mathbf {G} \mathbf {\phi }_j \mathbf {\phi }_i=\mathbf {G} \mathbf {\phi }_i \mathbf {\phi }_j$, so that the SMDs are symmetric $\mathbf {\Theta }^{(S)}_{ij}=\mathbf {\Theta }^{(S)}_{ji}$. This involves that the number of calculations for indexes $i\ne j$ is then halved in the case of SMDs as compared to MDs. The second reason lies in the fact that although the sizes of the systems to solve are comparable (the size of system (19) is $N+1$ and the size of system (20) is N), the computation of a SMD can be done with a standard operation in commercial FE software, whereas the computation of a MD cannot. Indeed, the non-intrusive computation of a SMD requires to solve a linear system $\mathbf {K}\mathbf {u} = \mathbf {f}$, where the applied force $\mathbf {f}$ is the right-hand side of Eq. (20) and the resultant displacement $\mathbf {u}$ is the SMD. Solving such linear system coincides with operating a simple linear static analysis on the structure with imposed force and unknown displacement. Conversely, the linear system to compute a MD is the one in Eq. (19). The solution of this system does not correspond to the standard operation one could easily perform in FE software. Consequently, to compute the MD, one needs not only to access the full stiffness and mass matrices but also to export them in an external code to be able to solve the linear system. When the structure is discretised with a large number of dofs, such operation can be memory and time-consuming when not infeasible.

2.3.2 Expression of MDs as function of the quadratic coefficients from the modal basis

In this section, the relation between MDs and SMDs and the coefficients of the quadratic tensor in modal basis $\mathbf {g}$ are derived. This relation will help to draw comparisons between the normal form method and the quadratic manifold method that will be introduced in the next section. For that purpose, the ijth MD in the modal basis, denoted as $\mathbf {\theta }_{ij}$, is introduced as

$$\begin{aligned} \quad \mathbf {\Theta }_{ij}=\mathbf {\Phi }\mathbf {\theta }_{ij}=\sum _{s=1}^N \mathbf {\phi }_s \theta ^s_{ij}, \end{aligned}$$

(21)

following the linear change of basis from physical to modal space, where the summation thus spans over all the modes of the structure, being $\mathbf {\Phi }$ the full eigenvector matrix. In the modal basis, the eigenvector $\mathbf {\phi }_i$ coincides with the ith vector basis $\mathbf {e}_i$, where the entries of $\mathbf {e}_i$ are all zero except 1 in position i, so that: $\mathbf {\phi }_i=\mathbf {\Phi }\mathbf {e}_i$.

The system of Eq. (19) can be now written in modal coordinates by premultiplying the first N rows by $\mathbf {\Phi }^T$ and by substituting the values of $\mathbf {\phi }_i$ and $\mathbf {\Theta }_{ij}$ with their values in modal coordinates. One finally obtains:

$$\begin{aligned} \left[ \begin{array}{cc} \mathbf {\Omega }^2-\omega _i^2 \mathbf {I} &{}\quad - \mathbf {e}_i \\ -\mathbf {e}_i^T &{} \quad 0 \end{array}\right] \left\{ \begin{array}{c} \mathbf {\theta }_{ij} \\ \frac{\partial \omega _i^2}{\partial R_j} \end{array}\right\} = \left\{ \begin{array}{c} -2\mathbf {g}_{ij} \\ 0 \end{array}\right\} , \end{aligned}$$

(22)

where the right-hand side has been simplified using the relationship $\mathbf {g}_{ij}=\mathbf {\Phi }^T \mathbf {G} \mathbf {\phi }_i \mathbf {\phi }_j$, demonstrated in Eq. (84a) of “Appendix C”.

System (22) is easier to understand when written by term:

$$\begin{aligned}&(\omega _s^2-\omega _i^2)\theta ^s_{ij}=-2g^s_{ij},\quad \text{ for } \quad s\ne i, \end{aligned}$$

(23a)

$$\begin{aligned}&\frac{\partial \omega _i^2}{\partial R_j}=2g^i_{ij}, \quad \text{ for } \quad s=i, \end{aligned}$$

(23b)

$$\begin{aligned}&\theta ^i_{ij}=0. \end{aligned}$$

(23c)

One can notice that the ijth modal derivatives are then directly proportional to the ijth component of the quadratic tensor in modal coordinates. This clearly shows that the ijth MD is able to retrieve a strong quadratic coupling occurring between slave mode s and the master modes i and j. The value of the modal derivative in physical coordinates can be now easily reconstructed from the preceding development and reads:

$$\begin{aligned} \mathbf {\Theta }_{ij} = \underset{s\ne i}{\sum _{s=1}^N} \mathbf {\phi }_s \;\dfrac{-2\;g^s_{ij}}{\omega _s^2-\omega _i^2}. \end{aligned}$$

(24)

If one follows a similar procedure for the case of static modal derivative, Eq. (20) is written in modal coordinates as $\mathbf {\Omega }^2\mathbf {\theta }^{(S)}_{ij}=-2\mathbf {g}_{ij}$ and the static modal derivative in physical coordinates is directly given as:

$$\begin{aligned} \mathbf {\Theta }^{(S)}_{ij}=\sum _{s=1}^N \mathbf {\phi }_s \;\dfrac{-2\;g^s_{ij}}{\omega _s^2}. \end{aligned}$$

(25)

In both cases, MDs and SMDs can be simply defined as a linear combination of modes weighted by a factor proportional to $g^s_{ij}$, the quadratic modal coupling coefficient. In the case of modal derivative, the method shows a divergent behaviour in case of 1:1 internal resonance between two eigenfrequencies, a feature that will be further commented in Sect. 2.4.1. One can also note that the weighting factors have larger values for the modes, the eigenfrequencies of which are closer to the eigenfrequency of the ith mode. On the other hand, for static modal derivatives, the weighting factors are simply proportional to the inverse of the squared eigenfrequencies and thus should decrease for higher modes. Note, however, that this fact can be severely compensated by the values of the quadratic coefficients, which scales according to the linear stiffness. Consequently, as shown, for example, in [49, 61] for thin and flat symmetric structures (beams and plates), the SMD is able to recover the most important couplings with in-plane modes.

As a conclusion, MDs and SMDs can be seen as a displacement field that takes into account the contribution of all quadratically coupled modes into one equivalent vector. From this perspective, the use of a reduced basis composed of MDs is equivalent to using a basis composed of all quadratically coupled modes, with the supplementary condition that the quadratic couplings makes appear new directions in phase space that are independent of the already selected mode. If the quadratic coupling is only dependent on modes already present in the reduced basis, then the new vector will not bring out new eigendirections.

2.3.3 Quadratic manifold

The quadratic manifold approach has been introduced in [19, 44] in order to extend the use of modal derivatives in the context of model order reduction and propose a nonlinear mapping from initial to reduced coordinates. The nonlinear mapping is quadratic in nature and does not account for nonlinear internal resonance as the normal form theory does. In this section, the derivation of reduced-order models using the quadratic manifold is given, following the previous results obtained in [19, 44]. A particular attention is paid on writing the differences one can await when using the quadratic manifold with MDs and SMDs, in comparison with the results provided by normal form theory in mind, thus giving rise to new developments. The coordinates describing the reduced-order models are denoted as $R_p$ for all the methods in order to compare more directly the equations. One has, however, to keep in mind that the meaning of these coordinates is not the same for each method.

Since the MDs are defined from a second-order Taylor expansion of the nonlinear eigenvalue problem, it is intuitive to use them in a quadratic nonlinear mapping. If one operates a Taylor expansion of the approximate solution $\mathbf {u}$ in the reduced coordinates $\mathbf {R}$ up to quadratic order, one finds:

$$\begin{aligned} \mathbf {u}(\mathbf {R})= & {} \mathbf {u}(\mathbf {0}) + \sum ^n_{i=1}\dfrac{\partial \mathbf {u}(\mathbf {R})}{\partial R_i}\bigg |_{\mathbf {0}}R_i \nonumber \\&+ \frac{1}{2} \sum ^n_{i=1}\sum ^n_{j=1} \left( \dfrac{\partial ^2 \mathbf {u}(\mathbf {R})}{\partial R_j \partial R_i} \right) \bigg |_{\mathbf {0}}R_iR_j + {\mathcal {O}}(|\mathbf {R}|^3),\nonumber \\ \end{aligned}$$

(26)

where n is the number of master modes retained for the ROMS, $\mathbf {R} = (R_1, \ldots , R_n)$. By extending the definition of linear eigenvectors to the nonlinear ones, the nonlinear eigenvector spans the tangent space of the displacement with respect to the reduced coordinates, so that:

$$\begin{aligned} \dfrac{\partial \mathbf {u}}{\partial R_i} = \tilde{\mathbf {\phi }}_i (\mathbf {R}). \end{aligned}$$

(27)

In Eq. (26), we can then substitute $\mathbf {u}(\mathbf {0}) = \mathbf {0}$ (the position at rest is at the origin of the coordinates), and

$$\begin{aligned} \dfrac{\partial \mathbf {u}(\mathbf {R})}{\partial R_i}\bigg |_{\mathbf {0}}&= \mathbf {\phi }_i, \end{aligned}$$

(28)

$$\begin{aligned} \dfrac{\partial ^2 \mathbf {u}(\mathbf {R})}{\partial R_i \partial R_j}\bigg |_{\mathbf {0}}&= \mathbf {\Theta }_{ij}. \end{aligned}$$

(29)

However, this series of operations would lead to an inconsistent formulation in the case of MDs due to their asymmetry, as already outlined in [19]. In fact, since $\mathbf {\Theta }_{ij} \ne \mathbf {\Theta }_{ji}$, it implies that the Schwarz’s identity $\partial ^2 \mathbf {u}/ \partial R_i \partial R_j \ne \partial ^2 \mathbf {u}/ \partial R_j \partial R_i$ is not fulfilled anymore. To overcome this issue, and given the independence of the quadratic mapping on the asymmetric part of each MD shown in [19], the correct strategy proposed in [19] is to express both the mapping and its tangent space by means of symmetrised MDs $\bar{\mathbf {\Theta }}_{ij}=(\mathbf {\Theta }_{ij}+\mathbf {\Theta }_{ji})/2$, leading to:

$$\begin{aligned} \mathbf {u}(\mathbf {R})&\approx \sum ^n_{i=1}\mathbf {\phi }_iR_i + \frac{1}{2}\sum ^n_{i=1}\sum ^n_{j=1} \bar{\mathbf {\Theta }}_{ij} R_iR_j = \mathbf {\phi }\mathbf {R}\nonumber \\&\quad +\frac{1}{2}\bar{\mathbf {\Theta }}\mathbf {R}\mathbf {R}, \end{aligned}$$

(30)

$$\begin{aligned} \tilde{\mathbf {\phi }}_i(\mathbf {R})&\approx \mathbf {\phi }_i + \sum ^n_{j=1}\bar{\mathbf {\Theta }}_{ij} R_j= \mathbf {\phi }_i+\bar{\mathbf {\Theta }}\mathbf {R}. \end{aligned}$$

(31)

Note that these expressions are used in order to define the reduced-order model, so the dimension n of $\mathbf {R}$ is much smaller than the dimension N of $\mathbf {u}$, $n\ll N$, since only the master coordinates of the ROM are present in $\mathbf {R}$. Consequently, $\mathbf {\phi }$ is the matrix of eigenvectors relative to the master coordinates and should be distinguished from the full matrix of eigenvectors $\mathbf {\Phi }$ used, for example, in (21). Finally, $\bar{\mathbf {\Theta }}$ is the third-order tensor gathering the MDs $\bar{\mathbf {\Theta }}_{ij}$.

For future comparison with the normal form method, it is useful to also define the quadratic mapping in modal coordinates:

$$\begin{aligned} \mathbf {X}(\mathbf {R}) \approx \mathbf {R}+\frac{1}{2}\bar{\mathbf {\theta }}\mathbf {R}\mathbf {R}, \end{aligned}$$

(32)

and by components:

$$\begin{aligned} X_k \approx R_k+\frac{1}{2}\sum ^n_{i=1}\sum ^n_{j=1} {\bar{\theta }}^k_{ij}R_iR_j. \end{aligned}$$

(33)

2.3.4 Reduced-order model obtained with quadratic manifold

The nonlinear mapping can then be used in order to derive the reduced-order equations by directly applying Eq. (30) to the original equations of motion, Eq. (1), and using a standard Galerkin projection. For that purpose, one has to compute the derivatives of Eq. (30) with respect to time and finally left-multiply Eq. (1) by $\tilde{\mathbf {\phi }}_i^T$. These derivations have already been proposed in [19, 44], and we refer the interested reader to these articles for details about the procedure. Here, we give the reduced-order dynamics obtained once the projection realised, as a function of the modal coupling coefficients $\mathbf {g}$ and $\mathbf {h}$, a derivation that is not given in [19, 44] and will allow drawing out more direct comparisons with the normal form approach.

The dynamics for each reduced coordinates $R_p$ finally reads, for $p=1,\ldots ,N$:

$$\begin{aligned} \begin{aligned}&\ddot{R}_p + \omega ^2_{p} R_p + \sum _{i,j=1}^n \left( \left( g^p_{ij} +\frac{\omega ^2_{p}}{2}\,{\bar{\theta }}^p_{ij}\right) \;R_i R_j \right. \\&\quad \left. +\,{\bar{\theta }}^p_{ij}\;({\dot{R}}_i {\dot{R}}_j +R_i\ddot{R}_j) +\; {\bar{\theta }}^j_{pi}\;(\omega ^2_{j}R_i R_j +R_i \ddot{R}_j ) \right) \\&\quad +\,\sum _{i,j,k=1}^n \left( \left( h^p_{ijk}+ \sum _{s=1}^n \left( {\bar{g}}^p_{is}\;{\bar{\theta }}^s_{jk}\right. \right. \right. \\&\quad \left. \left. \left. +\, {\bar{\theta }}^s_{pk}\; \left( g^s_{ij}+\frac{\omega ^2_{s}}{2}\,{\bar{\theta }}^s_{ij}\right) \right) \right) \; R_i R_j R_k + \sum _{s=1}^n \left( {\bar{\theta }}^s_{pk}\;{\bar{\theta }}^s_{ij} \right) \right. \\&\quad \left. ({\dot{R}}_i {\dot{R}}_j R_k+\ddot{R}_i R_jR_k ) \right) = 0, \end{aligned} \end{aligned}$$

(34)

where the following notations have been introduced for simplifying the expressions : ${\bar{g}}^p_{is}=\frac{{g}^p_{is}+{g}^p_{si}}{2}$. Note that this formula simplifies in the case of a symmetric quadratic tensor, which is generally the case in structural mechanics.

One can observe that the linear part is uncoupled, resulting from the fact that the first term of the quadratic manifold in Eq. (30) is the usual expansion onto the eigenmodes, thus implying, at linear order, uncoupled linear oscillators. The nonlinear terms can be compared to those obtained when using the normal form approach as nonlinear mapping, Eq. (9). In particular, one can observe that the normal form approach completely cancels all quadratic terms, provided that no second-order internal resonance is present, a key feature embedded in the derivation which makes the distinction between resonant and non-resonant terms. On the other hand, quadratic terms are always present in (34). A second comment is on the presence of terms depending on accelerations in (34), not present in the reduced-order dynamics given by the normal form approach.

The restriction to a single master dof is provided, so that one could draw out a term-by-term comparison between the reduced-order dynamics provided by the two methods. Assuming that only mode p is present as reduced coordinates and thus $R_i = 0$, for all $i\ne p$, Eq. (34) simplifies to:

$$\begin{aligned} \begin{aligned}&\ddot{R}_p + \omega ^2_{p} R_p + \left( g^p_{pp} +\frac{\omega ^2_{p}}{2}\,\theta ^p_{pp}\right) \;\\&R_p^2 + \theta ^p_{pp}\;({\dot{R}}_p^2 +R_p\ddot{R}_p) + \theta ^p_{pp}\;(\omega ^2_{p}R_p^2 +R_p \ddot{R}_p ) \\&\quad +\; h^p_{ppp}R_p^3 + \sum _{s=1}^n \left( {\bar{g}}^p_{ps}\;\theta ^s_{pp}\;R_p^3 \right. \\&\quad \left. + \theta ^s_{pp}\; \left( g^s_{pp}+\frac{\omega ^2_{s}}{2}\,\theta ^s_{pp}\right) \; \right. \\&\left. R_p^3 + \left( \theta ^s_{pp}\right) ^2\; ({\dot{R}}_p^2 R_p+\ddot{R}_p R_p^2 ) \right) = 0. \end{aligned} \end{aligned}$$

(35)

This equation can then be used for either MD or SMD, so that one can contrast the results obtained by using one of these two strategies (modal derivatives, be they static or dynamic) with the nonlinear change of coordinates provided by normal form theory, which is the aim of the next section.

2.4 Comparison of the methods and slow/fast approximation

This section aims at comparing the different nonlinear mappings used to derive reduced-order models on the different outcomes they provide: reduced-order dynamics, and prediction of typical nonlinear features such as hardening/softening behaviour, and dependence of mode shapes on amplitude. For that purpose, we restrict ourselves to a single master mode. Moreover, from now on, we introduce the symmetry property of the quadratic tensor $\mathbf {g}$, which results from the fact that the internal force derives from a potential, thus leading to $g^i_{jk}=g^i_{kj}$ and $g^i_{jk}=g^j_{ki}=g^k_{ij}$. Note, however, that, due to our initial choice of fully populated sums and tensors without assuming commutativity of the product, the symmetry property may appear a bit different from, for example, [34] when equal indexes are present. Indeed, in [34] one can read, for example, $g^p_{ps}=2g^s_{pp}$. This is the only consequence of the initial choice since in [34] one has $g^p_{sp}=0$ for $s>p$. In our case, the relationship reads $g^p_{ps}=g^s_{pp}$ and $g^p_{sp}=g^s_{pp}$.

By using such symmetry property, we can also simplify ${\bar{g}}^p_{ps}=g^p_{ps}=g^s_{pp}$, and substituting the value of the modal derivative in modal space $\theta ^s_{pp}=-2g^s_{pp}/(\omega ^2_{s}-\omega ^2_{p})$ when $s\ne p$ and $\theta ^p_{pp}=0$ in Eq. (35), one obtains:

$$\begin{aligned} \begin{aligned}&\ddot{R}_p + \omega ^2_{p} R_p + g^p_{pp}\;R_p^2 + h^p_{ppp}R_p^3 \\&\quad - \underset{s\ne p}{\sum _{s=1}^n } (g^s_{pp})^2 \frac{2}{\omega ^2_{s}-\omega ^2_{p}} \left( \frac{\omega ^2_{s}-2\omega ^2_{p}}{\omega ^2_{s}-\omega ^2_{p}}\;R_p^3\right. \\&\quad \left. - \frac{2}{\omega ^2_{s}-\omega ^2_{p}}\; ({\dot{R}}_p^2 R_p+\ddot{R}_p R_p^2 ) \right) = 0. \end{aligned} \end{aligned}$$

(36)

If the value of the SMD is used instead of the MD, then the reduced-order dynamics writes:

$$\begin{aligned} \begin{aligned}&\ddot{R}_p + \omega ^2_{p} R_p - g^p_{pp}\frac{2}{\omega ^2_{p}}\;\left( \omega ^2_{p}R_p^2+{\dot{R}}_p^2 +2R_p\ddot{R}_p\right) \\&\quad + h^p_{ppp}R_p^3 - \sum _{s=1}^n \left( g^s_{pp}\right) ^2\frac{2}{\omega ^2_{s}} \left( R_p^3 \right. \\&\quad \left. - \frac{2}{\omega ^2_{s}}\left( {\dot{R}}_p^2 R_p+\ddot{R}_p R_p^2 \right) \right) = 0. \end{aligned} \end{aligned}$$

(37)

For the explicit comparison, we rewrite the reduced-order dynamics derived with the normal form approach, Eq. (13), where the $A^p_{ppp}$ and $B^p_{ppp}$ terms have been expanded:

$$\begin{aligned} \begin{aligned}&\ddot{R}_p + \omega _p^2 R_p + h_{ppp}^p R_p^3 - \sum _{s=1}^n \left( {g}^s_{pp}\right) ^2 \frac{2}{\omega ^2_{s}} \\&\quad \left( \frac{\omega ^2_s-2 \omega ^2_p}{\omega ^2_s-4 \omega ^2_p}\;R_p^3 - \frac{2}{\omega ^2_s-4 \omega ^2_p}\;{\dot{R}}_p^2 R_p \right) = 0. \end{aligned} \end{aligned}$$

(38)

Note that the remark on the order of the truncations given at the end of Sect. 2.2 may be better understood from these single-mode reduced dynamics. Equation (38) is the ROM given by normal form, be the calculation of the nonlinear change of coordinate truncated at order two or at order three. Consequently, this equation gives the third-order reduced dynamics produced by truncating the normal form at second order. In the same line, Eqs. (36) and (37) are the third-order reduced-dynamics provided by the quadratic manifold approach. Hence, comparing the predictions given by these reduced dynamics is correct since the same order of asymptotic developments is at hand. The only difference one can estimate in the analysis thus relies in the nonlinear mapping, which can be pushed at third order easily in the normal form approach since the calculation has already been proposed in the past. This means that in the comparisons, the only difference will be on the geometry of the manifold in phase space and the reconstruction formula, but not on the reduced-order dynamics.

In order to have a better view on the reduced-order dynamics for each of the methods, the general nonlinear oscillator equation describing the dynamics on the reduced subspace can be written under the general form as:

$$\begin{aligned}&\ddot{R}_p + \omega _p^2 R_p + C_1 R_p^2 + C_2 \dfrac{{\dot{R}}_p^2}{\omega ^2_{p}}+ C_3 \dfrac{\ddot{R}_p R_p}{\omega ^2_{p}}\nonumber \\&\quad +\,C_4 R_p^3 + C_5 \dfrac{{\dot{R}}_p^2 R_p}{\omega ^2_{p}}+ C_6 \dfrac{\ddot{R}_p R_p^2}{\omega ^2_{p}} = 0, \end{aligned}$$

(39)

with $C_1$ to $C_6$ being different coefficients, whose values are summarised in Tables 1 and 2 for the three different methods.

Table 1 Table of coefficients of the reduced system given by the three methods

Table 2 Table of coefficients of the reduced system given by the three methods

As already remarked, only the normal form approach is able to cancel the quadratic nonlinearity and produce a parsimonious, cubic-order reduced dynamics, depending on two separate coefficients only. Using SMDs creates the larger number of coefficients, while only four are needed for MDs. Most importantly, the closeness of the results given by the three methods can be underlined in the case where a slow/fast decomposition can be assumed between the master mode p and the slave modes s. This case is often encountered in mechanical vibrations since one has often to deal with a large number of modes with very high eigenfrequencies. Let us assume that all the slave modes s are well separated from the master mode, so that for all s one has $\omega _s \gg \omega _p$. It is then very easy to verify on the coefficients given in Tables 1 and 2 that those provided by the normal form and the MD method tend to the values given by the SMD approach. More specifically, $C_4$ and $C_5$ from normal form exactly match those from the SMD, so that the only difference between the two reduced-order dynamics lies in the additional terms $C_1$, $C_2$, $C_3$ and $C_6$ for the SMD method. On the other hand, using the slow/fast approximation for the coefficients provided by the MD shows that $C_4$, $C_5$ and $C_6$ tend exactly to the values obtained with SMDs, the only difference being in the summation, where the p term is excluded in the MD approach, whereas it is not in the SMD, as a direct consequence from Eq. (23). Indeed, Eq. (23b) shows that for MD, the $g^p_{pp}$ term is taken into account in the amplitude–frequency relationship, and not in the reconstruction of the vector as given by Eq. (24). On the other hand, for SMD, the $g^p_{pp}$ term is taken into account in the vector defining the SMD, Eq. (25), but not in the frequency dependence on amplitude. This important difference between the two methods will have consequences that are commented further in the next sections, and the $g^p_{pp}$ will be denoted further as the self-quadratic term.

In order to better understand the observed differences on the reduced-order dynamics, a fair comparison has to be given not onto a term-by-term comparison, since the meaning of the reduced variables is not the same, but on the general predictions given by each reduction method on the most important nonlinear features. The next sections are thus devoted to comparing the prediction of the type of nonlinearity provided by each method (i.e. the first term in the amplitude–frequency relationship that dictates the hardening or softening behaviour), as well as the mode shape dependence on amplitude.

2.4.1 Hardening/softening behaviour

The generic reduced-order dynamics, Eq. (39), can be solved with a perturbation method in order to derive the type of nonlinearity predicted by each method. Keeping the general notation with the $C_i$ coefficients for the ease of reading, the general solution up to second order in amplitude reads:

$$\begin{aligned} R_p= & {} a_0 \cos \left[ \,\omega _p\,t\,( 1 + \varGamma a_0^2)\right] \nonumber \\&+\,a_0^2 \left( \frac{C_1 - C_2 - C_3}{6\,\omega _p^2}\cos \left[ \,2\,\omega _p\,t\,\left( 1 + a_0^2 \varGamma \right) \right] \right. \nonumber \\&\left. -\,\frac{C_1 + C_2 - C_3}{2\,\omega _p^2}\right) + {\mathcal {O}}(a_0^3), \end{aligned}$$

(40)

with $a_0$ being the amplitude and $\varGamma$ the general coefficient that dictates the hardening/softening behaviour. Indeed, one can introduce the nonlinear frequency $\omega _{\mathrm {NL}} = \omega _p( 1 + \varGamma a_0^2)$. If $\varGamma > 0$, then the system is hardening. The general expression for $\varGamma$ with all the $C_i$ coefficients writes:

$$\begin{aligned} \varGamma= & {} -\frac{1}{24\,\omega _p^4}\left( 10\,C_1^2 + 10\,C_1 C_2 + 4\,C_2^2 - 7\,C_2 C_3 \right. \nonumber \\&\left. + \,C_3^2 - 11\,C_1 C_3 \right) +\frac{1}{8\,\omega _p^2}\left( 3\,C_4 +\,C_5 -3\,C_6\right) .\nonumber \\ \end{aligned}$$

(41)

One can note in particular that with the normal form approach, one has $C_1=C_2=C_3=0$ since the method fully cancels the quadratic terms, so that there is no second harmonic term in the reduced-order dynamics and Eq. (40) reduces to its first term at order two. However, since quadratic terms are present in the nonlinear change of coordinates, this simplification does not imply that the second harmonic is not present in the reconstructed displacements, as it will be shown in the next section. Once again, these two last equations show that normal form approach produces a parsimonious representation of the reduced dynamics which is generally easier to read and interpret.

Replacing the values of the $C_i$ coefficients obtained for each method (MD, SMD or NF for normal form), one arrives at the prediction of the type of nonlinearity provided by each reduced-order model as:

$$\begin{aligned}&\varGamma _{\text {MD}}= - \dfrac{5}{12\;\omega ^2_{p}} \left( \dfrac{g^p_{pp}}{\omega _{p}}\right) ^2 + \dfrac{3}{8\;\omega ^2_{p}}\nonumber \\&\quad \left( h_{ppp}^p - \sum _{\underset{s\ne p}{s=1}}^n 2\left( \dfrac{g^s_{pp}}{\omega _{s}}\right) ^2 \left( 1+ \dfrac{\omega ^2_p(4\,\omega ^2_s-3\,\omega ^2_p)}{3(\omega ^2_s-\omega ^2_p)^2} \right) \right) , \end{aligned}$$

(42a)

$$\begin{aligned}&\varGamma _{\text {SMD}}= - \dfrac{5}{12\;\omega ^2_{p}} \left( \dfrac{g^p_{pp}}{\omega _{p}}\right) ^2 + \dfrac{3}{8\;\omega ^2_{p}}\nonumber \\&\quad \left( h_{ppp}^p - \sum _{\underset{s\ne p}{s=1}}^n 2\left( \dfrac{g^s_{pp}}{\omega _{s}}\right) ^2 \left( 1+\dfrac{4\,\omega ^2_p}{3\,\omega ^2_s} \right) \right) , \end{aligned}$$

(42b)

$$\begin{aligned}&\varGamma _{\text {NF}}= - \dfrac{5}{12\;\omega ^2_{p}} \left( \dfrac{g^p_{pp}}{\omega _{p}}\right) ^2 + \dfrac{3}{8\;\omega ^2_{p}}\nonumber \\&\quad \left( h_{ppp}^p - \sum _{\underset{s\ne p}{s=1}}^n 2\left( \dfrac{g^s_{pp}}{\omega _{s}}\right) ^2 \left( 1+ \dfrac{4\,\omega ^2_p}{3(\omega ^2_s-4\,\omega ^2_p)} \right) \right) . \end{aligned}$$

(42c)

One can note that the first terms of the prediction are the same, while the difference arises from the way the slave (or neglected) coordinates are taken into account in order to predict the type of nonlinearity. This feature is, however, key in order to give a correct prediction since there is a strong need to take properly into account the curvature of the manifolds in phase space; otherwise, incorrect predictions are given [57].

In order to give more insights into Eq. (42), let us first notice that in the summed terms, the first one is always the same since the different expressions all start with $1+\cdots$. Let us isolate this term and introduce the following notation :

$$\begin{aligned} {\mathcal {C}}^s_{\text {SC}} = 2\left( \dfrac{g^s_{pp}}{\omega _{s}}\right) ^2. \end{aligned}$$

(43)

One can notice that this correction term is the one obtained by using static condensation, as already shown, for example, in [48, 61], and thus the subscript $\text {SC}$. Denoting as ${\mathcal {C}}_{\text {MD}}$, ${\mathcal {C}}_{\text {SMD}}$ and ${\mathcal {C}}_{\text {NF}}$ the correction factors given by each method (i.e. the term in the summation), one can then simply compares all these terms to ${\mathcal {C}}^s_{\text {SC}}$ in order to have an expression depending only on the eigenfrequencies. Assuming that there is only one slave mode s in the summation in order to highlight the contribution brought by each term, the following ratios can be written:

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {MD}}}{{\mathcal {C}}_{\text {SC}}}\;\,=\;&1+\dfrac{4}{3} \dfrac{\rho ^2-3/4}{(\rho ^2-1)^2}, \end{aligned}$$

(44a)

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {SMD}}}{{\mathcal {C}}_{\text {SC}}}=\;&1+\dfrac{4}{3} \dfrac{1}{\rho ^2}, \end{aligned}$$

(44b)

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {NF}}}{{\mathcal {C}}_{\text {SC}}}\;\,=\;&1+\dfrac{4}{3} \dfrac{1}{\rho ^2-4}, \end{aligned}$$

(44c)

where $\rho =\omega _s/\omega _p$ has been introduced in order to highlight their behaviour with respect to the fulfilment of the slow/fast partition. These expressions clearly underline the fact that each method refines the correction factor of static condensation by an additional term. One can also observe that the refinement of the ${\mathcal {C}}^s_{\text {SMD}}$ comes from the inertia and velocity terms $C_5$ and $C_6$, whereas the term $C_4$ is exactly the one from static condensation. Consequently, using SMD without quadratic manifold could lead to erroneous predictions since inertial and velocity corrections could be missed. This remark should be particularly relevant in a case of geometric nonlinearity involving inertia, as, for example, in the case of a cantilever beam.

To better assess the quality of the predictions given by the three methods, Eq. (44) can be Taylor-expanded by using the slow/fast assumption $\omega _s \gg \omega _p$ for the slave modes s. This assumption allows introducing a small parameter $\omega _p/\omega _s$, or, equivalently, considering the expansion under the assumption $\rho \rightarrow \infty$. One then obtains:

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {MD}}}{{\mathcal {C}}_{\text {SC}}}\;\,=\;&1+\dfrac{4}{3}\dfrac{1}{\rho ^2}+ \sum _{i=2}^\infty \dfrac{3+i}{3\;\rho ^{2i}}, \end{aligned}$$

(45a)

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {SMD}}}{{\mathcal {C}}_{\text {SC}}}=\;&1+\dfrac{4}{3} \dfrac{1}{\rho ^2}, \end{aligned}$$

(45b)

$$\begin{aligned} \dfrac{{\mathcal {C}}_{\text {NF}}}{{\mathcal {C}}_{\text {SC}}}\;\,=\;&1+\dfrac{4}{3}\dfrac{1}{\rho ^2}+ \sum _{i=2}^\infty \dfrac{4^i}{3\;\rho ^{2i}}. \end{aligned}$$

(45c)

These formulas show in particular that all the methods predict the same first two terms in the expansion that assumes slow/fast partition, and the limit for $\rho \rightarrow \infty$ is the same for all methods, including static condensation, since the ratios tends to 1 in this case. This means that a formal equivalence in the prediction of the type of nonlinearity is obtained only in the limit case of $\omega _s \gg \omega _p$ for all the studied methods. Figure 1 illustrates this convergence and shows that it is obtained rapidly, indicating, in particular, that from the value $\omega _s / \omega _p \simeq 4$, all methods are almost converged in terms of type of nonlinearity, thus quantifying more properly the value from which the slow/fast partition is effective so that one can use the methods based on modal derivatives safely. In order to be a bit more quantitative, one can remark that the relative difference between ${\mathcal {C}}_{\text {MD}}$ and ${\mathcal {C}}_{\text {NF}}$ is equal to 5% for $\rho =3.25$ and 1% for $\rho =4.6$, so that the proposed bound $\omega _s / \omega _p \simeq 4$ has not to be understood as a strict one. Moreover, the error on $\varGamma$ will be smaller than the error on the correction factor ${\mathcal {C}}$, being $\varGamma$ composed of other terms that are not affected by the reduction method. The conclusion is that $\rho \in [3,4]$ can be understood as a transition region, and converged results thanks to slow/fast assumption can be faithfully obtained over 4, but below 3 caution has to be exercised.

Fig. 1

Evolution of the ratios $\dfrac{{\mathcal {C}}_{\text {MD}}}{{\mathcal {C}}_{\text {SC}}}$, $\dfrac{{\mathcal {C}}_{\text {SMD}}}{{\mathcal {C}}_{\text {SC}}}$ and $\dfrac{{\mathcal {C}}_{\text {NF}}}{{\mathcal {C}}_{\text {SC}}}$, defined in Eq. (44), as a function of the parameter $\rho =\omega _s/\omega _p$, from which the behaviour of the type of nonlinearity defined by the $\varGamma$ coefficients in Eq. (42) can be directly deduced. Dashed grey line is the (constant) value predicted by static condensation. (All curves are normalised with respect to this value.) Yellow curve: $\dfrac{{\mathcal {C}}_{\text {SMD}}}{{\mathcal {C}}_{\text {SC}}}$ predicted by static modal derivatives; orange curve: $\dfrac{{\mathcal {C}}_{\text {MD}}}{{\mathcal {C}}_{\text {SC}}}$ computed from modal derivatives; blue curve: $\dfrac{{\mathcal {C}}_{\text {NF}}}{{\mathcal {C}}_{\text {SC}}}$ given by normal form theory. (Color figure online)

Figure 1 shows also other interesting features on the behaviour of the type of nonlinearity. Besides the convergence of all curves in the limit $\rho \rightarrow \infty$, important differences occur in the regions where the methods have a singularity. The normal form approach displays a singular behaviour in the vicinity of the 1:2 internal resonance when $\omega _s \simeq 2\omega _p$. This fact is logical and has already been commented in numerous prior publications. Indeed, when such a resonance exists, then a strong coupling arises between the two modes, so that reducing the dynamics to a single master mode has no meaning anymore, and the minimal model should be composed at least by these two internally resonant modes. The divergence in the behaviour of ${\mathcal {C}}_{\text {NF}} / {\mathcal {C}}_{\text {SC}}$ reflects this fact, meaning that in this zone the definition of the type of nonlinearity is of no more use since another dynamical regime takes place. Previous publications also clearly underline that the prediction given by $\varGamma _{\text {NF}}$ in Eq. (42c) is correct [57], which has been confirmed with comparisons to direct simulations of the full-order model, and this prediction of the type of nonlinearity has then been used for continuous structures such as cables and shells [6, 38, 41, 55].

On the other hand, the prediction given by MD displays a divergence at the 1:1 resonance, when the slave and master modes have close eigenfrequencies, $\omega _s \simeq \omega _p$. This divergence does not rely on a firm theoretical result from dynamical systems. Indeed, even though in the case a 1:1 internal resonance exists so that the two modes need to be taken into account to study the coupled dynamics, uncoupled solutions still exist and the backbone curves of these uncoupled solutions can be computed, thus preserving the meaning of the $\varGamma$ coefficients defined in Eq. (42), see, e.g. [13, 31, 56]. Thus, the divergence of ${\mathcal {C}}_{\text {MD}} / {\mathcal {C}}_{\text {SC}}$ is interpreted as a failure of the method. Finally, for small values of $\rho$, one can observe that the SMD method shows a singular behaviour and will predict unreasonably stiff behaviour. On the other hand, MD method gives a finite value, which is a bit different from the correct one given by normal form approach. All these results underline that MD and SMD can be used safely only when the assumption $\omega _s > 4\omega _p$ is fulfilled; otherwise, unreliable predictions may be given by these two methods.

2.4.2 Drift and mode shapes

A second comparison on the global outcomes of the three method can be provided by contrasting the mode shape dependence on amplitude. Indeed, assuming a single mode motion with $R_s = 0$ for all $s\ne p$ (only master mode p participates to the vibration), allows recovering the amplitude dependence of the pth mode shape. At small amplitude, the three methods recover the usual eigenmode, but they then differ in the way they are taking into account the cross-couplings with slave modes. Let us denote as $\mathbf {u}_{\text {MD}}$, $\mathbf {u}_{\text {SMD}}$ and $\mathbf {u}_{\text {NF}}$ the physical displacement following single-mode motion for each of the three methods. Using the previous formula allows one to reconstruct

$$\begin{aligned} \mathbf {u}_{\text {MD}}(t)&=\mathbf {\phi }_p R_p(t) - \underset{s\ne p}{\sum _{s=1}^N } \frac{g^s_{pp}}{\omega _s^2-\omega _p^2} R_p^2(t) \mathbf {\phi }_s, \end{aligned}$$

(46)

$$\begin{aligned} \mathbf {u}_{\text {SMD}}(t)&=\mathbf {\phi }_p R_p(t) - \mathbf {\phi }_p \frac{g^p_{pp}}{\omega _p^2} R_p^2(t) \nonumber \\&\quad - \underset{s\ne p}{\sum _{s=1}^N} \frac{g^s_{pp}}{\omega _s^2} R_p^2(t) \mathbf {\phi }_s, \end{aligned}$$

(47)

$$\begin{aligned} \mathbf {u}_{\text {NF}}(t)&=\mathbf {\phi }_p R_p(t) - \mathbf {\phi }_p \frac{g^p_{pp}}{\omega _p^2} \frac{1}{3} \left( R_p^2(t) + \frac{2}{\omega _p^2}{\dot{R}}_p^2(t) \right) \nonumber \\&\quad - \underset{s\ne p}{\sum _{s=1}^N} \frac{g^s_{pp}}{\omega _s^2} \left( \frac{\omega _s^2-2\omega _p^2}{\omega _s^2-4\omega _p^2} R_p^2(t)\right. \nonumber \\&\quad \left. - \frac{2}{\omega _s^2-4\omega _p^2} {\dot{R}}_p^2(t) \right) \mathbf {\phi }_s. \end{aligned}$$

(48)

Comparing the mode shapes given by MD and SMD, one can already underline that the summed term given by MD reduces to that given by SMD if one considers the slow/fast assumption with $\omega _s \gg \omega _p$. However, a difference persists in the two methods since with SMD an added quadratic term, depending on mode p only, is present (second term in (47)). This comes again from the treatment of the self-quadratic $g^p_{pp}$ term in Eq. (23), already underlined in Sect. 2.3.2. Indeed, the $g^p_{pp}$ term for the MD method is not present in the reconstruction, but in the dependence of the nonlinear frequency with amplitude instead, while the SMD method distributes the influence of this $g^p_{pp}$ term on the spatial reconstruction, but not on the amplitude–frequency relationship. This explains why the prediction of the hardening/softening behaviour appears to be more general for the MD method than for the SMD. Comparing now with the normal form approach, one can see that NF reduction gives rise to velocity-dependent terms in these formula, a feature that is not present in the other method, which is a direct consequence of the fact that NF method takes into account both independent displacement and velocity variables as it should be from a dynamical system perspective.

Again, one can also observe that the summed term in (48) reduces (at first significant order) to that provided by SMD when the slow/fast assumption is at hand, showing that the SMD method provides the most simplified expressions.

From the general expressions given in (46)–(48), one can isolate the constant term (zeroth harmonic) which is produced by the quadratic nonlinearity, in order to compare more closely one term of this general expansion. This constant term is known as a drift since it corresponds to the fact that due to quadratic nonlinearity, the oscillations are no more centred around zero, and it has already been compared for different reduction methods, see, e.g. [35, 57]. One can then simply replace $R_p(t)$ by the expression given by Eq. (40), while the values of $R_p^2(t)$ and ${\dot{R}}_p^2(t)$ up to second order write:

$$\begin{aligned} R_p^2(t)=&\frac{a_0^2}{2}\left( 1+\cos [2\omega _{\text {NL}}t]\right) +{\mathcal {O}}\left( a_0^3\right) \end{aligned}$$

(49)

$$\begin{aligned} {\dot{R}}_p^2(t)=&\frac{a_0^2}{2}\omega ^2_{\text {NL}}\left( 1-\cos [2\omega _{\text {NL}}t]\right) +{\mathcal {O}}\left( a_0^3\right) , \end{aligned}$$

(50)

where the nonlinear frequency $\omega _{\text {NL}}=\omega _p( 1 + a_0^2 \varGamma )$ has been introduced. Isolating the constant term leads to the following expressions for the drift d predicted by each reduction method:

$$\begin{aligned} \mathbf {d}_{\text {MD}}=&\frac{a_0^2}{2} \left( - \frac{g^p_{pp}}{\omega _p^2}\mathbf {\phi }_p - \underset{s\ne p}{\sum _{s=1}^n } \frac{g^s_{pp}}{\omega _s^2-\omega _p^2} \mathbf {\phi }_s \right) , \end{aligned}$$

(51)

$$\begin{aligned} \mathbf {d}_{\text {SMD}}=&\frac{a_0^2}{2} \left( - \frac{g^p_{pp}}{\omega _p^2}\mathbf {\phi }_p - \underset{s\ne p}{\sum _{s=1}^n } \frac{g^s_{pp}}{\omega _s^2} \mathbf {\phi }_s \right) , \end{aligned}$$

(52)

$$\begin{aligned} \mathbf {d}_{\text {NF}}=&\frac{a_0^2}{2} \left( - \frac{g^p_{pp}}{\omega _p^2} \left( \frac{1}{3}+\frac{2}{3}\frac{\omega ^2_{\text {NL}}}{\omega _p^2} \right) \mathbf {\phi }_p \right. \nonumber \\&\left. - \underset{s\ne p}{\sum _{s=1}^n } \frac{g^s_{pp}}{\omega _s^2} \left( 1-\frac{2\left( \omega ^2_{\text {NL}}-\omega _p^2\right) }{\omega _s^2-4\omega _p^2} \right) \mathbf {\phi }_s \right) . \end{aligned}$$

(53)

One can observe that assuming slow/fast dynamics, the drift predicted by MD reduces to that given by SMD. On the other hand, one can also see that in order to retrieve the drift predicted by SMD from $\mathbf {d}_{\text {NF}}$, one has to assume that the deviation of the nonlinear frequency $\omega _{\text {NL}}$ is small as compared to the linear frequency so that $\omega _{\text {NL}} \simeq \omega _p$. Hence, the prediction of the mode shape dependence on amplitude given by SMD is reliable only in the case where the backbone curve does not depart severely from the linear resonance, which is a strong assumption.

In order to point out a last difference on the theoretical expressions which will have important consequences in the next sections, let us also follow the first harmonic of the solution in the reconstruction procedure. Using Eq. (40) to define the harmonic content of the master variable, and going back to the harmonic content of the modal coordinates $X_i$ defined using either the QM method, Eq. (33), or the normal form approach, Eqs. (11)–(12), one can easily follow the first harmonic and retrieve its expression in the modal coordinates. Since p is the master mode and at lowest order $X_p=R_p$, then the most important contribution is present in $X_p$ as compared to other $X_k$’s. Let us denote as $\left[ X_p^{(H1)}\right] _{\text {MD}}$ the first harmonic for the MD approach (and SMD and NF for the other two methods); these expressions write:

$$\begin{aligned} \left[ X_p^{(H1)}\right] _{\text {MD}}&= a_0 \cos (\omega _{\text {NL}} t) \left( 1+{\mathcal {O}}\left( a_0^4\right) \right) , \end{aligned}$$

(54a)

$$\begin{aligned} \left[ X^{(H1)}_p\right] _{\text {SMD}}&= a_0 \cos (\omega _{\text {NL}} t)\nonumber \\&\quad \left( 1-a_0^2\;\dfrac{2}{3}\left( \dfrac{g^p_{pp}}{\omega _p^2}\right) ^2 +{\mathcal {O}}\left( a_0^4\right) \right) , \end{aligned}$$

(54b)

$$\begin{aligned} \left[ X^{(H1)}_p\right] _{\text {NF}}&= a_0 \cos (\omega _{\text {NL}} t) \left( 1+{\mathcal {O}}(a_0^4)\right) . \end{aligned}$$

(54c)

They underline the importance of the treatment of the self-quadratic $g^p_{pp}$ term between MD and SMD method. Indeed, whereas the amplitude $a_0$ defined from (40) corresponds, for the MD and NF cases, to the amplitude of the first harmonic in $X_p$, this is not the case for the QM derived from SMD. In that case, the amplitude has an extra term implying the self-quadratic coupling term. Importantly, this term appears as a difference so that the amplitude of the first harmonic can tend to small values with increasing $a_0$. Whereas all the comparisons led in this section show that the methods tend to be equivalent under a slow/fast assumption, this last expression highlights the fact that, for the SMD method, the amplitude of the master mode can be very different from the amplitude of the initial coordinate. The consequence of this finding will be more clearly illustrated in the next sections on examples and will be key to understand why the SMD method can fail even under the slow/fast assumption.

3 Comparison on two-degree-of-freedom systems

In this section, the comparisons drawn out on the theoretical expressions are illustrated on two-dof systems, in order to highlight the main differences on simple cases. Two different models are selected. The first one is derived from the equations of motion of a beam and is selected in order to mimic the nonlinearities present in a flat symmetric system, where these simplifying assumptions help in letting the methods based on SMD work properly. The second example has important quadratic couplings and better accounts from the problems arising with curved structures such as arches and shells.

3.1 A two-dof model representing a flat symmetric structure

3.1.1 Presentation of the model

The particular nature of the nonlinear couplings in the case of flat symmetric structures such as beams and plates relies on the simplifying facts that flexural and in-plane modes are linearly uncoupled and their nonlinear couplings involve simple terms that can be easily traced from the von Kármán models. These simplifications have been used in numerous recent papers in order to explain why a number of methods for producing ROMs are able to predict very good results in this case, see, e.g. [12, 19, 60, 61]. In order to propose a simple two-dof system mimicking these particular relationships, the von Kármán model for slender beams is used and simplified to two vibration modes, one flexural and one longitudinal, in order to produce the simplified system, from which the coefficients can be related to meaningful quantities of the beam and in particular to its slenderness.

A non-prestressed beam of length L is thus considered, with a uniform rectangular cross section of area $S=bh$ (h being the thickness and b the width) and second moment of area $I=bh^3/12$, made in an homogeneous and isotropic material of Young’s modulus E and density $\delta$. Boundary conditions are clamped at $X=0$ and $X=L$.

The equations of motion for the transverse displacement W(X, T) and the longitudinal displacement U(X, T) (X and T being the dimensional space and time variables), assuming von Kármán theory, read [12, 36]:

$$\begin{aligned}&{\ddot{W}}+\dfrac{E I}{\delta S} W^{''''}-\dfrac{E}{\delta } \left( U^{'}W^{'}+\dfrac{1}{2}W^{'\,3}\right) ^{'}=0, \end{aligned}$$

(55a)

$$\begin{aligned}&{\ddot{U}} - \dfrac{E}{\delta } (U{''}+W^{'}W^{''})=0. \end{aligned}$$

(55b)

A particular feature of Eq. (55) is that the longitudinal displacements are only quadratically coupled with the transverse, as shown in (55b). On the other hand, the only nonlinear terms appearing on the equations of motion for the flexural term W are: (i) a quadratic coupling involving a product between one in-plane and one transverse component, and a cubic term with only transverse components, see Eq. (55a).

Following [12], the equations of motion can be made non-dimensional so that the resulting system depends only on two physically meaningful parameters: the slenderness ratio $\sigma =h/L$, and the wavelength $\beta$ appearing naturally when solving the eigenvalue problem. Indeed, focusing on the linear problem for the transverse motion, the eigenvalue problem $\phi ^{''''} = \omega ^2 \frac{\delta S}{EI} \phi$ is solved by using a combination of sine, cosine, hyperbolic sine and hyperbolic cosine functions of kx, with k-dimensional wavelength such that $k^4 = \frac{\delta S}{EI} \omega ^2$ and $\beta = kL$. Assuming clamped boundary conditions, the characteristic equation for $\beta$, from which the eigenfrequencies are deduced, reads: $\cos (\beta )\cosh (\beta )=1$. The reader is referred to “Appendix F” for the details of this classical derivation.

Introducing the thickness h as characteristic length, so that the non-dimensional displacements are as $w=W/h$ and $u=U/h$, normalising time using $t=T (\beta ^2/L^2\sqrt{EI/\delta S})$ and the space variable with the beam length, $x=X/L$; Eq. (55) is rewritten as follows:

$$\begin{aligned}&w_{,tt}+\dfrac{1}{\beta ^4} w_{,xxxx} -\dfrac{12}{\beta ^4 \sigma } \left( u_{,x}w_{,x}\right) _{,x} \nonumber \\&\quad -\dfrac{6}{\beta ^4} \left( {w_{,x}}^3\right) _{,x}=0, \end{aligned}$$

(56a)

$$\begin{aligned}&u_{,tt} - \frac{12}{\beta ^4 \sigma ^2}u_{,xx} -\frac{12}{\beta ^4 \sigma }w_{,x}w_{,xx}=0. \end{aligned}$$

(56b)

In order to derive a minimal two-dof system from these equations, we select the first flexural eigenmode and the most important longitudinal mode coupled to the first flexural. From previous studies, see, e.g. [12, 49, 61], it is known that the fourth in-plane mode is strongly coupled to the first flexural. Let us denote as $q_1$ the modal amplitude of the first transverse mode and $p_4$ the modal amplitude of the fourth in-plane mode (see “Appendix F” for the details). Using a standard Galerkin projection (see, e.g. [12]), Eq. (56) can be rewritten as

$$\begin{aligned}&\ddot{q}_1+ q_1+ \dfrac{2G}{\sigma }p_4 q_1 + D q_1^3=0, \end{aligned}$$

(57a)

$$\begin{aligned}&\ddot{p}_4+\dfrac{(4\pi )^212}{\beta ^4 \sigma ^2} p_4 + \dfrac{G}{\sigma }q_1^2=0, \end{aligned}$$

(57b)

where D and G are the nonlinear coupling coefficients arising from the Galerkin projection and involve integral on the length of products of derivatives of the mode shape functions, see [12] for the general calculation and “Appendix F” for the detailed expression of these two coefficients. One can note, in particular, that due to the choice of the non-dimensional time to arrive at Eq. (56), the eigenfrequency of the first flexural mode is 1, while the natural frequency of the fourth in-plane mode reads $\omega _2^2 = \dfrac{(4\pi )^212}{\beta ^4 \sigma ^2}$. Due to the normalisation selected (involving $\omega _1 = 1$ for the fundamental mode), the term in factor of $p_4$ in Eq. (57a) can be easily interpreted as the square of the ratio $\rho =\omega _2/\omega _1$, recovering the term introduced in Sect. 2.4.1. Thanks to its explicit expression, $\rho$ can now be directly related to the slenderness ratio:

$$\begin{aligned} \rho =\dfrac{4\pi \;\sqrt{12}}{\beta ^2}\;\dfrac{1}{\sigma }\approx 1.95 \dfrac{1}{\sigma }, \end{aligned}$$

(58)

so that the final two-dof system that will be used for the investigations reads:

$$\begin{aligned}&{\ddot{X}}_1+ X_1+ 2\,{\bar{G}}\,\rho \,X_1 X_2 + D X_1^3=0, \end{aligned}$$

(59)

$$\begin{aligned}&{\ddot{X}}_2+\rho ^2 X_2 + {\bar{G}}\,\rho X_1^2=0, \end{aligned}$$

(60)

where ${\bar{G}}=G\,\beta ^2/(4\pi \sqrt{12})$ has been introduced for the ease of reading. Also the notation for the variables has been changed with $X_1=q_1$ and $X_2=p_4$ for the sake of simplicity. A particular feature of this system is that the coupling between master and slave mode is purely quadratic. Consequently, the potential third-order tensors from the normal form approach are all vanishing. In this case, the two nonlinear mappings are thus exactly at the same order due to the very simplified shape of the starting equations.

This system is now investigated in order to see how the methods under study behave when reducing the system to its first (flexural) mode using different nonlinear mappings. The advantage of this formulation is that all coefficients are related to a physical problem so that some insights can be given to the results obtained with this simplistic model with regard to continuous problems. In particular, Sect. 2.4 underlines that all methods show a convergence on some properties when a slow/fast assumption is assumed, which has been quantified precisely on the type of nonlinearity as occurring for $\rho > 4$. Also, divergent behaviours have been underlined and explained for $\rho \simeq 1$ (case of MD) and $\rho \simeq 2$ (case of normal form). Consequently, the system will be studied for three different values close to these points, namely $\rho =1.25$, $\rho =2.5$ and $\rho =10$. Note that a beam is generally considered as slender if $\sigma \le 1/20$. Thanks to Eq. (58), this means that $\rho \ge 40$. The consequence of this remark is that in all slender beams the slow/fast assumption is very well fulfilled, and our study concerns specific cases occurring for very thick beams. Regarding the nonlinear coefficients G and D, they only depend on the modes selected in the expansion. In our study, we will always consider the first flexural and fourth axial, so that G and D are constants since they only depend on the non-dimensional shape functions of the selected modes. In the remainder of the study, we have selected $D=2.67$, ${\bar{G}}=0.63$.

3.1.2 Results

The comparisons between the different methods are drawn out on the geometry of the manifolds, as well as on the dynamics onto these manifolds, described by the frequency–amplitude relationship (backbone curve). All the solutions are computed thanks to a numerical continuation method using the asymptotic–numerical method, implemented in the software MANLAB, where the unknowns are represented thanks to the harmonic balance method [11, 15, 29]. After a convergence study, the number of harmonics retained in the computations is 7. In each case, the master mode is the fundamental one, $X_1$, and the slave mode, $X_2$. The dynamics onto the reduced subspaces is given by Eq. (36) when using the MD approach, Eq. (37) if one considers SMD instead, and Eq. (38) with the normal form method, with $R_1$ the master coordinates. For the reduced models, continuation is performed on the master coordinate in order to compute the frequency–amplitude relationships. From these values, the nonlinear mappings, given either by Eq. (8) for the normal form approach or by Eq. (33) for the QM method, allow to retrieve the initial modal amplitude $X_1$ and $X_2$. From all these data, one can plot either the geometry of the manifolds in phase space $(X_1,Y_1,X_2,Y_2)$, or the backbone curves.

Fig. 2

Comparison of manifolds in phase space for the first example, and for three different values of $\rho =\omega _2/\omega _1$. The exact NNM, represented in violet (full system solution: FS), is compared to the reduction manifolds obtained by QM from MDs (dark orange), SMDs (yellow) and normal form (blue). a $\rho =1.25$, b $\rho =2.5$, c $\rho =10$ with slow/fast assumption fulfilled. (Color figure online)

Figure 2 shows the geometry of the manifolds obtained for this first system, when one increases the values of $\rho$ so as to meet the slow/fast assumption. One can remark that the reduced subspaces produced by the quadratic manifold method do not show a dependence on the velocity. Increasing the values of $\rho$ it is observed that the real manifold obtained from the full system loses this velocity dependence so that this approximation is less and less wrong. On the other hand, the manifold produced by normal form has two important advantages: it is an invariant manifold of the full system by construction, and it has this velocity dependence, hence allowing for a correct prediction of the reduction subspace, whatever the value of $\rho$. As a matter of fact the only limitation of the normal form approach is that it relies on a Taylor expansion, so that for large amplitudes, the solution departs from the exact manifold. But in any case, the correct invariant subspace is approximated. As already remarked, due to the fact that only quadratic couplings are present between master and slave coordinates, the manifolds shown in Fig. 2 for the normal form are obtained thanks to the second-order expansion, the third-order terms being all equal to zero.

Figure 2a shows also that the quadratic manifold produced by MD encounters a problem near the 1:1 resonance, which is here underlined since $\rho$ has been selected close to 1. Comparison with a full-order solution clearly shows that this is a failure of the method. On the other hand, Fig. 2c shows that when the slow/fast assumption is verified, all methods converge to the same reduced subspace, in line with the theoretical results.

We now turn to the prediction given on the backbone curves. First of all, one can compare the values of the $\varGamma$ coefficients dictating the type of nonlinearity. Equation (42) has thus been rewritten for the present two-dof system and now read, as a function of the ratio $\rho =\omega _2/\omega _1$:

$$\begin{aligned} \varGamma _{\text {MD}}&= \frac{3D}{8}- \frac{{\bar{G}}^2(3\rho ^2-2)\rho ^2}{4(\rho ^2-1)^2}, \end{aligned}$$

(61a)

$$\begin{aligned} \varGamma _{\text {SMD}}&= \frac{3D}{8}- \frac{{\bar{G}}^2(3\rho ^2+4)}{4\rho ^2}, \end{aligned}$$

(61b)

$$\begin{aligned} \varGamma _{\text {NF}}&= \frac{3D}{8}- \frac{{\bar{G}}^2(3\rho ^2-8)}{4(\rho ^2-4)}. \end{aligned}$$

(61c)

Fig. 3

Values of the coefficient $\varGamma$ dictating the hardening/softening behaviour for the first two-dof system. Comparison of $\varGamma _{\text {MD}}$, $\varGamma _{\text {SMD}}$ and $\varGamma _{\text {NF}}$, given, respectively, by QM with MDs, with SMDs, and normal form, Eq. (61), and for varying $\rho =\omega _2/\omega _1$ ratio

These values are represented in Fig. 3, which shows important similarities with Fig. 1. Again the same divergent behaviours are observed, and the convergence of all methods for $\rho > 4$ is clearly observed. To be more quantitative, the relative difference between $\varGamma _{\text {MD}}$ and $\varGamma _{\text {NF}}$ is 5% for $\rho =2.95$ and 1% for $\rho =3.93$. On the other hand, the difference between $\varGamma _{\text {SMD}}$ and $\varGamma _{\text {NF}}$ is 5% for $\rho =3.06$, and 1% for $\rho =4.18$, underlining clearly that $\rho \in [3,4]$ has to be understood as a transition zone. For very small values of $\rho$, the quadratic manifold based on SMD will predict incorrect result with a softening behaviour. Also, after its failure at $\rho =1$, the MD method will also produce an incorrect prediction with a softening behaviour.

Fig. 4

First-mode backbone curves as a function of modal amplitudes $X_1$ (first row) and $X_2$ (second row) and for different values of $\rho = \omega _2/\omega _1$. Comparisons between the exact solution (FS: full system, violet), that predicted by QM with MDs (dark orange), SMDs (yellow) and normal form (NF, blue). (Color figure online)

Figure 4 shows the backbone curves obtained from the reduced-order dynamics and compared to that obtained from the full system. The comparison is drawn on the main modal amplitude $X_1$, which shows the largest values (first row), but also on the slave coordinate $X_2$ (second row). The first case selected, just after the 1:1 resonance with $\rho =1.25$, shows, as envisioned in Fig. 3, that the QM produced from MD can be very wrong in this case and predict at first order a softening behaviour. When $\rho =2.5$, the three methods predict a very similar behaviour and are almost undistinguishable. One can note that for large amplitude, the full system solution is less and less hardening. This is probably a consequence of the vicinity of the 2:1 internal resonance. Since $\omega _2 = 2.5\omega _1$ and the behaviour is hardening, the nonlinear frequency tends to approach the 2:1 ratio at higher amplitudes, which could explain this particular behaviour of the full system solution. Finally, for $\rho =10$, the three methods give the same predictions which are fully aligned with the full system.

The conclusion on this first example with simple nonlinearities is in the line of the theoretical results, since all methods tend to perform well in the limit of the slow/fast assumption, again estimated as a ratio of 4 between the eigenfrequencies of the master and slave mode. On the other hand, when this assumption is not fulfilled, the quadratic manifold is not reliable and can produce incorrect predictions, in contrary to the normal form approach, that gives a correct ROM up to the third order, whatever the link between slave and master coordinates. These results explain also why the application of modal derivatives on slender structures that are flat and symmetric produce accurate results. Indeed, slenderness is fulfilled when $\rho$ is larger than 40, and our numerical experiments show that the slow/fast assumption can be considered as valid as soon as $\rho > 4$.

Fig. 5

Comparison of manifolds in phase space for the second two-dof example, and for three different values of $\rho =\omega _2/\omega _1$. The exact NNM, represented in violet (full system solution: FS), is compared to the reduction manifolds obtained by QM from MDs (dark orange), SMDs (yellow), and normal form up to the second order (blue manifold in the first line, plots a–c) and third order (green manifolds, second line in plots d–f) are given. a–d $\rho =1.25$, b–e $\rho =2.5$, c–f $\rho =10$ with slow/fast assumption fulfilled. In all cases, $\omega _1=1$. (Color figure online)

3.2 A two-dof model representative of a shell structure

3.2.1 Equations of motion

In this section, a system composed of a mass connected to two springs representing geometric nonlinearity is selected. This system has been used in a number of studies so that numerous results are already present in the literature; the interested reader is referred to [57] for the derivation of the equation of motions specifying the behaviour of the springs, and to [9, 42, 53, 57] for different results already published on this example system. The equations of motion read:

$$\begin{aligned} \begin{aligned}&{\ddot{X}}_1+\omega _1^2X_1+\frac{\omega _1^2}{2}\left( 3X_1^2+X_2^2\right) + \omega _2^2X_1X_2\\&\quad +\frac{\omega _1^2+\omega _2^2}{2}X_1\left( X_1^2+X_2^2\right) =0,\\&{\ddot{X}}_2+\omega _2^2X_2+\frac{\omega _2^2}{2}\left( 3X_2^2+X_1^2\right) +\omega _1^2X_1X_2\\&\quad +\frac{\omega _1^2+\omega _2^2}{2}X_2\left( X_1^2+X_2^2\right) =0.\\ \end{aligned} \end{aligned}$$

(62)

As compared to the previous example, this system has all quadratic nonlinear terms present in the equations of motion, and all the nonlinear coefficients are expressed directly from the two eigenvalues $\omega _1$ and $\omega _2$, so that the problem has only two parameters. Note that this model is not derived from a continuous shell structure like the previous example which was derived from the von Kármán beam equations; however, it is known that curved structures display strong quadratic couplings that are found in this system. Moreover, the results will show that this model is sufficient to show important departures between the three tested methods, which are due to the way the quadratic terms are processed.

3.2.2 Results

As for the preceding example, comparisons are drawn out on the geometry of the manifolds and the backbone curves. Numerical continuation is used to solve out the different systems and compare their outcomes. The eigenfrequency ratio $\rho =\omega _2/\omega _1$ is also used and the same values, namely 1.25, 2.5 and 10, are selected to observe the differences between the methods when tending to fulfil the slow/fast assumption. In the computation, $\omega _1=1$ in all cases so that one simply has $\omega _2=\rho$.

Fig. 6

First-mode invariant manifolds cut on the $Y_1=0$ plane, evaluated with the quadratic manifold method (QM) (either with MD in dark orange and SMD in yellow) and normal form (NF) approach, where the distinction between NF up to second order (blue line) and third order (dashed green line) is reported, and compared to the numerical solution obtained with the full system (FS). In all cases $\omega _1=1$. (Color figure online)

Figure 5 shows the geometry of the manifolds in phase space, as compared to the exact invariant manifold defining the first NNM of the system. The comment on the velocity dependence, already raised in the previous example, still holds: While for small values of $\rho$ the quadratic manifolds are not able to catch the correct curvature in this direction, for large values of $\rho$ the velocity dependence vanishes. Note that in all the three figures, the manifold produced by the SMD method has a smaller range in amplitude. This maximal range used for the representation has been fixed from the frequency–amplitude relationships (Fig. 7, when the nonlinear frequency has decreased of ten per cent and reaches the value 0.9—a softening behaviour is at hand in the considered cases), so that all manifolds span the same frequency range, but correspond to different amplitudes. This underlines, in particular, that even if the correct manifold is approximated, which is the case for $\rho =10$, the amplitude–frequency relationship may be not.

Since the only difference between second- and third-order normal form can be appreciated from the nonlinear mapping and not the reduced dynamics, Fig. 5 illustrates the case. In the first line, the manifold produced by the second-order normal form (in blue) is contrasted to the other methods, while the third-order normal form is shown in the second line (in green). One can observe that the effect of retaining the cubic term is especially important for the smallest values of $\rho =1.25$, but then the differences between second and third orders are barely visible. Interestingly, this example also shows that the quadratic manifolds produced by MD and SMD do not tend to the same geometries, even under the assumption of slow/fast dynamics. This is appreciated in Fig. 5, but is more clearly evidenced in Fig. 6 where a section of the manifolds in space $(X_1,X_2)$ is shown, without the amplitude limit given by the frequency, used in the 3D plot.

Unlike Fig. 5, Fig. 6 is directly obtained from the manifolds expressions given by Eq. (33) for the QM approach, and (11)–(12) for the normal form approach, by simply prescribing the values of $R_1$ and compute the resulting $(X_1,X_2)$ values. More specifically, let us underline the main difference between the MD and SMD method in this case. Substituting Eq. (33) with (23), the reconstruction of $(X_1,X_2)$ from the QM method derived from MD reads:

$$\begin{aligned} X_1&=R_1, \end{aligned}$$

(63a)

$$\begin{aligned} X_2&=-\frac{g^2_{11}}{\omega _2^2 - \omega _1^2}R_1^2 = - \frac{\omega _2^2}{2\left( \omega _2^2 - \omega _1^2\right) }R_1^2. \end{aligned}$$

(63b)

On the other hand, using SMD in the QM leads to:

$$\begin{aligned} X_1&=R_1 - \frac{g^1_{11}}{\omega _1^2}R_1^2 = R_1 - \frac{3}{2}R_1^2, \end{aligned}$$

(64a)

$$\begin{aligned} X_2&=-\frac{g^2_{11}}{\omega _2^2 }R_1^2 = - \frac{1}{2}R_1^2. \end{aligned}$$

(64b)

One can first notice that for this specific example, the manifold produced with the SMD method does not depend on the parameters $(\omega _1,\omega _2)$. Consequently, the cut of this manifold in $(X_1,X_2)$ plane in Fig. 6a–c for different values of $\rho$, is always the same. The second comment is on the slow/fast approximation: even though the value given for $X_2$ tends to be the same under the slow/fast assumption $\omega _2 \gg \omega _1$, this is not the case for $X_1$. This is a major difference between the two methods, so that a persistent error on the manifold is done when using SMD, whereas MD tends to the solution provided by the NF and full system when $\rho$ increases. The last interesting comment is on the fact that the manifold produced by SMD shows a constant folding point. Indeed, $X_1$ from Eq. (64a) cannot exceed the value of 1/6 (achieved at $R_1=1/3$) after which the quadratic term in Eq. (64a) is larger than the linear one.

Fig. 7

First-mode backbone curves for the second two-dof example with quadratic nonlinearity, as a function of modal amplitude $X_1$ (first row), and $X_2$ (second row), and for different values of $\rho = \omega _2/\omega _1$. Comparisons between the exact solution (FS: full system, violet) predicted by QM with MDs (dark orange), SMDs (yellow) and normal form (NF, blue, NF third order, dashed green). (Color figure online)

This is a direct consequence of the different treatment of the self-quadratic coupling term $g^1_{11}$, already underlined at the end of Sect. 2.4.2, leading to the fact that even under the slow/fast assumption, the QM built on SMD can lead to erroneous results. This point is further commented on the backbone curves comparison. First, Eq. (42) is written for this specific system, leading to the following predictions, as a function of the ratio $\rho =\omega _2/\omega _1$:

$$\begin{aligned} \varGamma _{\text {MD}}&= -\frac{16\rho ^4 -27\rho ^2+12}{16(\rho ^2-1)^2}, \end{aligned}$$

(65)

$$\begin{aligned} \varGamma _{\text {SMD}}&= -1, \end{aligned}$$

(66)

$$\begin{aligned} \varGamma _{\text {NF}}&= -\frac{\rho ^2-3}{\rho ^2-4}. \end{aligned}$$

(67)

In line with the constant manifold found with SMD, the method also predicts a constant type of nonlinearity, independent of the variations of the eigenfrequencies $(\omega _1,\omega _2)$. Assuming slow/fast partition, $\rho \rightarrow \infty$, then all three methods tend to predict the same $\varGamma$ coefficient dictating the hardening/softening behaviour. However, as underlined at the end of Sect. 2.4.2, the amplitude of the first harmonics for each method is different. Since in this case $g^1_{11}\ne 0$, a direct consequence of (54b) is that the backbone curves for the SMD method will show a saturation effect, the amplitude $X_1$ being unable to overcome a maximum value. This phenomenon is clearly visible in Fig. 7, depicting the backbone curves obtained for the three selected values of $\rho$. The constant value of $\varGamma _{\text {SMD}}$ has for direct consequence that the backbone predicted by the SMD quadratic manifold is almost unchanged with respect to variations of $\rho$. When the slow/fast assumption is fulfilled for $\rho =10$, as shown in Fig. 7c, the backbone predicted by SMD QM is in line with those predicted by the other methods at small amplitude level. However, at higher amplitude the SMD backbone moves away from the others and saturates to a limit value for all cases, since the amplitude is differently computed as shown in Eq. (54b). On the other hand, the backbone predicted by the MD method tends to the correct values under the slow/fast approximation, while the normal form approach always produces a correct prediction. More specifically, the prediction for the master $X_1$ component given by the normal form is the same if one considers a quadratic or cubic normal form expansion, see Eq. (11). On the other hand, the slave component $X_2$ is affected by the order and this is illustrated in Fig. 7d–f, where one can observe that, as for the manifold approximation in phase space, the third-order terms bring about a better estimate.

Fig. 8

Mesh and deformation along the mode under study for three different test cases: a a flat beam, b a shallow arch and c a non-shallow arch. For each test case, the mode under study corresponds to the first flexural mode in the plane $y-z$ where the curvature is imposed if present

4 Comparison on continuous structures

4.1 Presentation of the test cases

This section aims at drawing a comparison between the different methods when applied on a continuous structure discretised with three-dimensional finite elements. In order to investigate how the results obtained in the previous section are confirmed in the general case, three beams are considered and shown in Fig. 8. They have been selected in order to fulfil different assumptions that have been highlighted on the two-dof examples in order to achieve correct predictions from the ROMs. The first case, as shown in Fig. 8a, is a slender flat symmetric beam. The two other examples, as shown in Fig. 8b and c, are arches: the first one is shallow, while the third one is non-shallow. Adding curvature has two important effects. First, flexural and in-plane modes are no longer linearly uncoupled. Second, the curvature renders the restoring force asymmetric and an important quadratic nonlinearity appears between the bending modes. This example illustrates the fact that the slow/fast assumption is not enough to guarantee that the method based on static modal derivatives will converge. The curvature will be used in order to play it on the slow/fast assumption as well as on the values of the quadratic coupling terms.

In all three cases, the boundary conditions are clamped; the material parameters are selected as homogeneous linear elastic with Young modulus $E =124$ GPa, Poisson ratio $\nu = 0.3$ and density $\delta = 4400~\hbox {kg m}/\text {s}^2$. In each cases, an equal thickness h and width b are selected: $h=b=5~\hbox {cm}$. For the flat beam, the length is $L=0.7$ m. The arches have been built from a portion of a circle. For the shallow arch, the radius of curvature is set as 250 cm, for an angular span of $2\pi /15$, resulting in a curvilinear length of $20\pi /3 \simeq 1.05$ m. The height of the static deflection at centre is 5.5 cm, i.e. almost equal to the thickness. For the non-shallow arch, the radius of curvature is set as 50 cm, for an angular span of $2\pi /3$, resulting in the same curvilinear length of $20\pi /3 \simeq 1.05$ m, but with a static deflection of 25 cm, i.e. five times the thickness. All beams are discretised with three-dimensional hexahedral 20-node finite elements. The flat beam uses 60 elements (4 in the section and 15 in the length), resulting in a total number of 1287 dofs. The two arches have 96 solid elements (4 in the section and 24 in the length) and 2097 dofs. A relative coarse mesh has been selected in order to have a limited number of degrees of freedom so that all the methods can be handled easily. Indeed, the key point here is not to look for converged and refined results on a large frequency range, but to compare the different reduction methods on the same test examples. Moreover, as already shown in [61], using 3D elements leads to couplings with very high-frequency thickness modes, so that truncations and convergence are difficult to observe in general.

In the three cases, the nonlinear behaviour of the first flexural mode in the curvature plane is investigated. The mode shape is shown in Fig. 8. In the case of the flat beam, it corresponds to the first mode and its eigenfrequency is 545.60 Hz. As already underlined in Sect. 3.1, the most important coupling arises with the fourth in-plane mode, whose eigenfrequency is 15.19 kHz, so that the ratio $\rho$ between the most important slave mode and the master mode is in this case equal to 27.83. Consequently, the slow/fast assumption and our criterion $\rho \ge 4$ are perfectly fulfilled. This example can be seen as an extension of the first two-dof example, with the distinctive feature that now many more modes are coupled to the first bending, all of them being of higher frequencies than the fourth axial. Also, the nonlinear coupling terms have in this case a simplified form, following the general discussion given in Sect. 3.1. In the case of the arches, for the shallow arch the first flexural mode corresponds to the second mode of the structure and its eigenfrequency is equal to 372.28 Hz, and for the non-shallow arch, the first flexural mode corresponds to the fourth mode of the structure and its eigenfrequency is equal to 1003.99 Hz. Contrary to the case of flat symmetric structures, the curvature renders the restoring force asymmetric and an important quadratic nonlinearity appears between the bending modes. Investigating the important couplings between the linear modes of the curved beams shows that the first bending mode is strongly coupled with the third one. While the ratio between the first and third bending modes is 5.4 in the case of the flat beam, it decreases when the curvature increases. Consequently, for the case of the shallow arch, this ratio is equal to 3.44 (eigenfrequency of third bending equal to 1283.33 Hz), and 1.66 for the non-shallow arch (eigenfrequency of third bending equal to 1665.11 Hz). These two examples have thus been built as an extension of the second two-dof example. For the shallow arch, the slow/fast assumption is almost fulfilled (3.44 is slightly smaller than the criterion we proposed with a limit value at 4), but now important quadratic couplings are present and in particular the self-quadratic term $g^p_{pp}$. Finally, the case of the non-shallow arch allows testing a case where the slow/fast assumption does not hold, and important self-quadratic terms are present.

4.2 Amplitude–frequency relationships

The methods are compared on their ability to predict the backbone curves. A reference solution is computed thanks to a numerical continuation on all the degrees of freedom of the structure, using a code with parallel implementation of harmonic balance method and pseudo-arc-length continuation algorithm [8]. In this computation, a small amount of mass proportional damping is added under the form $\zeta \omega _p \mathbf {M}$ so that a frequency-response function (FRF) is computed, in the vicinity of the eigenfrequency of the master mode (first flexural). The values of $\zeta$ are 0.18%, 0.27% and 0.1% for the flat beam, shallow, and non-shallow arches, respectively. The forcing is located in the central node of each mesh in the y-direction in order to excite the first flexural mode. The force amplitude is chosen in order to have a displacement amplitude at resonance comparable to the thickness so that its values are 5 kN, 1.5 kN and 2.5 kN for the flat beam, shallow and non-shallow arches, respectively. It must be noticed that in the case of curved structures the value of amplitude of vibration equal to the thickness has not been achieved and the reported FRFs excite a maximum amplitude of approximately half of it. In fact, due to the long computational time that the full model FRF requires, approximately 1 day for each FRF, and due to its high chances to undergo internal resonance with higher modes, these values have been selected in order to stay in the limit of one-mode approximation without exciting more complex dynamics. However, with this level of amplitudes, the nonlinearity is sufficiently important so that its effect is clearly visible on the backbone curves.

The ROMS are built using QM or NNM approach, and their backbone curves are computed in the same manner than in the previous section, assuming a single master mode. For the normal form approach, the third-order coefficients have not been included in the computation. Indeed, the third-order tensors require the computation of huge number of coupling coefficients from the modal basis expressions, which would need for an important number of precomputation steps. This choice has also been guided by the fact that comparing the two methods at the same order of accuracy is more meaningful. The FRF of the ROMS has not been computed since taking into account the damping of the slave modes is important to achieve good results. If the normal form theory has been developed for that purpose, see, e.g. [53] where the effect of a small amount of damping of the slave modes on the FRF of the master mode is reported, the inclusion of the damping for the modal derivatives has not been derived theoretically yet. Hence, it appears that a better comparison is given on the backbone curves only, and the FRF of the full model with a small amount of damping is used to underline whether the nonlinearity is correctly addressed by the methods.

Fig. 9

Comparison of backbone curves obtained from QM with MDs (dark orange), SMDs (yellow) and normal form approach (blue), for the three tested structures: a flat beam, b shallow arch, c non-shallow arch. Non-dimensional amplitude of flexural displacement (along y, non-dimensionalised with respect to the thickness) of the central node of each beam as a function of $\omega /\omega _1$ where $\omega _1$ refers to the eigenfrequency of the first flexural mode studied shown in Fig. 8. The backbone curves are contrasted to the FRF obtained on the full system (FS, violet) with numerical continuation and a small amount of damping, see text. The vertical grey dashed lines represent the frequency at which the saddle-node bifurcation of the forced response occurs. (Color figure online)

Figure 9 shows the numerical results obtained for the three cases. The case of the flat beam is the one having the most assumptions fulfilled (slow/fast separation and no self-quadratic terms). Consequently, the three methods match very well and are all able to retrieve correctly the nonlinearity of the full model with a very good accuracy. In the case of the shallow arch, the slow/fast assumption is almost fulfilled (since being a little bit below the proposed criterion $\rho \ge 4$), and important self-quadratic coupling appears due to the curvature. The main consequence is that the QM built from SMD is not able anymore to predict the correct type of nonlinearity. As already found for the second two-dof example, it overpredicts the softening behaviour and makes appear again the saturation phenomenon in the amplitude of the backbone. On the other hand, both QMs based on MD and normal form methods give a correct prediction. For the non-shallow arch, the slow/fast assumption does not hold anymore. The consequence is that the MD method does not predict the correct nonlinearity. This example again illustrates clearly that: (i) as soon as important self-quadratic terms appear (case of arches and shells), then the SMD method is not reliable anymore, whatever the slow/fast assumption is fulfilled or not; (ii) the MD can still give correct result but only if the criterion $\rho \ge 4$ for the slow/fast assumption is fulfilled. As soon as $\rho$ gets under this value, then the solution starts departing from the full-order model and becomes unreliable when $\rho \le 2$.

4.3 Nonlinear mode shapes

The different approximations made by the three methods are finally contrasted on the mode shape dependence on amplitude, illustrating the equations given in Sect. 2.4.2. Recalling Eqs. (46)–(48), it is possible to see that, for each method, the contributions to the nonlinear mode shape can be divided into (i) a deformation along the master p mode and (ii) a deformation that contains all the coupled modes but the pth. In order to make the figures more illustrative, and since the amplitude of the deformation along pth mode generally gives the dominant contribution, it is decided to compare the outcomes of the methods only on (ii) part of the solution. Also, since the normal form approach constructs the solution both with displacements and velocities, to draw a better comparison the focus will be on the time step where the reduced variable $R_p(t)$ reaches its maximum and minimum values (i.e. a turning point such that ${\dot{R}}_p(t)=0$).

Fig. 10

Comparisons of the additional terms perturbing the linear mode shape solutions (deformation orthogonal to $\mathbf {\phi }_p$) for the case of the flat beam, computed at the saddle-node bifurcation point marked in Fig. 9a, fixing the frequency at which they have been computed. a Full model solution, representation of the axial component of displacement $\mathbf {u}^{\perp ,sym}_{\text {FS}}.\mathbf {e}_z$ of the centre line nodes, b normal form: $\mathbf {u}^{\perp }_{\text {NF}}(t^*).\mathbf {e}_z$, c modal derivative: $\mathbf {u}^{\perp }_{\text {MD}}(t^*).\mathbf {e}_z$, d static modal derivative $\mathbf {u}^{\perp }_{\text {SMD}}(t^*).\mathbf {e}_z$

Under these assumptions, let us define as $\mathbf {u}^{\perp }$ the component of the nonlinear mode shape $\mathbf {u}$ that is orthogonal to $\mathbf {\phi }_p$. From Eqs. (46)–(48), it reads, for the three different methods:

$$\begin{aligned} \mathbf {u}^{\perp }_{\text {MD}}(t^*)=&-R_p^2(t^*) \underset{s\ne p}{\sum _{s=1}^N } \frac{g^s_{pp}}{\omega _s^2-\omega _p^2} \mathbf {\phi }_s, \end{aligned}$$

(68a)

$$\begin{aligned} \mathbf {u}^{\perp }_{\text {SMD}}(t^*)=&-R_p^2(t^*) \underset{s\ne p}{\sum _{s=1}^N} \frac{g^s_{pp}}{\omega _s^2} \mathbf {\phi }_s, \end{aligned}$$

(68b)

$$\begin{aligned} \mathbf {u}^{\perp }_{\text {NF}}(t^*)=&-R_p^2(t^*) \underset{s\ne p}{\sum _{s=1}^N} \frac{g^s_{pp}}{\omega _s^2} \left( \frac{\omega _s^2-2\omega _p^2}{\omega _s^2-4\omega _p^2} \right) \mathbf {\phi }_s, \end{aligned}$$

(68c)

where $t^*$ is the time instant where $R_p$ is either maximum or minimum.

In order to compare to the full-order solution, the deformation must be first filtered out from its component along the pth mode. One can thus define $\mathbf {u}^{\perp }_{\text {FS}}$ for the full system as:

$$\begin{aligned} \mathbf {u}^{\perp }_{\text {FS}}(t)= \mathbf {u}_{\text {FS}}(t)-\dfrac{\mathbf {\phi }^\text {T}_p \mathbf {u}_{\text {FS}}(t)}{\mathbf {\phi }^\text {T}_p \mathbf {\phi }_p}\mathbf {\phi }_p. \end{aligned}$$

(69)

Finally, given the quadratic nature of the deformation computed from the reduction methods based on second-order expansions (and clearly underlined by the dependence in $R_p^2$ in Eq. (68)), the third-order component should be also filtered out from $\mathbf {u}^{\perp }_{\text {FS}}$ for a closer comparison. In order to cancel the odd harmonics of the full-order solution, we thus define $\mathbf {u}^{\perp ,sym}_{\text {FS}}$ as the symmetric part of $\mathbf {u}^{\perp }_{\text {FS}}$ with respect to amplitude:

$$\begin{aligned} \mathbf {u}^{\perp ,sym}_{\text {FS}}=\dfrac{1}{2}\left( \mathbf {u}^{\perp }_{\text {FS}}(t^{\mathrm{max}})+ \mathbf {u}^{\perp }_{\text {FS}}(t^{\mathrm{min}})\right) . \end{aligned}$$

(70)

This value will be used as reference and compared to the prediction of the ROMS given by Eq. (68).

Fig. 11

Comparisons of the additional terms perturbing the linear mode shape solutions (deformation orthogonal to $\mathbf {\phi }_p$) for the case of the shallow arch, computed at the saddle-node bifurcation point marked in Fig. 9b, fixing the frequency at which they have been computed. a Full model solution, representation of the displacement $\mathbf {u}^{\perp ,sym}_{\text {FS}}$ of the centre line in the $y-z$ plane; vertical axis $\mathbf {u}^{\perp ,sym}_{\text {FS}}.\mathbf {e}_y$ (transverse direction) and horizontal axis $\mathbf {u}^{\perp ,sym}_{\text {FS}}.\mathbf {e}_z$ (axial direction), b normal form: $\mathbf {u}^{\perp }_{\text {NF}}(t^*)$, c modal derivative: $\mathbf {u}^{\perp }_{\text {MD}}(t^*)$, d static modal derivative $\mathbf {u}^{\perp }_{\text {SMD}}(t^*)$. Grey lines: position of the centre line of the beam at rest. Solution amplified of factor 15

Fig. 12

Comparisons of the additional terms perturbing the linear mode shape solutions (deformation orthogonal to $\mathbf {\phi }_p$) for the case of the non-shallow arch, computed at the saddle-node bifurcation point marked in Fig. 9c, fixing the frequency at which they have been computed. a full model solution, representation of the displacement $\mathbf {u}^{\perp ,sym}_{\text {FS}}$ of the centre line in the $y-z$ plane; vertical axis $\mathbf {u}^{\perp ,sym}_{\text {FS}}.\mathbf {e}_y$ (transverse direction) and horizontal axis $\mathbf {u}^{\perp ,sym}_{\text {FS}}.\mathbf {e}_z$ (axial direction), b normal form: $\mathbf {u}^{\perp }_{\text {NF}}(t^*)$, c Modal derivative: $\mathbf {u}^{\perp }_{\text {MD}}(t^*)$, d Static modal derivative $\mathbf {u}^{\perp }_{\text {SMD}}(t^*)$. Grey lines: position of the centre line of the beam at rest. Solution amplified of factor 50

Figure 10 shows the comparison between the $\mathbf {u}^{\perp }$ defined by Eq. (70) for the full-order system and those produced by the reduced-order models, Eq. (68), for the case of the flat beam. Importantly enough, since the nonlinear couplings are with in-plane modes, the contributions of the $\mathbf {u}^{\perp }$ along the axial z direction are shown in Fig. 8, since the most important contributions are along this direction. As it could be awaited from the previous analyses, Fig. 10 clearly shows that the three ROMs are all able to recover the correct spatial dependence of the contributions of coupled modes to the fundamental flexural NNM. Also, this contribution is mostly conveyed by the fourth in-plane mode, being the most importantly coupled to the fundamental flexural mode. Note that the amplitude used to construct this figure is the one corresponding to the upper saddle-node bifurcation point in the FRF of the full-order system, as shown by the grey vertical line in Fig. 9a. At that point, the backbones and the FRF meet so that it can be used safely for a correct comparison. It also corresponds to an amplitude of one time the thickness for the mode shape.

In the case of the shallow arch, some differences are appearing due to the self-quadratic coupling term, creating a deficiency in the prediction given by the SMD. This is underlined in the nonlinear mode shape dependence in Fig. 11, where in this case, since the most important coupling is between bending modes, the contributions of the different $\mathbf {u}^{\perp }$ are represented along the transverse y direction. The amplitude used for the figure is illustrated in Fig. 9b with a grey line and still corresponds to the upper saddle-node bifurcation point in the FRF of the full-order system. One can observe in Fig. 11 that, in the line of the results found on the nonlinear amplitude–frequency relationships, normal form and MD methods are able to retrieve the correct spatial dependence for the contribution of the slave modes. On the other hand, the treatment of the self-quadratic term by the SMD approach prevents the correct prediction of this spatial dependence.

The case of the non-shallow arch is shown in Fig. 12, for an amplitude of motion marked by the grey line in Fig. 9c. Following the observation on the frequency, one can notice that only the normal form approach is able to retrieve the correct spatial dependence. On the other hand, SMD method fails because of the incorrect treatment of the self-quadratic term, while QM constructed from MD does not produce the correct result since the slow/fast assumption is not fulfilled anymore.

5 Conclusion

In this contribution, a detailed comparison of different methods proposed in the recent years in order to define nonlinear mappings with the aim of providing accurate reduced-order models for geometrically nonlinear structures, has been made. The quadratic manifold proposed from the definitions of modal derivatives has thus been contrasted to the normal form theory, related to the definition of nonlinear normal modes as invariant manifolds in phase space. While the quadratic manifold only contains the displacements as unknowns, the normal form approach takes into account displacements and velocities, thus giving a more complete link to the geometry in phase space. Secondly, the quadratic manifold is defined up to the second order, while current expressions of normal form are up to order three and can be continued to higher orders easily. Thirdly, normal form theory relies of firm mathematical theorems, ensuring a clean conceptual framework, while modal derivatives appear as an ad hoc, yet efficient, method used in the vibration community.

The main outcomes of this article are the following. First, the theoretical derivations of the quadratic manifold using either MD or SMD have been fully made explicit. These calculations have highlighted the fact that both methods do not handle the quadratic terms in the same manner, especially the self-quadratic coupling terms arising between the master coordinates. This difference has been found to have important consequences on the global predictions of the methods. Secondly, detailed comparisons between the three methods have been fully analysed on the mathematical expressions: nonlinear change of coordinates, reduced-order dynamics, and main predictive outcomes of the methods such as type of nonlinearity, drift and mode shape dependence on amplitude. To illustrate the results, two two-dof systems have been used as starting example, and the results found from these have been extended to a continuous structure: a clamped–clamped beam with varying curvatures.

A main result of our investigations is that the results predicted by the QM approach with MDs converge to those provided by the normal form approach, only in the case where a slow/fast assumption between master and slave coordinates holds. This result is fully in the line of general theorems provided in [17, 60] and thus further illustrates the general findings given in these papers where a more general framework including damping is given, together with an exact result that do not rely on asymptotic expansion. A first quantification of the limit value for the slow/fast assumption to hold has been provided, based on the predicted values for the type of nonlinearity, showing that a small gap is needed: $\omega _s > 4 \omega _p$, thus justifying a posteriori the good results found by previously published papers using this method. However, the different treatment of the quadratic nonlinearity (and more specifically the self-quadratic coupling term) between MD and SMD leads to the fact that even with a slow/fast assumption, the QM built from SMD can lead to erroneous predictions, as soon as an important self-quadratic coupling term is present. This result has important implications when one wants to build ROMs for slender curved structures such as arches and shells. This specific feature has been clearly highlighted on the two-dof system and found in the more general case of a non shallow arch. On the other hand, the robustness of the normal form approach has been underlined in each case.

These results argue for the use of the tools from dynamical system theory to derive safe and robust ROMS: invariant manifold, normal form theory and spectral submanifold. A limitation could be the use of these methods in the context of FE models where the need of computing, possibly in a non-intrusive manner, the nonlinear coefficients might be a difficult task, see, e.g. all the literature related to the STEP method (stiffness evaluation procedure, see, e.g. [12, 33, 34]). However, recent developments show that the coefficients can be directly computed, for the case of spectral submanifold [60], or for the case of normal form in a non-intrusive manner [62], so that this limitation does not hold anymore.

Appendices

Appendix A: Definition of dot products for tensors

The products between quadratic nonlinear tensors $\mathbf {G}$ and $\mathbf {g}$ and the displacement vectors $\mathbf {u}$ and $\mathbf {X}$, in physical and modal basis, respectively, are defined as follows:

$$\begin{aligned} \mathbf {G}\mathbf {u}\mathbf {u} =&\sum _{i=1}^N \sum _{j=1}^N \mathbf {G}_{ij}{u}_{j}{u}_{i}, \end{aligned}$$

(71a)

$$\begin{aligned} \mathbf {g}\mathbf {X}\mathbf {X} =&\sum _{i=1}^N \sum _{j=1}^N \mathbf {g}_{ij}{X}_{j}{X}_{i}, \end{aligned}$$

(71b)

where $\mathbf {G}_{ij}$ and $\mathbf {g}_{ij}$ are the N-dimensional vector of coefficients $G^p_{ij}$ and $g^p_{ij}$, for $p=1\ldots N$. Note that in the course of this paper, the summation is kept complete so that all the terms in Eqs. (71a)–(71b) are present. This is a choice of representation, the other choice (often realised in literature) consisting in symmetrising the tensor, given the fact that products ${u}_{j}{u}_{i}$ can commute. In this case, a quadratic tensor is made symmetric such that for example $\mathbf {G}_{ij}=0$ when $j<i$, so that the summations can be written for $j\ge i$ only. Here, we consider full tensors of coefficients without using their potential symmetry.

Similarly, the cubic nonlinear tensors $\mathbf {H}$ and $\mathbf {h}$, with current terms $H^p_{ijk}$ and $h^p_{ijk}$, contract to an N-dimensional vector when multiplied with three displacement vectors. The cubic terms thus explicitly write, in indicial components, with $\mathbf {u}$ and $\mathbf {X}$ the two associated displacement vectors:

$$\begin{aligned} \mathbf {H}\mathbf {u}\mathbf {u}\mathbf {u} =&\sum _{i=1}^N \sum _{j=1}^N \sum _{k=1}^N \mathbf {H}_{ijk}{u}_{k}{u}_{j}{u}_{i}, \end{aligned}$$

(72a)

$$\begin{aligned} \mathbf {h}\mathbf {X}\mathbf {X}\mathbf {X} =&\sum _{i=1}^N \sum _{j=1}^N \sum _{k=1}^N \mathbf {h}_{ijk}{X}_{k}{X}_{j}{X}_{i}, \end{aligned}$$

(72b)

where again $\mathbf {H}_{ijk}$ is the N-dimensional vector with entries $H^p_{ijk}$, for $p=1,\ldots ,N$. The innermost product defined above coincides with a matrix product performed on the last index of the tensors; a more extensive definition of this notation is provided:

$$\begin{aligned} (\mathbf {G} \mathbf {u})_i =&\sum _{j=1}^N \mathbf {G}_{ij} {u}_{j}, \end{aligned}$$

(73a)

$$\begin{aligned} (\mathbf {H} \mathbf {u})_{ij} =&\sum _{k=1}^N \mathbf {H}_{ijk} {u}_{k}, \end{aligned}$$

(73b)

$$\begin{aligned} (\mathbf {H} \mathbf {u} \mathbf {u})_i =&\sum _{j=1}^N \sum _{k=1}^N \mathbf {H}_{ijk} {u}_{j} {u}_{k}. \end{aligned}$$

(73c)

Same definition holds for the equivalent products in modal basis.

Appendix B: Normal form coefficients

In this appendix, the nonlinear coefficients of the mapping given by the normal form theory are given in detail. For the sake of brevity, only the coefficients appearing in the simplified case where a single master mode p is selected in the reduced model, are recalled. The interested reader can refer to [53, 57] for most complete expressions covering all the cases, including also damping. We begin with the second-order coefficients, $a^s_{pp}$, $b^s_{pp}$ and $\gamma ^s_{pp}$ coefficients, with p the master mode and s a slave mode :

$$\begin{aligned} a^s_{pp}&=g^s_{pp}\dfrac{2\,\omega ^2_p-\omega ^2_s}{-\omega ^2_s \left( 4\,\omega ^2_p-\omega ^2_s\right) }, \end{aligned}$$

(74a)

$$\begin{aligned} b^s_{pp}&=g^s_{pp}\dfrac{2}{-\omega ^2_s\left( 4\,\omega ^2_p-\omega ^2_s\right) }, \end{aligned}$$

(74b)

$$\begin{aligned} \gamma ^s_{pp}&=g^s_{pp}\dfrac{2}{4\,\omega ^2_p-\omega ^2_s}. \end{aligned}$$

(74c)

The third-order nonlinear mapping coefficients are equal to zero in the specific case when $s= p$. Their full expressions for $s\ne p$ read:

$$\begin{aligned} r^s_{ppp}&= \dfrac{ \left( A^s_{ppp}+h^s_{ppp}\right) \left( 7 \omega _p^2-\omega _s^2\right) + 2 B^s_{ppp}\left( \omega _p^4\right) }{\left( \omega _s^2-\omega _p^2\right) \left( \omega _s^2-9\omega _p^2\right) }, \end{aligned}$$

(75a)

$$\begin{aligned} u^s_{ppp}&= \dfrac{ 6\left( A^s_{ppp}+h^s_{ppp}\right) + B^s_{ppp}\left( 3\omega _p^2-\omega _s^2\right) }{\left( \omega _s^2-\omega _p^2\right) \left( \omega _s^2-9\omega _p^2\right) }, \end{aligned}$$

(75b)

$$\begin{aligned} \mu ^s_{ppp}&= \dfrac{ 6\left( A^s_{ppp}+h^s_{ppp}\right) + B^s_{ppp}\left( 3\omega _p^2-\omega _s^2\right) }{\left( \omega _s^2-\omega _p^2\right) \left( \omega _s^2-9\omega _p^2\right) }, \end{aligned}$$

(75c)

$$\begin{aligned} \nu ^s_{ppp}&= \dfrac{ 3\left( A^s_{ppp}+h^s_{ppp}\right) \left( 3\omega _p^2-\omega _s^2\right) + 2 B^s_{ppp}\left( \omega _p^2\omega _s^2\right) }{\left( \omega _s^2-\omega _p^2\right) \left( \omega _s^2-9\omega _p^2\right) }, \end{aligned}$$

(75d)

where

$$\begin{aligned} A^s_{ppp}&= \sum _{l=1}^N 2 g^s_{pl} a^l_{pp}, \end{aligned}$$

(76a)

$$\begin{aligned} B^s_{ppp}&= \sum _{l=1}^N 2 g^s_{pl} b^l_{pp}, \end{aligned}$$

(76b)

with $\mathbf {g}_{pl}$ being the vector of quadratic coupling between the master mode p and a generic mode l of the structure.

Appendix C: Linear change of coordinates from physical to modal basis

The nonlinear force vector in physical basis reads:

$$\begin{aligned} \mathbf {F}(\mathbf {u})= \mathbf {K}\mathbf {u}+ \mathbf {G} \mathbf {u} \mathbf {u} + \mathbf {H} \mathbf {u} \mathbf {u} \mathbf {u}. \end{aligned}$$

(77)

The nonlinear force vector in modal basis is in the form:

$$\begin{aligned} \mathbf {f}(\mathbf {X}) = \mathbf {\Omega }^2\mathbf {X}+ \mathbf {g} \mathbf {X} \mathbf {X} + \mathbf {h} \mathbf {X} \mathbf {X} \mathbf {X} , \end{aligned}$$

(78)

where the assumption of mass normalised eigenvectors is used to retrieve the squared eigenfrequencies on the diagonal of the matrix $\mathbf {\Omega }^2$.

The transformation from physical to modal basis uses the full linear eigenvector matrix $\mathbf {\Phi }$ and reads:

$$\begin{aligned} \mathbf {f}(\mathbf {X}) = \mathbf {\Phi }^T\mathbf {F}(\mathbf {\Phi }\mathbf {X}). \end{aligned}$$

(79)

Expanding the right-hand side (RHS) term, it reads:

$$\begin{aligned} \mathbf {f}(\mathbf {X})= & {} \mathbf {\Phi }^T\mathbf {K} \;\mathbf {\Phi }\mathbf {X} + \mathbf {\Phi }^T\mathbf {G} \;\mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X} \nonumber \\&+ \mathbf {\Phi }^T \mathbf {H}\; \mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X} . \end{aligned}$$

(80)

The relation between linear stiffness matrix in physical and modal coordinates can easily be found by comparing the linear terms in Eqs. (78) and (80), allowing one to retrieve the classical formula:

$$\begin{aligned} \mathbf {\Omega }^2=\mathbf {\Phi }^T\mathbf {K}\mathbf {\Phi }. \end{aligned}$$

(81)

To relate the quadratic and cubic tensors in physical basis to those in modal basis, it is necessary to expand the term $\mathbf {\Phi }\mathbf {X}$ into the sum of all eigenvectors multiplied by their modal amplitudes:

$$\begin{aligned} \mathbf {\Phi }\mathbf {X}=\sum ^N_{i=1} \mathbf {\phi }_i X_i, \end{aligned}$$

(82)

and substitute this sum into Eq. (80). By doing so, one obtains:

$$\begin{aligned}&\mathbf {\Phi }^T \mathbf {G}\;\mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X}\nonumber \\&\quad =\; \mathbf {\Phi }^T \mathbf {G} \;\sum _{i=1}^N \mathbf {\phi }_i X_i \; \sum _{j=1}^N \mathbf {\phi }_j X_j \nonumber \\&\quad = \sum _{i=1}^N \sum _{j=1}^N \mathbf {\Phi }^T \mathbf {G}\mathbf {\phi }_i \mathbf {\phi }_j \;X_i X_j, \end{aligned}$$

(83a)

$$\begin{aligned}&\mathbf {\Phi }^T\mathbf {H}\; \mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X}\; \mathbf {\Phi }\mathbf {X}\nonumber \\&\quad = =\; \mathbf {\Phi }^T \mathbf {H} \;\sum _{i=1}^N \mathbf {\phi }_i X_i \; \sum _{j=1}^N \mathbf {\phi }_j X_j \; \sum _{k=1}^N \mathbf {\phi }_k X_k \nonumber \\&\quad = \sum _{i=1}^N \sum _{j=1}^N \sum _{k=1}^N \mathbf {\Phi }^T \mathbf {H}\mathbf {\phi }_i \mathbf {\phi }_j\mathbf {\phi }_k \; X_i X_j X_k , \end{aligned}$$

(83b)

where this simplification comes from the rearrangement of summations made in order to isolate the modal amplitudes. Finally, comparing the RHS of Eq. (83) with the definition of products given in Eqs. (71b) and (72b) of “Appendix A” leads to:

$$\begin{aligned} \mathbf {g}_{ij}&=\mathbf {\Phi }^T \mathbf {G}\mathbf {\phi }_i \mathbf {\phi }_j, \end{aligned}$$

(84a)

$$\begin{aligned} \mathbf {h}_{ijk}&=\mathbf {\Phi }^T \mathbf {H}\mathbf {\phi }_i \mathbf {\phi }_j \mathbf {\phi }_k. \end{aligned}$$

(84b)

These two equations allow expressing the quadratic and cubic coefficients of the modal basis from those computed in the physical basis. Please note that the obtained formula directly depends on the choice of the representation used for the coefficients. Since we have selected to keep full-order tensors of coefficients without exploiting the symmetries arising from the fact that the usual product is commutative, the obtained formula are as in (84a)–(84b). If one chooses to use symmetric tensors for the coefficients, then, for the quadratic term, the relationship would have read $\mathbf {g}_{ij}=2\mathbf {\Phi }^T \mathbf {G}\mathbf {\phi }_i \mathbf {\phi }_j$ and $\mathbf {g}_{ji}=\mathbf {0}$.

Appendix D: First- and second-order derivatives of the nonlinear force vector

Given the definition of the nonlinear force tensor in physical basis, one can show that:

$$\begin{aligned} \left( \dfrac{\partial \mathbf {F}}{\partial \mathbf {u}}\right) ^r_s&= \dfrac{\partial F^r}{\partial u_s}= \text {K}_{rs} + \sum _{j=1}^N \; ({G}^r_{sj}+{G}^r_{js}) {u}_{j}\\&\quad + \sum _{j=1}^N \sum _{k=1}^N\; ({H}^r_{sjk}+{H}^r_{jsk}+{H}^r_{jks}) {u}_{j} {u}_{k}\\&=\text {K}^r_{s} + \sum _{j=1}^N \; \;2\;{G}^r_{sj} {u}_{j} \\&\quad + \sum _{j=1}^N \sum _{k=1}^N\; \;3\;{H}^r_{sjk} {u}_{j} {u}_{k}, \end{aligned}$$

where this simplification is derived from the symmetry of the quadratic and cubic tensors [34]. In compact form, we can write:

$$\begin{aligned} \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}}=\mathbf {K}+2\mathbf {Gu}+3\mathbf {Hu}\mathbf {u}, \end{aligned}$$

(85)

and similarly, for the second-order derivatives:

$$\begin{aligned} \left( \dfrac{\partial ^2 \mathbf {F}}{\partial \mathbf {u}\partial \mathbf {u}}\right) ^r_{sp}&= \dfrac{\partial F^r}{\partial u_s \partial u_p} = ({G}^r_{sp}+{G}^r_{ps}) \\&\quad + \sum _{k=1}^N\; ({H}^r_{spk}+{H}^r_{psk} +{H}^r_{pks}\\&\quad +{H}^r_{skp}+{H}^r_{ksp}+{H}^r_{kps}) {u}_{k}\\&= 2\;{G}^r_{sp} + \sum _{j=1}^N \sum _{k=1}^N\; 6\;{H}^r_{spk} {u}_{k} \end{aligned}$$

and in compact form:

$$\begin{aligned} \dfrac{\partial ^2 \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}\partial \mathbf {u}}= 2\mathbf {G}+6\mathbf {Hu}. \end{aligned}$$

(86)

Appendix E: Derivation of modal derivatives

In this appendix, we derive the Taylor expansion of the nonlinear eigenproblem defined by Eq. (17) and recalled here for the sake of completeness:

$$\begin{aligned} \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} - {\tilde{\omega }}_i^2(\mathbf {u}) \mathbf {M}\right) \tilde{\mathbf {\phi }}_i(\mathbf {u}) = \mathbf {0}, \end{aligned}$$

(87)

Assuming small motions in the vicinity of the position at rest, this nonlinear eigenproblem can be expanded up to second order as:

$$\begin{aligned} \begin{aligned}&\left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} - {\tilde{\omega }}_i^2(\mathbf {u}) \mathbf {M}\right) \tilde{\mathbf {\phi }}_i(\mathbf {u})\\&\quad =\left( \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} - {\tilde{\omega }}_i^2(\mathbf {u}) \mathbf {M}\right) \tilde{\mathbf {\phi }}_i(\mathbf {u})\right) \bigg |_{\mathbf {u}=\mathbf {0}}\\&\qquad + \sum ^n_{j=1}\left( \dfrac{\partial }{\partial R_j} \left( \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} - {\tilde{\omega }}_i^2(\mathbf {u}) \mathbf {M}\right) \tilde{\mathbf {\phi }}_i(\mathbf {u}) \right) \right) \bigg |_{\mathbf {u}=\mathbf {0}}R_j\\&\qquad +{\mathcal {O}}(|\mathbf {R}|^2), \end{aligned} \end{aligned}$$

(88)

where the second term has been expanded along the coordinates $R_j$ used for the reduced basis, and the expansion has been written up to ${\mathcal {O}}(|\mathbf {R}|^2)$ terms, or equivalently ${\mathcal {O}}(|\mathbf {u}|^2)$ terms. As the constant term of the expansion, one retrieves the linear eigensystem of Eq. (5), since the Jacobian of the nonlinear force vector at $\mathbf {u}(\mathbf {R}=\mathbf {0})$ coincides with the linear stiffness matrix $\mathbf {K}$:

$$\begin{aligned} \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}}\bigg |_{\mathbf {u}=\mathbf {0}}=\mathbf {K}. \end{aligned}$$

(89)

Consequently, the first term of the expansion allows recovering the ith eigenvalue $\omega _i$ as well as the ith eigenvector $\mathbf {\phi }_i$.

To verify Eq. (17) up to first order, not only the constant term, but also all the linear terms in $R_j, \forall j=1,\ldots ,N$ must be zero. By expanding the jth term, one obtains the condition:

$$\begin{aligned}&\left( \left( \dfrac{\partial }{\partial R_j} \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}}\right) \right) \bigg |_{\mathbf {0}} - \dfrac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} \mathbf {M} \right) \mathbf {\phi }_i \nonumber \\&\quad + \left( \mathbf {K} - \omega _i^2 \mathbf {M} \right) \dfrac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} =\mathbf {0}. \end{aligned}$$

(90)

In this equation, the sought modal derivatives are the vector $\frac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j}\bigg |_{\mathbf {0}}$ and the other unknown of the system is the value $\frac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}}$. Moreover, by noticing that:

$$\begin{aligned}&\left( \dfrac{\partial }{\partial R_j} \left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}} \right) \right) \bigg |_{\mathbf {0}}= \left( \dfrac{\partial }{\partial \mathbf {u}}\left( \dfrac{\partial \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}}\right) \dfrac{\partial \mathbf {u}}{\partial R_j}\right) \bigg |_{\mathbf {0}}\nonumber \\&\quad = \left( \dfrac{\partial ^2 \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}\;\partial \mathbf {u}} \right) \bigg |_{\mathbf {0}} \mathbf {\phi }_j, \end{aligned}$$

(91)

one can write Eq. (90) as:

$$\begin{aligned}&\left( \dfrac{\partial ^2 \mathbf {F}(\mathbf {u})}{\partial \mathbf {u}\;\partial \mathbf {u}} \right) \bigg |_{\mathbf {0}} \mathbf {\phi }_j \mathbf {\phi }_i - \dfrac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} \mathbf {M} \mathbf {\phi }_i \nonumber \\&\quad + \left( \mathbf {K}-\omega _i^2 \mathbf {M}\right) \dfrac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} =\mathbf {0}. \end{aligned}$$

(92)

The first term can be further simplified by recalling the definition of the nonlinear force vector and the value of the second derivatives of it given in Eq. (86), leading to:

$$\begin{aligned}&2\mathbf {G} \mathbf {\phi }_j \mathbf {\phi }_i - \dfrac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} \mathbf {M} \mathbf {\phi }_i \nonumber \\&\quad + \left( \mathbf {K}-\omega _i^2 \mathbf {M}\right) \dfrac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} =\mathbf {0}. \end{aligned}$$

(93)

This is now an undetermined system of equation in the unknowns $\frac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}}$ and $\frac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}}$. To solve this system, the additional equation of mass normalisation must be introduced.

Following a similar approach, i.e. expanding in Taylor series the nonlinear mass normalisation equation, one obtains:

$$\begin{aligned}&\tilde{\mathbf {\phi }}_i(\mathbf {R})^T \mathbf {M} \tilde{\mathbf {\phi }}_i(\mathbf {R})-1= \left( \tilde{\mathbf {\phi }}_i(\mathbf {R})^T \mathbf {M} \tilde{\mathbf {\phi }}_i(\mathbf {R})\right) \bigg |_{\mathbf {0}} \nonumber \\&\quad -\,1 + \sum ^n_{i=1} \left( \dfrac{\partial }{\partial R_j} \left( \tilde{\mathbf {\phi }}_i(\mathbf {R})^T \mathbf {M} \tilde{\mathbf {\phi }}_i(\mathbf {R}) \right) \right) \bigg |_{\mathbf {0}} R_j\nonumber \\&\quad +\,{\mathcal {O}}(|\mathbf {R}|^2) . \end{aligned}$$

(94)

The constant term is verified by the linear eigenvectors $\mathbf {\phi }_i$, whereas the linear terms must be equal to zero. The linear term in $R_j$ becomes the required complement to Eq. (93). By expanding the derivatives in $R_j$, it reads:

$$\begin{aligned} \dfrac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}}^T \mathbf {M} \mathbf {\phi }_i+\mathbf {\phi }_i^T \mathbf {M} \dfrac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j} \bigg |_{\mathbf {0}} =0. \end{aligned}$$

(95)

In the usual case of symmetric mass matrix, the LHS only reduces to one of its terms as they are equal. In light of this, the system that permits to evaluate the modal derivatives reads:

$$\begin{aligned} \left[ \begin{array}{cc} \mathbf {K}-\omega _i^2 \mathbf {M} &{}\quad -\mathbf {M} \mathbf {\phi }_i \\ -\mathbf {\phi }_i^T \mathbf {M} &{} \quad 0 \end{array}\right] \left\{ \begin{array}{c} \frac{\partial \tilde{\mathbf {\phi }}_i(\mathbf {R})}{\partial R_j}|_{\mathbf {0}} \\ \frac{\partial {\tilde{\omega }}_i^2(\mathbf {R})}{\partial R_j}|_{\mathbf {0}} \end{array}\right\} = \left\{ \begin{array}{c} 2\mathbf {G} \mathbf {\phi }_j \mathbf {\phi }_i \\ 0 \end{array}\right\} .\nonumber \\ \end{aligned}$$

(96)

Appendix F: Derivation of the two-dof system from von Kármán beam equations

In this appendix, we give the detailed calculation for obtaining all the coefficients of the two-dof model used in Sect. 3.1, from the von Kármán beam model. As described in Sect. 3.1.1, we only refer to the coupling between the first flexural mode and the fourth axial mode of a clamped–clamped beam. The eigenfunctions and eigenvalues of these two modes can be found solving the linear eigenproblem for, respectively, flexural and longitudinal vibrations.

Recalling Eq. (55) and focusing only on their linear part read:

$$\begin{aligned}&{\ddot{W}}(X,T)+\dfrac{E I}{\delta S} W^{''''}(X,T)=0, \end{aligned}$$

(97a)

$$\begin{aligned}&{\ddot{U}}(X,T) - \dfrac{E}{\delta } U{''}(X,T)=0. \end{aligned}$$

(97b)

Equation (97a) can be solved assuming:

$$\begin{aligned} {W}(X,T)=Q(T)\varPhi (X), \end{aligned}$$

(98)

where $\varPhi (X)$ has to respect the clamped–clamped boundary conditions, namely $\varPhi (0)=\varPhi (L)=\varPhi ^{'}(0)=\varPhi ^{'}(L)=0$. The first three conditions are respected by the eigenfunction $\varPhi (X)$ of arbitrary amplitude A:

$$\begin{aligned} \varPhi (X)= & {} A (\cos {k X}-\cosh {k X}) (\sin {k L}-\sinh {k L}) \nonumber \\&- A (\sin {k X}-\sinh {k X}) (\cos {k L}-\cosh {k L}),\nonumber \\ \end{aligned}$$

(99)

whereas this condition gives rise to the wavelength equation:

$$\begin{aligned} \cos {k L}\cosh {k L}=1. \end{aligned}$$

(100)

The first value of $k>0$ that verifies the transcendental Eq. (100) is the dimensional wavelength $k_1$ of the first flexural mode. Its eigenfrequency is then obtained by solving Eq. (97a) and reads:

$$\begin{aligned} \omega _{1f}^2=\dfrac{E I}{\delta S}k_1^4. \end{aligned}$$

(101)

As for the fourth longitudinal mode, a similar approach is followed. Imposing the separation of variables on U:

$$\begin{aligned} U(X,T) = P(T)\varPsi (X), \end{aligned}$$

(102)

and imposing the clamped–clamped boundary conditions $\varPsi (0)=\varPsi (L)=0$, one can find the eigenfunction that has now a simpler form being Eq. (97b) a second-order differential equation in X. The first boundary condition is verified by:

$$\begin{aligned} \varPsi (X) = B \sin {\kappa X}, \end{aligned}$$

(103)

and the second one by the wavelength equation:

$$\begin{aligned} \sin {\kappa L}=0. \end{aligned}$$

(104)

The dimensional wavelength of the fourth axial mode is fourth value of $\kappa >0$ that respects Eq. (104) equal to $\kappa _4=4 \pi /L$. The eigenfrequency of the fourth axial mode is then obtained from Eq. (97b) and reads:

$$\begin{aligned} \omega _{4a}^2=\dfrac{E}{\delta }\dfrac{(4 \pi )^2}{L^2}. \end{aligned}$$

(105)

Before operating the reduction that will produce the two-dof system of ODEs, it is convenient to make Eq. (97) non-dimensional with the following identities:

$$\begin{aligned}&X = x L,\\&T = t T_0,\\&W = w h,\\&U = u h,\\&Q = q h,\\&P = p h, \end{aligned}$$

and by introducing the additional quantities:

$$\begin{aligned} c&= \sqrt{E/\delta },\\ \beta&= k_1 L,\\ \sigma&= h/L,\\ T_0&= 1/\omega _{1f} = \sqrt{12} L^2/(c h \beta ^2), \end{aligned}$$

where the rectangular section assumption is used in this equation to simplify $I/S=h^2/12$. It is now possible to rewrite Eq. (55) with respect to the new variables x, t, w, u as:

$$\begin{aligned}&\dfrac{h}{T_0^2} \left( w_{,tt} + \dfrac{1}{\beta ^4} w_{,xxxx} - \dfrac{12}{\beta ^4 \sigma } \left( u_{,x}w_{,x}\right) _{,x}\right. \nonumber \\&\quad \left. - \dfrac{6}{\beta ^4} \left( {w_{,x}}^3\right) _{,x} \right) =0, \end{aligned}$$

(106a)

$$\begin{aligned}&\dfrac{h}{T_0^2} \left( u_{,tt} - \frac{12}{\beta ^4 \sigma ^2}u_{,xx} - \frac{12}{\beta ^4 \sigma } w_{,x}w_{,xx} \right) =0, \end{aligned}$$

(106b)

which coincides with Eq. (56) multiplied by the nonzero factor $h/T_0^2$.

This system of equations can be reduced to a system of ODEs using the equations for w and u now in their non-dimensional form:

$$\begin{aligned}&w(x,t) = q_1(t)\varPhi _{1f}(x), \end{aligned}$$

(107)

$$\begin{aligned}&u(x,t) = p_4(t)\varPsi _{4a}(x), \end{aligned}$$

(108)

with:

$$\begin{aligned} \varPhi _{1f}(x)&= \alpha _1 (\cos {\beta x}-\cosh {\beta x}) (\sin {\beta }-\sinh {\beta }) \\&\quad - (\sin {\beta x}-\sinh {\beta x}) (\cos {\beta }-\cosh {\beta })\\ \varPsi _{4a}(x)&= \alpha _4\sin {4 \pi x}, \end{aligned}$$

by projecting Eq. (106a) on the shape function $\varPhi _{1f}$ and Eq. (106b) on the shape function $\varPsi _{4a}$.

The Galerkin projection of Eq. (106) leads to:

$$\begin{aligned}&q_{1,tt} + q_1 - \dfrac{2}{\sigma }G q_1\,p_4 +D q_1^3 =0, \end{aligned}$$

(109a)

$$\begin{aligned}&p_{4,tt} + \frac{12 (4\pi )^2}{\beta ^4 \sigma ^2} p_4 - \dfrac{1}{\sigma }C q_1^2 =0, \end{aligned}$$

(109b)

with the coefficients G, D, C being equal to:

$$\begin{aligned} G&=- \dfrac{6}{\beta ^4} \left( \dfrac{ \int _0^1\varPhi \left( \varPsi _{,x}\varPhi _{,x}\right) _{,x} \mathrm{d}x}{\int _0^1\varPhi ^2 \mathrm{d}x} \right) , \end{aligned}$$

(110)

$$\begin{aligned} C&=-\frac{12}{\beta ^4} \left( \dfrac{ \int _0^1\varPsi \left( \varPhi _{,x}\varPhi _{,xx}\right) \mathrm{d}x}{\int _0^1\varPsi ^2 \mathrm{d}x} \right) , \end{aligned}$$

(111)

$$\begin{aligned} D&=- \dfrac{6}{\beta ^4} \left( \dfrac{ \int _0^1\varPhi \left( {\varPhi _{,x}}^3\right) _{,x} \mathrm{d}x}{\int _0^1\varPhi ^2 \mathrm{d}x} \right) . \end{aligned}$$

(112)

If the arbitrary amplitudes $\alpha _1$ and $\alpha _4$ are chosen to have mass normalised eigenfunctions:

$$\begin{aligned}&\alpha _1: \int _0^1\varPhi ^2 \mathrm{d}x=1,\\&\alpha _4: \int _0^1\varPsi ^2 \mathrm{d}x=1, \end{aligned}$$

the quadratic coupling coefficients are symmetric $G=C=1.23$, the cubic coefficient is $D=2.67$ and Eq. 57 is finally retrieved.