Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques

Abstract

This paper aims at reviewing nonlinear methods for model order reduction in structures with geometric nonlinearity, with a special emphasis on the techniques based on invariant manifold theory. Nonlinear methods differ from linear-based techniques by their use of a nonlinear mapping instead of adding new vectors to enlarge the projection basis. Invariant manifolds have been first introduced in vibration theory within the context of nonlinear normal modes and have been initially computed from the modal basis, using either a graph representation or a normal form approach to compute mappings and reduced dynamics. These developments are first recalled following a historical perspective, where the main applications were first oriented toward structural models that can be expressed thanks to partial differential equations. They are then replaced in the more general context of the parametrisation of invariant manifold that allows unifying the approaches. Then, the specific case of structures discretised with the finite element method is addressed. Implicit condensation, giving rise to a projection onto a stress manifold, and modal derivatives, used in the framework of the quadratic manifold, are first reviewed. Finally, recent developments allowing direct computation of reduced-order models relying on invariant manifolds theory are detailed. Applicative examples are shown and the extension of the methods to deal with further complications are reviewed. Finally, open problems and future directions are highlighted.

Appendices

Symmetry of the quadratic and cubic tensors

This appendix is devoted to the demonstration of the symmetry properties of the nonlinear tensors of coefficients ${\mathbf {G}}$, ${\mathbf {H}}$, ${\mathbf {g}}$ and ${\mathbf {h}}$, that appear in the equations of motion written in terms of FE coordinates, Eq. 6), or with modal coordinates, Eq. (11).

The FE coordinates are first considered. Let us denote by ${\mathbf {k}}({\mathbf {X}})={\mathbf {K}}{\mathbf {X}}+{\mathbf {f}}_\text {nl}({\mathbf {X}})$ the internal force vector. Using indicial notations and Einstein summation convention, it can be written explicitly, for $i,j,l,s=1,\ldots N$:

$$\begin{aligned} k_s=K_{si}X_i+G^s_{ij}X_iX_j+H^s_{ijl}X_iX_jX_l, \end{aligned}$$

(51)

where $k_s$ is the s-th component of the internal force vector ${\mathbf {k}}$, $K_{si}$ are the components of the stiffness matrix ${\mathbf {K}}$; while $G^s_{ij}$, $H^s_{ijl}$ are the quadratic and cubic coefficients defined in Eq. (8). In a 3D finite element context, the physical displacement vector ${\mathbf {u}}({\mathbf {y}})$ is interpolated on a family of shape functions ${\mathbf {N}}_i$, such that $u_\alpha ({\mathbf {y}})=N_{\alpha i}({\mathbf {y}})X_i(t)$, $\alpha =1,2,3$. Using Voigt notations, the Green–Lagrange strain tensor, Eq. (5), and its variation can be written, for $\alpha =1,\ldots 6$ (see [67]):

$$\begin{aligned}&E_\alpha =B^{(1)}_{\alpha i}X_i+\frac{1}{2}B^{(2)}_{\alpha ij}X_iX_j,\quad \nonumber \\&\delta E_\alpha =B^{(1)}_{\alpha i}({\mathbf {y}})\delta X_i+B^{(2)}_{\alpha ij}X_i\delta X_j, \end{aligned}$$

(52)

with $B^{\alpha (1)}_{i}$ of size $6\times N$ and $B^{\alpha (2)}_{ij}$ of size $6\times N\times N$ are two discretised gradients operators, defined by:

$$\begin{aligned}&B^{(1)}_{\alpha i}=\begin{bmatrix} N_{1i,1} \\ N_{2i,2} \\ N_{3i,3} \\ N_{2i,3}+N_{3i,2} \\ N_{1i,3}+N_{3i,1} \\ N_{1i,2}+N_{2i,1} \end{bmatrix}\nonumber \\&\quad B^{(2)}_{\alpha ij}=\begin{bmatrix} N_{li,1 }N_{lj,1}\\ N_{li,2 }N_{lj,2}\\ N_{li,3 }N_{lj,3}\\ N_{li,2 }N_{lj,3}+N_{li,3 }N_{lj,2}\\ N_{li,3 }N_{lj,1}+N_{li,1 }N_{lj,3}\\ N_{li,2 }N_{lj,1}+N_{li,1 }N_{lj,2} \end{bmatrix} \end{aligned}$$

(53)

where each row of the matrices corresponds to the corresponding value of index $\alpha =1,\ldots 6$ and $N_{\alpha i,\beta }=\partial N_{\alpha i}/\partial y_\beta$ is a space derivative of the shape functions.

Then, the virtual work of the internal forces can be written:

$$\begin{aligned} \delta {\mathcal {W}}_\text {int}&=\int _\varOmega C_{\alpha \beta }E_\alpha \delta E_\beta \,{{\text {d}}}\varOmega = k_s\delta X_s \end{aligned}$$

(54)

$$\begin{aligned}&=\Bigg [\underbrace{\int _\varOmega C_{\alpha \beta }B_{\alpha i}^{(1)}B_{\beta s}^{(1)}\,{{\text {d}}}\varOmega }_{K_{si}} \,X_i \nonumber \\&\quad + \underbrace{\int _\varOmega C_{\alpha \beta }\left( B_{\alpha i}^{(1)}B_{\beta js}^{(2)} +\frac{1}{2}B_{\alpha s}^{(1)}B_{\beta ij}^{(2)}\right) \,{{\text {d}}}\varOmega }_{G^s_{ij}} \,X_i X_j \nonumber \\&\quad + \underbrace{\int _\varOmega \frac{1}{2}C_{\alpha \beta }B_{\alpha ij}^{(2)}B_{\beta ls}^{(2)}\,{{\text {d}}}\varOmega }_{H^s_{ijl}} \,X_i X_jX_l \Bigg ]\delta X_s\nonumber \\ \end{aligned}$$

(55)

where $C_{\alpha \beta }$ is the elasticity tensor in Voigt notations. The above equation defines the tensor components $K_{si}$, $G^s_{ij}$ and $H^s_{ijl}$ of the internal forces.

The elasticity tensor is symmetric ($C_{\alpha \beta }=C_{\beta \alpha }$) and Eq. (53) shows that $B^{(2)}_{\alpha ij}=B^{(2)}_{\alpha ji}$: on can invert the two latin subscripts. This leads to a symmetric stiffness matrix and allows any change of the order of the indices for the cubic coefficients ($4!=24$ possibilities if $i\ne j\ne l\ne s$):

$$\begin{aligned} H^s_{ijl} = H^l_{sij} = H^j_{lsi} = H^i_{jls} = H^s_{jil}=\ldots \end{aligned}$$

(56)

For the quadratic coefficients, no symmetry appears. However, the coefficients $G^s_{ij}$ and $G^s_{ji}$ refer to the same monomial $X_iX_j$ in Eq. (51). Due to the commutativity property of the usual product, one then understands that only the summation of these two quadratic coefficient $G^s_{ij}$ and $G^s_{ji}$ matters. This leads several authors (see, e.g. [69, 190, 287]) to adopt a so-called upper triangular form for those tensors for which Eq. (51) is rewritten as

$$\begin{aligned} k_s= & {} \sum _{i=1}^N K_{si}X_i + \sum _{i=1}^N \sum _{j=i}^N {\hat{G}}^s_{ij} X_i X_j \nonumber \\&+ \sum _{i=1}^N \sum _{j=i}^N \sum _{l=j}^N {\hat{H}}^s_{ijl} X_i X_j X_l. \end{aligned}$$

(57)

With this selection, an unequivocal representation of the monomials is given, and the coefficients are attributed such that only those with increasing indices ($l\ge j \ge i$) are nonzero, while the other ones ($l\le j \le i$) are set to zero. For the quadratic coefficients, this leads to, for all s:

$$\begin{aligned} {\hat{G}}^s_{ii}= & {} G^s_{ii}, \quad \forall i, \end{aligned}$$

(58)

$$\begin{aligned} {\hat{G}}^s_{ij}= & {} G^s_{ij}+G^s_{ji}=\int _\varOmega C_{\alpha \beta }\left( B_{\alpha i}^{(1)}B_{\beta js}^{(2)}\right. \nonumber \\&\left. +B_{\alpha j}^{(1)}B_{\beta is}^{(2)}+B_{\alpha s}^{(1)}B_{\beta ij}^{(2)}\right) \,{{\text {d}}}\varOmega ,\quad \forall j>i, \end{aligned}$$

(59)

$$\begin{aligned} {\hat{G}}^s_{ji}= & {} 0,\quad \forall j>i. \end{aligned}$$

(60)

In the above Eq. (59), any change of the order of the indices s, i, j is allowed, which leads to the following properties:

$$\begin{aligned} {\hat{G}}^i_{ij}= & {} 2{\hat{G}}^j_{ii},\quad {\hat{G}}^j_{ij}=2{\hat{G}}^i_{jj},\quad {\hat{G}}^l_{ij}={\hat{G}}^j_{il},\nonumber \\ {\hat{G}}^l_{ij}= & {} {\hat{G}}^i_{lj},\quad \forall i<j<l \end{aligned}$$

(61)

Analog expressions are obtained for the cubic coefficients ${\hat{H}}^s_{ijl}$ (see [69, 190]).

In the present article, we found convenient to use the standard form (51) instead of the upper triangular form (57) for all our demonstrations and we enforced the symmetry on the quadratic coefficients by redefining them as the symmetric part of their upper triangular counterparts: $G^s_{ij}=G^s_{ji}={\hat{G}}^s_{ij}/2$, $\forall j\ne i$. This leads to allows any change of the order of the indices also for the quadratic coefficients:

$$\begin{aligned} G^s_{ij}=G^s_{ji}=G^i_{js}=G^i_{sj}=G^j_{is}=G^j_{si}. \end{aligned}$$

(62)

All the above reasoning equally applies to the modal coefficients $g^s_{ij}$ and $h^s_{ijl}$ of Eq. (11), since, according to Eq. (12a), (12b), rewritten in full indicial form with Einstein notation, one has:

$$\begin{aligned}&g^p_{mn}=G^s_{ij}\phi _{sp}\phi _{im}\phi _{jn},\quad h^p_{mnq}=H^s_{ijl}\phi _{sp}\phi _{im}\phi _{jn}\phi _{lq},\nonumber \\&\qquad \forall p,m,n,q=1,\ldots N, \end{aligned}$$

(63)

where $\phi _{jn}$ refers to the j-th component of the n-th eigenvector. Consequently, the same symmetry as in Eqs. (62) and (56) applies for the modal coefficients:

$$\begin{aligned}&g^s_{ij}=g^s_{ji}=g^i_{js}=g^i_{sj}=g^j_{is}=g^j_{si}, \end{aligned}$$

(64)

$$\begin{aligned}&h^s_{ijl} = h^l_{sij} = h^j_{lsi} = h^i_{jls} = h^s_{jil}=\ldots \end{aligned}$$

(65)

Again, any permutation of the indices s, i, j, l is possible.

All the above symmetries are also a consequence of the existence of a potential energy [69, 190]:

$$\begin{aligned} {\mathcal {V}}=\frac{1}{2}\int _\varOmega C_{\alpha \beta }E_\alpha ({\mathbf {X}}) E_\beta ({\mathbf {X}})\,{{\text {d}}}\varOmega , \end{aligned}$$

(66)

which is a quartic polynomial in ${\mathbf {X}}$. Then, the i-th component of the internal force vector can be directly derived from the elastic energy as:

$$\begin{aligned} k_i({\mathbf {X}})= \dfrac{\partial {\mathcal {V}}}{\partial X_i}. \end{aligned}$$

(67)

Using Schwarz’s theorem, one can write:

$$\begin{aligned} \dfrac{\partial ^2 {\mathcal {V}}}{\partial x_i\partial x_j} = \dfrac{\partial ^2 {\mathcal {V}}}{\partial x_j\partial x_i}\quad \Rightarrow \quad \dfrac{\partial k_i}{\partial x_j} = \dfrac{\partial k_j}{\partial x_i}. \end{aligned}$$

(68)

Then, identifying the coefficients of identical monomials in the above last equations with the upper triangular form for $k_i$ leads to prove all the symmetry properties of Eq. (61) and their analogs for the cubic coefficients.

Classification of nonlinear terms

This section is devoted to give more details of the terminology used throughout the text to classify the different nonlinear coupling terms and monomials appearing in the dynamics. It is based on previous developments reported, e.g. in [280, 287], and all the textbooks dealing with normal form theory, where the reader can find more details. Systems with geometric nonlinearity are essentially driven by a large assembly of nonlinearly coupled oscillators, thus generating a very large number of coupling terms (of the order of $N^4$ terms). However, all the terms does not play the same role and it is important to identify the contributions of each monomials. For this discussion on the terminology, we use the equations of motion written in modal space, Eq. (11), and more specifically, $\forall \; p \, \in \, [1,\ldots ,N]$:

$$\begin{aligned}&\ddot{x}_p + \omega _p^2 x_p + \sum _{i=1}^{N} \sum _{j \ge i }^{N} g_{ij}^p x_i x_j \nonumber \\&\quad + \sum _{i=1}^{N} \sum _{j\ge i}^{N} \sum _{k \ge j}^{N} h_{ijk}^p x_ix_jx_k = 0. \end{aligned}$$

(69)

For the discussion herein, assume that m is the master mode such that most of the energy is contained within $x_m$. Invariant-breaking terms have already been commented in the main text, they are the monomials on p-th oscillator equation with the form $g^p_{mm} x_m^2$ and $h^p_{mmm} x_m^3$. As soon as $x_m \ne 0$, then all $x_p$ where this terms exist will also be excited and will thus have a nonzero amplitude. These terms break the invariance of the linear eigensubspaces and can be directly tracked in the equations defining the geometry of the invariant manifolds, as underlined in Sects. 3.2 and 3.3.

In order to go ahead in this classification, the link with internal resonance must be properly understood. Let us start by underlining that any nonlinear term can be interpreted as a forcing term for the corresponding oscillator equation. Continuing with the same example, the term $g^p_{mm} x_m^2$ is a forcing term on oscillator p. Interestingly, at the lowest order of approximation, $x_m \sim {{\,\mathrm{e}\,}}^{\pm i \omega _m t}$, such that $x_m^2$ will create a forcing with frequency components $2\omega _m$ and 0. From this, we can conclude that for oscillator p, if $\omega _p \simeq 2\omega _m$ then the forcing term $g^p_{mm} x_m^2$ will be a resonant forcing term, exciting component p in the vicinity of its eigenfrequency thus creating large amplitude response. One then call the monomial $g^p_{mm} x_m^2$ a resonant monomial. On the other hand, as long as $\omega _p \ne 2\omega _m$, then the forcing term is non-resonant, and the corresponding monomial is non-resonant.

This simple example generalizes thanks to normal form theory. The resonance relationship are then linked with internal resonance between the eigenfrequencies of the system, all of them being connected to a specific order of nonlinearity, conducting to so-called second-order internal resonance and third-order internal resonance. Due to the fact that for linear conservative system, the eigenspectrum is purely imaginary $\{ \pm i \omega _r \}$, some third-order relationships are always fulfilled, all those of the form:

$$\begin{aligned} \forall {r,p}=[1,\ldots ,N]: \quad +i\omega _r = +i\omega _p -i \omega _p + i\omega _r.\nonumber \\ \end{aligned}$$

(70)

These resonances are called trivial resonance and their associated monomials are called trivially resonant monomials. The main consequence in terms of normal form is that all these monomials cannot be cancelled from the normal form of the system. The normal form is not linear but stay nonlinear with only these trivially resonant monomials in case of no other internal resonances between the eigenfrequencies, following the general results from Poincaré and Poincaré-Dulac’s theorems.

From these developments, we can derive the following classification:

• trivially resonant monomials: for the m-th oscillator, all the terms $x_m^3$, $x_m x_p^2$, $\forall p=1,\ldots N$, are trivially resonant monomials, corresponding to the trivial resonance relationships (70). They cannot be cancelled from the normal form and stay in the resulting equations as ascertained in Eq. (29). Note in particular than none of these are invariant-breaking, recovering the fact that the dynamics is well expressed in an invariant-based span of the phase space for Eq. (29), which is not the case for the equations of motion in modal space, Eq. (69).

• resonant monomial: a resonant monomial is the nonlinear term connected to an internal resonance between the eigenfrequencies of the system. For any internal resonance, a few monomials exist which can be tracked from the simple interpretation of nonlinear term as forcing.

• non-resonant monomial: a nonlinear term that is not connected to an internal resonance.

When an internal resonance exist, the corresponding resonant monomials create what is generally called a strong, resonant coupling, and they convey the energy exchange between the coupled oscillators, leading to more complex form of the dynamics with bifurcation in a larger phase space. Otherwise, the coupling is termed weak or non-resonant.

Parametrisation of invariant manifold

In this appendix, more details on the parametrisation method for the computation of invariant manifold of vector fields in the vicinity of a fixed point, are given. The presentation follows strictly the one given by Haro et al., using their notations and developments, so that all the credit of the presentation reported here is given to [83]. Here, a simple summary (with slight simplifications in the presentation) is provided and the interested reader is referred to [83] for more details.

The unknowns ${\mathbf {W}}$ and ${\mathbf {f}}$ are searched as polynomial expansions of order k. They are computed order by order, under the form

$$\begin{aligned} {\mathbf {W}}({\mathbf {s}})&= {\mathbf {z}}_{\star } + {\mathbf {L}}{\mathbf {s}}+ \sum _{k \ge 2} {\mathbf {W}}_k ({\mathbf {s}}), \end{aligned}$$

(71a)

$$\begin{aligned} {\mathbf {f}}({\mathbf {s}})&= \varvec{\varLambda }_L {\mathbf {s}}+ \sum _{k \ge 2} {\mathbf {f}}_k ({\mathbf {s}}), \end{aligned}$$

(71b)

where ${\mathbf {L}}$ is the restriction of the matrix of eigenvectors to the master modes of interest, and $\varvec{\varLambda }_L$ the linear diagonalized part of the dynamics restricted to the master modes contained in ${\mathbf {L}}$. ${\mathbf {W}}_k ({\mathbf {s}})$ represents an n-dimensional vector of d-variate homogeneous polynomials of degree k, starting at second order (quadratic terms), while ${\mathbf {f}}_k$ is n-dimensional.

By replacing (71) into the invariance Eq. (17) and identifying terms with the same power k, on obtains the so-called order-k homological equation as

$$\begin{aligned} {{\text {D}}}{\mathbf {F}}(z_{\star }) {\mathbf {W}}_k ({\mathbf {s}}) - {\mathbf {L}}{\mathbf {f}}_k ({\mathbf {s}}) - {{\text {D}}}{\mathbf {W}}_k ({\mathbf {s}}) \varvec{\varLambda }_L {\mathbf {s}}= -{\mathbf {E}}_k ({\mathbf {s}}),\nonumber \\ \end{aligned}$$

(72)

where the order-k error term ${\mathbf {E}}_k ({\mathbf {s}})$ has been introduced as:

$$\begin{aligned} {\mathbf {E}}_k ({\mathbf {s}}) = \left[ {\mathbf {F}}({\mathbf {W}}_{<k} ({\mathbf {s}})) \right] _k - \left[ {{\text {D}}}{\mathbf {W}}_{<k} ({\mathbf {s}}) {\mathbf {f}}_{<k} ({\mathbf {s}}) \right] _k, \end{aligned}$$

(73)

and where the shortcut notation $\left[ \; .\; \right] _k$ refers to the selection of k-th order terms only, while ${\mathbf {W}}_{<k}$ refers to all orders strictly smaller than k.

More insight can be given by projecting onto the modal coordinates, then allowing separating contributions due to master and slave coordinates. Denoting as ${\mathbf {P}}$ the vector of eigenfunctions, one can introduce:

$$\begin{aligned} \varvec{\xi }_k ({\mathbf {s}}) = {\mathbf {P}}^{-1} {\mathbf {W}}_k ({\mathbf {s}}), \end{aligned}$$

(74)

the coefficients of the nonlinear mapping expressed in the modal basis, as well as

$$\begin{aligned} \varvec{\eta }_k ({\mathbf {s}}) = {\mathbf {P}}^{-1} {\mathbf {E}}_k ({\mathbf {s}}), \end{aligned}$$

(75)

the expression of the error-k vector in the modal basis.

The vector $\varvec{\xi }_k$ can be split as follows:

$$\begin{aligned} \varvec{\xi }_k ({\mathbf {s}}) = \left[ \begin{array}{c} \varvec{\xi }_k^L ({\mathbf {s}}) \\ \varvec{\xi }_k^N ({\mathbf {s}}) \end{array} \right] , \end{aligned}$$

(76)

where the first d lines $\varvec{\xi }_k^L$ is the tangent part, related to the original linear matrix ${\mathbf {L}}$ containing the master mode coordinates, and the last $n-d$ lines, $\varvec{\xi }_k^N$, refers to the normal part (slave coordinates). The normal part of Eq. (72) is now called the normal co-homological equation and reads

$$\begin{aligned} \varvec{\varLambda }_N \varvec{\xi }_k^N ({\mathbf {s}}) - D \varvec{\xi }_k^N ({\mathbf {s}}) \varvec{\varLambda }_L {\mathbf {s}}= \varvec{\eta }_k^N ({\mathbf {s}}), \end{aligned}$$

(77)

where the second member vector $\varvec{\eta }$ has been split following the same notation. One must first solve this equation as only depending on one unknown $\varvec{\xi }_k^N$. The remaining part is called the tangent co-homological equation and reads

$$\begin{aligned} \varvec{\varLambda }_L \varvec{\xi }_k^L ({\mathbf {s}}) - D \varvec{\xi }_k^L ({\mathbf {s}}) \varvec{\varLambda }_L {\mathbf {s}}- {\mathbf {f}}_k ({\mathbf {s}}) = \varvec{\eta }_k^L ({\mathbf {s}}), \end{aligned}$$

(78)

In order to fully analyze the co-homological equations and express their solutions, let us first denote as $\varvec{\xi }_k ({\mathbf {s}}) = [\xi _k^1 ({\mathbf {s}}), \ldots , \xi _k^n ({\mathbf {s}})]^t$ the n components of the vector $\varvec{\xi }_k ({\mathbf {s}})$, which all are homogeneous polynomial of degree k. The same notation is used for the second vector of unknowns, ${\mathbf {f}}_k ({\mathbf {s}}) = [f_k^1 ({\mathbf {s}}), \ldots , f_k^d ({\mathbf {s}})]^t$, which is d-dimensional and also composed of homogeneous polynomial of degree k. The known vector $\varvec{\eta }_k ({\mathbf {s}}) = [\eta _k^1 ({\mathbf {s}}), \ldots , \eta _k^n ({\mathbf {s}})]^t$ is developed as well following the same indicial notation. The normal part of the cohomological Eq. (77) can now simply be rewritten component by component, for $i=d+1, \ldots , n$:

$$\begin{aligned} \lambda _i \xi _k^i ({\mathbf {s}}) - D\xi _k^i ({\mathbf {s}}) \varvec{\varLambda }_L {\mathbf {s}}= \eta _k^i ({\mathbf {s}}). \end{aligned}$$

(79)

Let us denote as $\xi _{{\mathbf {m}}}^i$ the coefficient of the monomial term associated with the i-th line, $i=d+1, \ldots , n$, and to the vector of integers ${\mathbf {m}}$ such that ${\mathbf {m}}= [m_1, \ldots , m_d]^t$, where all $m_j$ are integers and $|{\mathbf {m}}| = m_1 + m_2 + \cdots + m_d = k$ is the order k of the polynomials considered. Saying things differently, $\xi _{{\mathbf {m}}}^i$ is the coefficient of the monomial term $s_1^{m_1}s_2^{m_2}\ldots s_d^{m_d}$, of order k, and each $\varvec{\xi }_k^i ({\mathbf {s}})$ is composed of the summations of all possible combinations of these order-k monomial terms. Since Eq. (77) is diagonal with respect to $\xi _{{\mathbf {m}}}^i$, one can write, for $i=d+1, \ldots , n$ and for $|{\mathbf {m}}|=k$:

$$\begin{aligned} \left( \lambda _i - {\mathbf {m}}\varvec{\lambda }_L \right) \xi _{{\mathbf {m}}}^i = \eta _{{\mathbf {m}}}^i, \end{aligned}$$

(80)

with the shortcut notation ${\mathbf {m}}\varvec{\lambda }_L = m_1 \lambda _1 + \cdots + m_d \lambda _d$.

A cross-resonance occurs if there exist pairs $({\mathbf {m}},i)$ such that $\lambda _i = {\mathbf {m}}\varvec{\lambda }_L$. If the system has no cross-resonance, then an explicit solution for the unknown coefficient $\xi _{{\mathbf {m}}}^i$ is found as

$$\begin{aligned} \xi _{{\mathbf {m}}}^i = \frac{\eta _{{\mathbf {m}}}^i}{\lambda _i - {\mathbf {m}}\varvec{\lambda }_L}. \end{aligned}$$

(81)

If a cross-resonance exist, then a strong coupling exist between one slave coordinate and the set of master coordinates, and there is an obstruction in solving the normal cohomological equation. This means that a strong nonlinear coupling exist between one master and one slave coordinate such that the initial choice is not good, and the remedy consists in enlarging the number of master coordinates by considering the slave resonant modes as masters.

Following the same notations, for the tangent co-homological equation, one arrives at, for all $i=1, \ldots , d$:

$$\begin{aligned} \lambda _i \xi _k^i ({\mathbf {s}}) - D\xi _k^i ({\mathbf {s}}) \varvec{\varLambda }_L {\mathbf {s}}- f_k^i ({\mathbf {s}})= {\tilde{\eta }}_k^i ({\mathbf {s}}), \end{aligned}$$

(82)

with the shortcut notation ${\tilde{\varvec{\eta }}}_k^L ({\mathbf {s}}) = \varvec{\eta }_k^L ({\mathbf {s}})$. Again, this last equation can be explicited in terms of the unknowns which are each of the coefficients of the monomial terms. Using the same notation, one arrives at the following, for all $i=1, \ldots , d$

$$\begin{aligned} \left( \lambda _i - {\mathbf {m}}\varvec{\lambda }_L \right) \xi _{{\mathbf {m}}}^i - f_{{\mathbf {m}}}^i= {\tilde{\eta }}_{{\mathbf {m}}}^i. \end{aligned}$$

(83)

The pairs $({\mathbf {m}},i) \in {\mathbb {N}}^d \times {1,\ldots ,d}$, with $|{\mathbf {m}}| \ge 2$ such that $\lambda _i = {\mathbf {m}}\varvec{\lambda }_L$ create an internal resonance, referring to a nonlinear resonance relationship between the eigenvalues of the master coordinates. When such an internal resonance exist, then there is an obstruction to the linearisation of the reduced dynamics.

Since in (83) two unknowns are present, namely the coefficient of the monomials of both the nonlinear change of coordinate $\xi _{{\mathbf {m}}}^i$ and the reduced dynamics $f_{{\mathbf {m}}}^i$, the solution to this equation is not unique and can also be given even if there exist some resonances. This explains why there exist many different ways of solving the problem, thus leading to the different styles of parametrisation. The two main style of solutions are given as the graph style and the normal form style.

The graph style is simply obtained by stating that from order $k=2$, all the corrections contained in the master coordinates $\varvec{\xi }_k^L ({\mathbf {s}})$, are vanishing: $\varvec{\xi }_k^L ({\mathbf {s}})=0$. This means that the master coordinates are only linearly related to the original ones. With this assumption, one can then replace in Eq. (83) to arrive at the terms allowing one to write the reduced dynamics as, for $i=1,\ldots d$, $|m|=k$:

$$\begin{aligned} f_{{\mathbf {m}}}^i= -{\tilde{\eta }}_{{\mathbf {m}}}^i, \quad \xi _{{\mathbf {m}}}^i = 0. \end{aligned}$$

(84)

With this choice, one recovers the classical technique promoted from centre manifold theorem giving as initial guess a functional relationship (graph) between slave and master coordinates.

In the normal form style, the idea is to simplify as much as possible the reduced-order dynamics, by keeping only the resonant monomials, and discarding all other non-essential terms for the dynamical analysis. This leads to a more complex calculation, and a full nonlinear mapping between original coordinates and reduced ones, and at the end one arrives at a normal form for the reduced vector fields ${\mathbf {f}}$. The drawback is that calculations are a bit more involved (which is particularly true when there are numerous internal resonances to handle). The advantage is that the parametrisation is able to go over the foldings of the manifold.

To this end, one solves Eq. (83) following the rules (depending on the presence of internal resonance or not):

$$\begin{aligned} \text{ If } \; \lambda _i \ne {\mathbf {m}}\varvec{\lambda }_L, \quad f_{{\mathbf {m}}}^i&=0, \quad \xi _{{\mathbf {m}}}^i =\frac{{\tilde{\eta }}_{{\mathbf {m}}}^i}{\lambda _i - {\mathbf {m}}\varvec{\lambda }_L}, \end{aligned}$$

(85a)

$$\begin{aligned} \text{ If } \; \lambda _i = {\mathbf {m}}\varvec{\lambda }_L, \quad f_{{\mathbf {m}}}^i&= - {\tilde{\eta }}_{{\mathbf {m}}}^i, \quad \xi _{{\mathbf {m}}}^i = 0. \end{aligned}$$

(85b)

The formulas given for this invariant manifold computation with normal form style can be extended to the case $d=n$, and one then strictly recovers the usual full normal form of the original system [83].

Comparison of reduced dynamics

The aim of this appendix is to demonstrate the equivalence between the reduced dynamics given by the graph style as derived using the Shaw and Pierre approach, Eq. (25), to the one obtained thanks to real normal form approach, Eq. (32). In both case, reduction to a single master mode is used. Theoretically speaking, the two methods compute the same manifold and should thus provide the same dynamics. However, their formulations are differing since they are not expressed with the same variables. Let us start from the dynamics obtained using the graph style, rewritten here for the ease of reading:

$$\begin{aligned}&\ddot{x}_m + \omega _m^2 x_m + g^m_{mm} x_m^2 + x_m \left( \underset{s\ne m}{\sum _{s=1}^N} 2\,g^m_{ms} g^s_{mm}\right. \nonumber \\&\quad \left. \left[ \frac{2\omega _m^2 - \omega _s^2}{\omega _s^2 (\omega _s^2 - 4 \omega _m^2)}x_m^2 + \frac{2}{\omega _s^2 (\omega _s^2 - 4 \omega _m^2)}y_m^2 \right] \right) \nonumber \\&\qquad + h^m_{mmm} x_m^3 = 0. \end{aligned}$$

(86)

Using standard symmetry relationships on the quadratic coefficients, namely $g^m_{sm} = g^m_{ms} = g^s_{mm}$, the equation can be rewritten as

$$\begin{aligned}&\ddot{x}_m + \omega _m^2 x_m + g^m_{mm} x_m^2 + h^m_{mmm} x_m^3 + \left( \underset{s\ne m}{\sum _{s=1}^N} (g^s_{mm})^2 \frac{2}{\omega _s^2} \right. \nonumber \\&\quad \left. \left[ \frac{2\omega _m^2 - \omega _s^2}{ (\omega _s^2 - 4 \omega _m^2)}x_m^3 + \frac{2}{ (\omega _s^2 - 4 \omega _m^2)}x_m y_m^2 \right] \right) = 0.\nonumber \\ \end{aligned}$$

(87)

On the other hand, the reduced dynamics given by normal form writes:

$$\begin{aligned}&{\ddot{R}}_m + \omega _m^2 R_m + h_{mmm}^m R_m^3 + \sum _{s=1}^n ({g}^s_{mm})^2 \frac{2}{\omega ^2_{s}} \nonumber \\&\quad \left( \frac{2 \omega ^2_m - \omega ^2_s}{\omega ^2_s-4 \omega ^2_m}\;R_m^3 + \frac{2}{\omega ^2_s-4 \omega ^2_m}\;{\dot{R}}_m^2 R_m \right) \;=\;0.\nonumber \\ \end{aligned}$$

(88)

Comparing term by term the two equations, one can observe two main differences: the presence of the quadratic term $g^m_{mm} x_m^2$, and the summation which excludes the term $s=m$ in the first equation. The nonlinear relationship between the modal and the normal variables, in this case of a single master coordinate, reads:

$$\begin{aligned} x_m = R_m - \frac{1}{3\omega _m^2}g^m_{mm} R_m^2 - \frac{2}{3\omega _m^2}g^m_{mm} S_m^2 + {{\mathcal {O}}}(R_m^4,S_m^4),\nonumber \\ \end{aligned}$$

(89)

where the shortcut notation $S_m = {\dot{R}}_m$ is used. This equation is simply obtained from Eq. (28), assuming a single-mode motion, and replacing $a^m_{mm}$, $b^m_{mm}$ and $\gamma ^m_{mm}$ with their exact analytical values given in [287]. Note that this expansion is valid up to fourth-order, since the cubic coefficients are vanishing: $r^m_{mmm}=u^m_{mmm}=0$. Consequently, no cubic terms are present in (89). Replacing (89) in (87), and denoting as $T=\ddot{x}_m + \omega _m^2 x_m + g^m_{mm} x_m^2$ the term that will produce extra quadratic and cubic terms, one arrives at:

$$\begin{aligned} T&= {\ddot{R}}_m + \omega _m^2 R_m + \frac{2g^m_{mm}}{3} R_m^2 - \frac{2g^m_{mm}}{3\omega _m^2} S_m^2 \nonumber \\&\quad - \frac{2g^m_{mm}}{3\omega _m^2} \left( {\dot{R}}_m^2 + R_m {\ddot{R}}_m \right) - \frac{4g^m_{mm}}{3\omega _m^4} \left( {\dot{S}}_m^2 + S_m {\ddot{S}}_m \right) \nonumber \\&\quad - \frac{2(g^m_{mm})^2}{3\omega _m^2} R_m^3 -\frac{4(g^m_{mm})^2}{3\omega _m^4} R_m {\dot{R}}_m^2, \end{aligned}$$

(90)

where the first line gathers linear and quadratic terms and the second the cubic terms. One can first observe that the linear terms will produce the same as those in (88), a general property of identity-tangent nonlinear mapping. In order to obtain the equivalence between the two formulations, the goal is to show that the quadratic terms are vanishing. This is easily achieved by assuming that, at lower order, ${\dot{R}}_m = i\omega _m R_m$, ${\ddot{R}}_m = -\omega _m^2 R_m$ and so on. Replacing all the combinations, one obtains that the quadratic terms exactly cancels. Also the cubic terms appearing in (90) are exactly the one obtains from the summation in Eq. (88) for $s=m$. Consequently, the two equations are strictly equivalent up to the third order.

Implicit static condensation

This appendix aims at giving a few more details on the static condensation method, underlining the link that the ICE method shares with explicit condensation and the role of invariant-breaking terms to produce the curvatures of the stress manifold. This section use explanations reported in [255], and the interested reader is referred to this paper for more details.

Let us first underline the role of invariant-breaking terms in the construction of the stress manifold. For the sake of simplicity, let us assume modal expansion for the equations of motion, Eq. (13), and that only $x_m$ is selected as master mode. The method consists in applying a static body force of the form ${\mathbf {f}}_\text {e} = \beta _m {\mathbf {M}}\varvec{\phi }_m$ for several values of $\beta _m\in {\mathbb {R}}$, to compute with the FE code the corresponding displacement ${\mathbf {X}}(\beta _m)$ and the modal coordinates $x_i(\beta _m)$ by modal expansion. It results in solving the following system:

$$\begin{aligned}&\omega _m^2 x_m + \sum _{i=1}^{N} \sum _{j=1}^{N} g_{ij}^m x_i x_j \nonumber \\&\quad + \sum _{i=1}^{N} \sum _{j=1}^{N} \sum _{k=1}^{N} h_{ijk}^m x_ix_jx_k = \beta _m, \end{aligned}$$

(91a)

$$\begin{aligned}&\forall s \ne m,\quad \omega _s^2 x_s + \sum _{i=1}^{N} \sum _{j=1}^{N} g_{ij}^s x_i x_j \nonumber \\&\quad + \sum _{i=1}^{N} \sum _{j=1}^{N} \sum _{k=1}^{N} h_{ijk}^s x_ix_jx_k = 0, \end{aligned}$$

(91b)

in which the forcing is aligned with the m-th eigenvector, resulting in a zero forcing of the other oscillators, for $s\ne m$. Because of this last property and the implicit function theorem, Eq. (91b) leads to the existence of a static nonlinear relationship between the slave coordinates $x_s$ and the master one $x_m$, expressed formally as:

$$\begin{aligned} x_s = c_s (x_m),\quad \forall s\ne m. \end{aligned}$$

(92)

Replacing in Eq. (91a), the reduced-order dynamics simply reads:

$$\begin{aligned} \ddot{x}_m + \underbrace{\omega _m^2 x_m + \sum _{i=1}^{N} \sum _{j=1}^{N} g_{ij}^m c_i (x_m) c_j(x_m) + \sum _{i=1}^{N} \sum _{j=1}^{N} \sum _{k=1}^{N} h_{ijk}^m c_i (x_m) c_j(x_m) c_k(x_m) }_{\beta _m(x_m)}= 0, \end{aligned}$$

(93)

in which $\beta _m(x_m)$ can be identified as a polynomial in $x_m$ to obtain the reduced order dynamics at any order. From the above developments, it is clear that the dynamics of Eq. (93) is equivalent to the one of the full model with all the slave modal coordinates $x_s$ statically condensed into the master dynamics.

The general explicit solution to Eq. (92) with closed formulation is generally out of reach such that the $c_s$ functions are known implicitly. Nevertheless, they can be searched for as polynomial expansions and the first terms can be found. In particular, the quadratic term of the development is easy to find and is sufficient to derive the third-order dynamics, it reads [255]:

$$\begin{aligned} c_s (x_m) = -\frac{g^s_{mm}}{\omega _s^2} x_m^2 + {\mathcal {O}}(x_m^3), \end{aligned}$$

(94)

which, substituted into Eq. (93) with a subsequent truncation up to third order leads to Eq. (37).

Moreover, Eq. (91b) shows that in the nonlinear terms, the important invariant-breaking terms $g^s_{mm} x_m^2$ and $h^s_{mmm} x_m^3$ have a large magnitude, and these terms create a nonzero static response for $x_s$. These couplings make the manifold depart from the linear eigensubspace to create the stress manifold.