Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements

Abstract

Non-intrusive methods have been used since two decades to derive reduced-order models for geometrically nonlinear structures, with a particular emphasis on the so-called STiffness Evaluation Procedure (STEP), relying on the static application of prescribed displacements in a finite-element context. We show that a particularly slow convergence of the modal expansion is observed when applying the method with 3D elements, because of nonlinear couplings occurring with very high frequency modes involving 3D thickness deformations. Focusing on the case of flat structures, we first show by computing all the modes of the structure that a converged solution can be exhibited by using either static condensation or normal form theory. We then show that static modal derivatives provide the same solution with fewer calculations. Finally, we propose a modified STEP, where the prescribed displacements are imposed solely on specific degrees of freedom of the structure, and show that this adjustment also provides efficiently a converged solution.

1 Introduction

Geometrically nonlinear effects appear generally in thin structures such as beams, plates and shells, when the amplitude of the vibration is of the order of the thickness [26, 37]. The von Kármán family of models for beams, plates and shells allows one to derive explicit partial differential equations (PDE) [1, 8, 18, 37, 39], showing clearly that a coupling between bending and longitudinal motions causes a nonlinear restoring force of polynomial type in the equations of motion. This geometric nonlinearity is then at the root of complex behaviours, that also need dedicated computational strategies in order to derive quantitative predictions. On the phenomenological point of view, structural nonlinearities give rise to numerous nonlinear phenomena that have been analysed in a number of studies: frequency dependence on amplitude [14, 20, 25], hardening/softening behaviour [45, 46], hysteresis and jump phenomena [24, 36], mode coupling through internal resonances [9, 24, 39], bifurcations and loss of stability [10, 44], chaotic and turbulent vibrations [3, 5]. On the computational point of view, nonlinear couplings break the invariance property of the linear eigenmodes. Consequently, deriving reduced-order models (ROM) is no longer a straightforward problem and care has to be taken in order to find out a ROM that is capable of describing the dynamics of the whole system without losing accuracy.

When the structure under study is discretized with the finite element (FE) method, the problem of deriving accurate ROM is more stringent since the user cannot rely on a PDE in order to unfold an ad hoc mathematical method for building the ROM. Moreover, because of the intrinsic nature of the geometrical nonlinearities, all degrees of freedom of the FE model are nonlinearly coupled and substructuring techniques are not suitable, on the contrary to localized nonlinearities occurring frequently in contact and friction problems [17, 52]. Consequently, for geometric nonlinearity, the computation of the nonlinear coupling coefficients is as important as finding out a correct reduced basis, and sometimes the two problems are interwoven. Moreover, a number of recent studies highlighted the importance of using so-called indirect or non-intrusive methods, where the idea is to use the standard operations that any FE code is used to perform in order to build the ROM, hence avoiding the need to enter in the code and write new lines computing the needed quantities. On the other hand, direct methods also exist, for which there is a need to implement the computations at the level of the element [32, 47]. The STEP (STiffness Evaluation Procedure) is a non-intrusive method and has been first introduced by Muravyov and Rizzi [23]. In its first version as described in [23], it allows computation of the nonlinear coupling coefficients of the discretized problem in the modal basis from a series of static computations with prescribed modal displacements. It has then been used in a number of contexts [8, 19, 21, 22, 28], and is generally connected to the modal basis. However one has to understand that per se, STEP is just an evaluation technique, a computational non-intrusive method, that can be used with other inputs than those from the modal basis.

The STEP, although being largely applied in numerous cases, is known to suffer from a number of problems, making it not as so simple as its formulation could let one think. A first one lies in the amplitude of the prescribed displacement one has to impose in order to excite sufficiently the nonlinearity. As shown in [8], there is a clear range of amplitude for which the method works properly, between a minimal value where nonlinearity is not sufficiently excited and a maximum value for which other nonlinear effects are appearing. Another problem is connected to the use of the modal basis with a STEP computation, and must be interpreted as a drawback of using the modal basis for nonlinear computation, but is not directly linked with the STEP calculation. The problem is that of the slow convergence linked to the loss of invariance of eigenspaces, and the numerous couplings between low frequency bending modes and high frequency longitudinal modes, as underlined in a number of papers, see e.g. [8, 22] and references therein. Consequently, numerous methods have been proposed in order to overcome this limitation: dual modes [16], POD modes [29], discrete empirical interpolation [41], modal derivatives [12, 51], quadratic manifold [13, 30], just to name a few.

Incidentally, the majority of papers where the STEP has been applied to thin structures, report results obtained with structures discretized with beam, plate or shell finite elements, as being mostly interested in slender structures. A few examples with block elements can be found in the literature, see e.g. [27, 50], where dual modes have been used in order to achieve convergence. Using 3D finite elements can be interesting in practice, sometimes mandatory. In engineering applications, structures are often defined with a 3D geometrical model, for which a 3D FE discretization is straightforward. In other cases, some particular physical effects (piezoelectricity for instance) are sometimes implemented in existing FE codes only with 3D elements. Preliminary studies of the authors reveal a number of difficulties when blindly applying the STEP with modal basis to 3D finite elements, with more stringent inaccuracies and problems than those encountered with plate or shell elements. In particular, an unexpected slow convergence was observed when using the modal basis, and very high frequency modes appear to be involved.

The objective of this paper is to diagnose properly the issues one can encounter when applying the STEP with a modal basis to a structure discretized with 3D finite elements and present methods to overcome the problems. In the course of the paper, we will show that the problem is intrinsically related to the use of the modal basis, and that the STEP can be used with other inputs than the modal basis, in order to get better results. We restrict our attention to thin structures that are symmetric in the thickness direction, such as straight beams or plates, for which there is a bending / longitudinal uncoupling at the linear level that greatly simplifies the understanding of the phenomena and enables to obtain reference results.

The paper is organized as follows. Section 2 is dedicated to the framework of the study, the equations of motion and a brief recall of the STEP. Section 3 addresses on test examples the issues of the STEP applied to 3D FE. It is shown that a nonlinear coupling of bending modes with very high frequency modes involving deformations in the thickness of the structures occurs, and thus called thickness modes. Those modes are the result of 3D deformation effects that are not present when using beam or plate models. This unexpected coupling is different from the traditional bending-longitudinal coupling and the numerical examples show that they are of prime importance to achieve a converged solution. If one can compute all the coupled modes, two strategies are given to overcome the large dimension of the reduced basis: static condensation and the reduction to a single nonlinear normal modes. Both method shows that when a single master mode drives the dynamics, the ROM can still be composed of a single nonlinear oscillator. However, finding out all the coupled high-frequency modes is generally out of reach for complex structures. In Sect. 4, we then show that using a static condensation of a single modal derivative allows to retrieve the same converged result, in a more direct and efficient way. In addition, a modified version of the STEP is proposed, to directly embed the coupling with thickness modes. It consists in applying the prescribed displacement only on selected degrees of freedom and let the other free. In Sect. 5, the physical mechanisms of those nonlinear couplings are explained and Sect. 6 presents numerical results to validate the proposed numerical methods, able to overcome the convergence issues of the classical STEP.

2 Modelling

2.1 Reduced order model and STEP

An elastic mechanical structure discretized by the finite element method and having geometrical nonlinearities is considered. The time-dependent displacement vector ${{{\varvec{x}}}}$ gathers all the degrees of freedom of the model (displacements/rotations at each nodes) and is N-dimensional. The equation can be written:

$$\begin{aligned} {{{\varvec{M}}}}\ddot{{{{\varvec{x}}}}} + {{{\varvec{C}}}}\dot{{{{\varvec{x}}}}} + {{{\varvec{f}}}}({{{\varvec{x}}}}) = {{{\varvec{f}}}}_e, \end{aligned}$$

(1)

where ${{{\varvec{M}}}}$ and ${{{\varvec{C}}}}$ are the $N\times N$ dimensional mass and damping matrices, ${{{\varvec{f}}}}({{{\varvec{x}}}})$ is the internal force vector, ${{{\varvec{f}}}}_e(t)$ is the external force vector and the overdot represents the classical differentiation with respect to time t: $\dot{\bullet }={\text {d}}\bullet /{\text {d}}t$. Note that since we are more interested in the computation of the nonlinear restoring force, a linear viscous damping model has been used. In the present case of geometrically nonlinear structures, the internal force vector encompasses only polynomial terms up to order three and involves the displacement vector ${{{\varvec{x}}}}$ only [11, 19, 22, 47]. For the present study, we restrict our attention to this case, but extensions to more complex cases involving for examples velocity terms can also be handled.

It is convenient to split the internal force vector into a linear part and a purely nonlinear part. Assuming that the equilibrium point, the structure at rest, is given by ${{{\varvec{x}}}} = {{{\varvec{0}}}}$, the tangent stiffness matrix ${{{\varvec{K}}}}$ classically writes:

$$\begin{aligned} {{\varvec{K}}}=\left. \frac{\partial {{\varvec{f}}}}{\partial {{\varvec{x}}}}\right| _{{{\varvec{x}}}={{\varvec{0}}}}, \end{aligned}$$

(2)

so that the nonlinear internal force vector is defined as

$$\begin{aligned} {{\varvec{f}}}_\text {nl}({{\varvec{x}}}) = {{\varvec{f}}}({{\varvec{x}}}) - {{\varvec{K}}}{{\varvec{x}}}, \end{aligned}$$

(3)

and the equations of motion reads:

$$\begin{aligned} {{\varvec{M}}}\ddot{{{\varvec{x}}}} + {{{\varvec{C}}}}\dot{{{{\varvec{x}}}}} +{{\varvec{K}}}{{\varvec{x}}} + {{\varvec{f}}}_\text {nl}({{\varvec{x}}}) = {{\varvec{f}}}_\text {e} , \end{aligned}$$

(4)

where the geometrically nonlinear part of the problem is concentrated in the nonlinear internal force vector ${{\varvec{f}}}_\text {nl}({{\varvec{x}}})$.

The modal basis can be used as a first step in order to make the linear terms diagonal. The eigenmodes are the family of couples eigenfrequency-eigenvectors $(\omega _k,{{{\varvec{\phi }}}}_k)$, $k=1,\ldots , N$, solutions of the undamped, free and linearized Eq. (4):

$$\begin{aligned} ({{{\varvec{K}}}}-\omega _k^2{{{\varvec{M}}}}){{{\varvec{\phi }}}}_k={{{\varvec{0}}}}. \end{aligned}$$

(5)

Assuming the modal expansion for the displacement vector:

$$\begin{aligned} {{{\varvec{x}}}}(t)=\sum _{k=1}^{N} {{{\varvec{\phi }}}}_k X_k(t), \end{aligned}$$

(6)

where $X_k(t)$ is the modal amplitude, and using Galerkin projection allows one to rewrite the equations of motion in the modal space, for all $k=1,\ldots , N$, as:

$$\begin{aligned}&\ddot{X}_k+ 2\zeta _k\omega _k\dot{X}_k+ \omega _k^2X_k+ \sum _{i=1}^{N}\sum _{j=i}^{N} \alpha _{ij}^k X_iX_j\nonumber \\&\quad + \sum _{i=1}^{N}\sum _{j=i}^{N}\sum _{l=j}^{N} \beta _{ijl}^k X_i X_j X_l = Q_k, \end{aligned}$$

(7)

with

$$\begin{aligned} Q_k={{{{\varvec{\phi }}}}_k^{{\text {T}}}}{{{\varvec{f}}}}_e/m_k,\qquad m_k={{{{\varvec{\phi }}}}_k^{{\text {T}}}}{{{\varvec{M}}}}{{{\varvec{\phi }}}}_k, \end{aligned}$$

(8)

where $m_k$ is the modal mass, and where a modal damping (of factor $\zeta _k$) has been assumed uncoupled (an assumption valid for small damping, even with non-proportional ${{{\varvec{C}}}}$ matrix [7]). The nonlinear part of this reduced order is written with quadratic and cubic polynomial terms with coefficients $\alpha _{ij}^k$ and $\beta _{ijl}^k$. It is an exact expansion in the case of 3D FE [47] or beam/plate/shell FE based on von Kármán strain/displacement law [19], whereas it is truncated in the case of geometrically exact theories [38].

In Eq. (7), the linear parameters $\omega _k$, ${{\varvec{\phi }}}_k$ and $Q_k$ are obtained by the modal analysis of Eq. (5), available in any finite element code. The main issue is thus to compute the nonlinear coupling coefficients $\alpha _{ij}^k$ and $\beta _{ijl}^k$. The STEP (STiffness Evaluation Procedure) when used with the modal basis as first introduced by Muravyov and Rizzi [23], is a non-intrusive (or indirect) method, allowing one to get these coupling coefficients from standard computations available in any FE code. It relies on imposing prescribed static displacements having the shapes of selected eigenmodes, with a given amplitude. From a clever choice of the modes and the amplitudes, a simple algebra allows one to retrieve all the coefficients from the internal force vector given by the FE code, the key idea being to impose plus/minus the displacement with selected combinations of modes. The reader can find detailed explanation on this calculation in a number of papers, including the original one [23], as well as the improvement proposed in [28] using the tangent stiffness matrix.

To illustrate the method, we only show the computation of coefficients $\alpha _{pp}^k$ and $\beta _{ppp}^k$; for the general case the interested reader is referred to [23, 28]. The following static displacements are prescribed to the structure:

$$\begin{aligned} {{\varvec{x}}}_p =\pm \lambda {{{\varvec{\phi }}}}_{p}\quad \Rightarrow \quad \left\{ \begin{array}{l} q_p=\lambda , \\ q_j=0\quad \forall j\ne p, \end{array}\right. \end{aligned}$$

(9)

where $\lambda$ refers to an amplitude coefficient of the eigenvector ${{{\varvec{\phi }}}}_{p}$ whose value has to be chosen so as to activate the geometrical nonlinearities. A value of h/20, where h is the thickness of the plate, is recommended in [8]. Since the modes are orthogonal, imposing a displacement along mode p is equivalent to consider that only the modal coordinate $q_p$ is not vanishing, as detailed in the second part of Eq. (9). Since ${{\varvec{x}}}_p$ is time independent, introducing Eq. (9) into Eqs. (4) and (7) leads to, for all $k=1,\ldots , N$:

$$\begin{aligned} \lambda ^2 \alpha _{pp}^k + \lambda ^3 \beta _{ppp}^k&= {{{{\varvec{\phi }}}}_k^{{\text {T}}}} {{\varvec{f}}}_\text {nl}(\lambda {{\varvec{\phi }}}_p)/m_k, \end{aligned}$$

(10a)

$$\begin{aligned} \lambda ^2 \alpha _{pp}^k - \lambda ^3 \beta _{ppp}^k&= {{{{\varvec{\phi }}}}_k^{{\text {T}}}} {{\varvec{f}}}_\text {nl}(-\lambda {{\varvec{\phi }}}_p)/m_k. \end{aligned}$$

(10b)

Hence the unknown quadratic and cubic coefficients are found easily as

$$\begin{aligned} \alpha _{pp}^k&= \frac{1}{2\lambda ^2}\frac{{{{\varvec{\phi }}}}_k}{m_k}\left( {{\varvec{f}}}_\text {nl}(\lambda {{\varvec{\phi }}}_p) + {{\varvec{f}}}_\text {nl}(-\lambda {{\varvec{\phi }}}_p) \right) , \end{aligned}$$

(11a)

$$\begin{aligned} \beta _{ppp}^k&= \frac{1}{2\lambda ^3}\frac{{{{\varvec{\phi }}}}_k}{m_k}\left( {{\varvec{f}}}_\text {nl}(\lambda {{\varvec{\phi }}}_p) - {{\varvec{f}}}_\text {nl}(-\lambda {{\varvec{\phi }}}_p) \right) . \end{aligned}$$

(11b)

Similar algebraic manipulations with more modes involved in the prescribed displacement then allows one to get the full family of quadratic and cubic coefficients. The non intrusive nature of the method appears clearly: in the FE software, one prescribes the displacement field ${{\varvec{x}}}_p$ and computes the corresponding external force vector ${{\varvec{f}}}_\text {e}={{\varvec{f}}}({{\varvec{x}}}_p)$ with successive linear and nonlinear static computations. Then, ${{\varvec{f}}}_\text {nl}({{\varvec{x}}}_p)$ is obtained with Eq. (3).

2.2 The case of flat structures

In this article, we restrict the analysis to flat structures, such as straight beams or plates with a boundary of arbitrary shape. The thickness of the plate can be non-constant, but the geometrical and material distribution must be symmetric with respect to the middle line/plane of the structure.

If a plate theory, with a Kirchhoff–Love kinematics for instance, is applied to this structure, the displacement field of any point of the structure is described by the displacement field of the middle plane. There is a membrane/bending decoupling in linear elasticity and two families of modes are obtained: bending modes, for which only the transverse component of the middle plane displacement are non-zero, and membrane modes, for which the transverse component of the middle plane displacement field is zero. Analogous properties are valid for a beam theory, with a longitudinal/transverse decoupling on the middle line.

In the present case of a 3D structure which has the shape of a beam / plate, it is also possible to split the eigenmodes into two families in the same manner. The first family includes the bending modes, analogous to the ones obtained in the plate theory. Their frequencies are in the lower part of the spectrum, while their deformed shapes are dominated by transverse displacements. They, up to the accuracy of the plate theory, have the same displacement field in the middle plane/line, with no longitudinal displacement. The second family gathers all the other modes, denoted by non-bending (NB) modes, that appear at higher frequencies in the spectrum. Some of them are analogous to the longitudinal modes of the plate theory, with the same transverse displacement field in the middle plane/line. Other modes are also present, linked to 3D effects and thus with no counterpart in the beam/plate theory, with mode shapes dominated by thickness deformations. For any mode of this second family, the displacement field in the middle plane / line has only longitudinal components, and thus no transverse component. Examples of NB modes will be shown throughout the paper, especially in Table 2.

Let us decompose the displacement vector ${{{\varvec{X}}}}$ by denoting as $q_r$, $r=1, \ldots , N_B$ the bending coordinates, and $p_s$, $s=N_B + 1, \ldots , N$ the non-bending coordinates: ${{{\varvec{X}}}} = {[q_1, \ldots , q_{N_B}, p_{N_B+1}, \ldots , p_{N}]^{{\text {T}}}}$. Then, Eq. (7) can be rewritten for each coordinates [8, 13] and involves quadratic and cubic coupling terms between the $q_r$ and the $p_s$. We restrict ourselves to the case of a transverse low frequency excitation, for which the external forces remain normal to the middle plane of the plate. As a consequence, the dynamics is dominated by the bending modes, which are the only ones that receive external excitation. In this case, Eq. (7) can be simplified. First, all quadratic $\alpha _{ij}^r$ coefficients involving two bending coordinates i, j vanishes, in order to fulfil the symmetry of the restoring force [8, 39]. In addition, one can assume that the bending coordinates, which are directly excited by the external forcing, are considered of the order magnitude of a small parameter $\varepsilon$: $q_r=O(\varepsilon )$, for all $r\in \{1,\ldots , N_B\}$. On the other hand, NB coordinates, since they are not directly excited and shall vibrate at a lower order of magnitude, are assumed to scale as $\varepsilon ^2$: $p_s=O(\varepsilon ^2)$ for all $s\in \{N_B + 1,\ldots , N\}$. Plugging these two scaling in Eq. (7) and keeping only the leading order, one arrives to, for the bending coordinates, $\forall \; r=1, \ldots , N_B$:

$$\begin{aligned}&\ddot{q}_r+ 2\zeta _r\omega _r\dot{q}_r+ \omega _r^2q_r+\sum _{i=1}^{N_B}\sum _{l=N_B+1}^{N} \alpha _{il}^r q_ip_l\nonumber \\&\quad + \sum _{i=1}^{N_B}\sum _{j=i}^{N_B}\sum _{k=j}^{N_B} \beta _{ijk}^r q_i q_j q_k +O(\varepsilon ^4) = Q_r, \end{aligned}$$

(12)

and for the non-bending coordinates, $\forall \; s=N_B+1, \ldots , N$:

$$\begin{aligned} \ddot{p}_s+ 2\zeta _s\omega _s\dot{p}_s+ \omega _s^2p_s+\sum _{i=1}^{N_B}\sum _{j=i}^{N_B} \alpha _{ij}^s q_iq_j +O(\varepsilon ^3) = 0. \end{aligned}$$

(13)

In the above equation, one can notice that because of the transverse excitation, the second member of Eq. (13) is zero.

Equations (12,13) approximate the dynamics of the structure in the present case of a 3D FE model. If an analytic von Kármán model was used, those equations would be exact, with the terms $O(\varepsilon ^3)$ and $O(\varepsilon ^4)$ identically vanishing [8]. Those equations also show that the in-plane vibrations are quadratically coupled to the bending coordinates, while in the equations of motion of the transverse modes, only two nonlinear terms have to be taken into account: a quadratic coupling involving a product between a transverse and an in-plane coordinate, and a cubic term involving three transverse modes. This very specific form of equations renders the case of flat structures easier to solve than general shell problems that encompass all the possible nonlinear couplings as stated in Eq. (7).

2.3 Static condensation and nonlinear normal modes

Since the non-bending (NB) modes have natural frequencies very large as compared to those of the directly excited bending modes, $\omega _s\gg \omega _r$, the dynamical part of Eq. (13) can be neglected. The nonlinearities being more simple in this case, one can directly express the non-bending coordinate as function of the bending ones as:

$$\begin{aligned} p_s = -\sum _{i=1}^{N_B}\sum _{j=i}^{N_B} \frac{\alpha _{ij}^s}{\omega _s^2} q_iq_j. \end{aligned}$$

(14)

Substituting Eq. (14) into Eq. (12), one can rewrite the dynamics of the structure as a closed system involving only bending coordinates as

$$\begin{aligned} \ddot{q}_r+ 2\zeta _r\omega _r\dot{q}_r+ \omega _r^2q_r+ \sum _{i=1}^{N_B}\sum _{j=i}^{N_B}\sum _{k=j}^{N_B} \varGamma _{ijk}^r q_i q_j q_k = Q_r, \end{aligned}$$

(15)

where the cubic $\varGamma _{ijk}^r$ coefficients appear. Their general expression is derived in “Appendix A”.

If one is interested in deriving a reduced-order model for a single bending mode (say the master mode with label p), taking into account all the other non-bending mode, then Eq. (15) can be used and the leading nonlinear cubic term simply reduces to:

$$\begin{aligned} \varGamma _{ppp}^p = \beta _{ppp}^p - \sum _{s=N_B+1}^{N} \mathcal {C}^{ps}_{ppp} \end{aligned}$$

(16)

where the correction factors have been introduced and read:

$$\begin{aligned} \mathcal {C}^{ps}_{ppp}= \frac{\alpha _{ps}^p\alpha _{pp}^s}{\omega _s^2}=\frac{2(\alpha _{pp}^s)^2}{\omega _s^2}. \end{aligned}$$

(17)

These expressions show that the cubic term $\beta _{ppp}^p$ of the standard modal expansion must be corrected by the summation of all NB modes quadratically coupled to the master one.

The quadratic coefficients $\alpha ^p_{ij}$ have some symmetry relationships, provided that the nonlinear stiffness derives from a potential [23]. In particular, the following relationship holds: $\alpha _{ps}^p = 2 \alpha _{pp}^s$, which leads to the second equation (17). Recalling Eq. (10), the evaluation of $\alpha _{pp}^s$ requires only the computation of the nonlinear force ${{\varvec{f}}}_\text {nl}$ when the displacement along the pth linear mode is prescribed, whereas the calculation of the $\alpha _{ps}^p$ coefficients for each sth mode would require as many evaluation of ${{\varvec{f}}}_\text {nl}$ as the number of non-bending modes. Nevertheless, it requires the projection of nonlinear force onto each non-bending eigenvector $\phi _s$, and thus their computation, which is costly in practice since a large number of them is required to reach convergence (see Sect. 3.2).

From a physical perspective, following Eq. (10), coefficient $\alpha _{pp}^s$ can be seen as the projection onto the sth non-bending mode of the quadratic stiffness forces arising in the structure when a displacement along the linear pth mode is imposed. Interestingly, this quadratic coefficient is related to a monomial $q_p^2$ on the sth oscillator equation for $p_s$. These terms are recognized as invariant-breaking terms (see e.g. [42, 46]), in the sense that as soon as energy is given to the master mode p, all s modes having these important invariant-breaking terms will no longer be vanishing. These invariant-breaking terms are responsible for the loss of invariance of the linear eigenspaces, and they are found back naturally as correction factors when applying static condensation. They are also key in the formulation of invariant manifolds in order to define NNMs in phase space [33, 46].

In parallel to the static condensation emphasised here, one can use the reduction formulae given by the normal form approach, restricting the motion to a single Nonlinear Normal Mode (NNM) [42, 46, 47]. In this case, the reduced order model is directly constructed from Eq. (7). The main advantage as compared to the above described static condensation is that there is no need to assume the particular structure of the equations obtained for flat structures (Eqs. 12,13), thus generalizing the results to arches and shells. Considering only the NNM label p, the reduced-order model reads:

$$\begin{aligned}&\ddot{q_p} + \omega _p^2 q_p + \left( \sum _{s=N_B+1}^{N} -\dfrac{\alpha ^p_{ps}\alpha ^s_{pp}}{\omega _s^2}\left( \dfrac{\omega _s^2-2\omega _p^2}{\omega _s^2-4\omega _p^2}\right) + \beta _{ppp}^p\right) q_p^3\nonumber \\&\quad + \left( \sum _{s=N_B+1}^{N} \dfrac{\alpha ^p_{ps}\alpha ^s_{pp}}{\omega _s^2}\left( \dfrac{2}{\omega _s^2-4\omega _p^2}\right) \right) q_p \dot{q_p}^2 \; = \; 0 \; . \end{aligned}$$

(18)

Once again, one can observe that the correction brought to the cubic term $\beta ^p_{ppp}$ is solely given by the quadratic invariant-breaking terms. If one has been able to compute all the quadratic $\alpha ^p_{ij}$ coefficients appearing in Eq. (18), then the model can be used to simulate the dynamics. Also, it is worth mentioning that since the NB modes have high frequencies, we can assume that $\omega _s \gg \omega _p$ (which is equivalent to neglect the membrane inertia). Then the term in factor of $q_p^3$ in Eq. (18) exactly reduces to the one obtained with static condensation in Eq. (15). On the other hand, the term in factor of $q_p \dot{q_p}^2$ has no counterpart in static condensation, but is an order of magnitude smaller since it scales as $1/\omega _s^4$, so that both models are almost equivalent when a slow/fast decomposition can be assumed. This extends the results of [4], in which the term $q_p^3$ is chosen as the leading term for experimental identification purposes. A complete comparison of static condensation and nonlinear normal modes is also provided in [34] in the context of clarifying the implicit condensation and expansion method. Finally, one can note that the formula used in Eq. (18) have been obtained thanks to a normal form approach on the conservative system [46], but they can be extended in order to take into account the damping of the slave modes in the master coordinate ROM, hence accounting for a finer prediction of the losses [43], a feature that once again is not possible with static condensation.

3 STEP convergence with 3D elements

3.1 Test examples and direct computation of coefficients with the STEP

Fig. 1

Mesh used in the FE computations for the first two test cases on thin and thick beams

Table 1 Nonlinear dimensionless coefficients $\alpha _{il}^r$, $\alpha _{ij}^s$ and $\beta _{ijk}^r$ of the clamped beam, with non-zero and zero Poisson’s ratios ($\nu =0$ and $\nu = 0.3$). The modes considered in the coefficient indexes are the first two bending modes ($i=1,2$) and the second axial mode ($l=N_B+2$). The maximum of displacement amplitudes has been fixed at h/20 for these computations

In order to properly point out the convergence issues faced by using the modal basis as input prescribed displacements for the STEP, we consider the simple case of a clamped-clamped thin beam, shown in Fig. 1a, with length, thickness and width equal to $L=1~\text {m}$, $h=1~\text {mm}$, $b=50~\text {mm}$. The Young’s modulus is chosen as $E=210~\text {GPa}$. This particular geometry has been chosen to be thin (the thickness to length ratio is $10^{-3}$) to compare the results of the STEP computation to analytic values, obtained from a beam model with Euler-Bernoulli kinematics and von Kármán assumptions, see [8] where these comparisons have been more fully addressed. Table 1 presents the computations of nonlinear modal coupling coefficients obtained by the classical STEP with 3D and shell elements, compared to the analytical values. The two meshes used here consist of four node DKQ shell elements and twenty-node brick elements (HEX20), respectively. The computations are realized with the open software Code_Aster [6]. 100 elements in length, 4 elements in the width have been used for both meshes, with 2 elements in the thickness for the 3D mesh. This sufficiently refined mesh ensures that there is no convergence issue for the computations of the nonlinear coefficients.

The results clearly highlights the fact that using blindly 3D elements in a STEP computation with the modal basis leads to individual values of coupling coefficients that are far from their reference, analytical values. In particular, the cubic coefficients are largely overestimated and a strong dependence to the Poisson’s ratio is found with the 3D elements: the $\beta _{iii}^r$ are exactly twice the expected result with 3D elements and a zero Poisson’s ratio, but they become almost three times overestimated with $\nu =0.3$. On the other hand, using shell elements allows recovering the exact analytical result if selecting $\nu =0$, whereas a 10 % error is found for the same STEP computation with $\nu = 0.3$. These results clearly demonstrate that the modal basis as input for the STEP can be used safely with 2D elements but its extension to 3D elements is very problematic and should lead to large errors. As already noticed in the introduction, the problem comes from the fact the one uses eigenmodeshape functions as projection basis, but not from the calculation procedure itself.

3.2 Condensation of the cubic coefficient and frequency-response curves

Fig. 2

Nondimensional correction factor $\mathcal {C}^{1s}_{111}/\beta ^1_{111}$ (Eq. 17), associated to the first bending mode ($p=1$) and to all the other modes ($s=2,\ldots N$) of the thick beam testcase, over the nondimensional mode number s/N. The Poisson ratio is $\nu =0.3$

In the previous section, we showed that a direct calculation of individual coefficients leads to different values as compared to analytical results. However, of main importance is the prediction of the global behaviour of the structure, in a dynamical regime where modes are nonlinearly coupled and interacting together. In this section, we show how the static condensation presented in Sect. 2.3 can help to understand how the modes are coupled in order to define the hardening/softening behaviour of bending modes.

In order to shed light on the couplings arising between the modes, another test case is chosen. It is a thick beam, with the same length $L=1~\text {m}$ but with a square cross section with $h=b=30~\text {mm}$, as shown in Fig. 1b. The cross section was chosen square to be able to easily observe the 3D deformations of the cross section. The material properties are $E=210~\text {GPa}$ for the Young’s modulus and $\rho =7800~\text {kg/m}^3$ for the density. A coarse mesh of 15 HEX20 elements along the axis and $2\times 2$ in the cross section is chosen, to obtain full model with a reduced number of degrees of freedom (1287). The analyses on this beam are run in the software CodeAster [6]. For the sake of simplicity, we restrict attention to the convergence of the effective cubic coefficient of the first bending mode, $\varGamma _{111}^1$.

Fig. 3

Nondimensional correction factor $\mathcal {C}^{1s}_{111}/\beta ^1_{111}$ (Eq. 17), associated to the first bending mode ($p=1$) and to all the other modes ($s=2,\ldots N$) of the thick beam testcase, over the sorted nondimensional mode number s/N, where the imposed order is by decreasing correction factor. The Poisson ratio is $\nu =0.3$

Figure 2a shows the behaviour of the correction factor $\mathcal {C}^{1s}_{111}$ (defined in Eq. (17)), normalized by the cubic coefficient $\beta _{111}^1$, as a function of the mode number s, for the thick beam having 1287 dofs. This plot shows that the number of modes that are coupled to the first bending mode by invariant-breaking terms is very large, and uniformly distributed along the frequency spectrum. In order to facilitate the readings, the modes for which the correction factor is below 10$^{-15}$ have been sorted as negligible. In this family of modes that are not important, one find backs all the odd axial modes, for symmetry reason. On the other hand, all even axial modes are strongly coupled to the first one. The most surprising result is that if one does considers only axial modes, then only a few portion of the couplings will be revealed and taken into account. Indeed, Fig. 2a shows that there is a very large number of modes having very large frequencies, and still being strongly coupled to the master bending mode.

In order to check the independence of this behaviour from the mesh refinement, a second mesh of 20 elements on the axis and 3 $\times$ 3 elements on the section has been defined on the same beam geometry. Fig. 2b reports a very similar behaviour for this second test case, where the distribution of coupled modes is uniform along the whole set of modes.

Figure 3 shows the same data than Fig. 2, but with now the modes sorted by decreasing correction factor. It is possible to observe that, by choosing 10$^{-15}$ as a threshold for the significance of each contribution, a small percentage of modes, around 20%, is actually relevant. Consequently, the number of relevant modes depends on the mesh: the more refined it is, the more relevent modes are needed to reach convergence. Moreover, as seen on Fig. 2, these modes are spread over the entire spectrum, which would need the computation of all the eigenmodes of the structure, an operation impossible in practice for a complex structure with a larger number of dofs.

Table 2 Order of appearance in the basis (#), eigenfrequencies, correction factors and shapes of the most relevant modes coupled with the first bending mode in y-direction for the thick clamped-clamped beam. The colors scale the modulus of the displacement field and arrows display the axial displacement for axial mode 34. The Poisson ratio is $\nu =0.3$

In order to gain insight into these coupled modes, Table 2 shows the associated mode shapes, sorted according to the importance of their contribution in the correction factor, thus following Fig. 3a. The table shows the first nine eigenmode shapes, recalling in the first column their number of appearance when the modes are sorted according to the eigenfrequencies. One can observe that only one of these modes is a pure axial mode: the fourth one appearing in table 2, also being the 34th by order of increasing frequencies. All the other ones involve important deformation in the thickness of the beam. They are thus called thickness modes, their presence being the direct consequence of 3D effects. The second column of Table 2 displays the eigenfrequencies, showing that they all are high-frequency modes. The last column shows the deformation of the section, showing the importance of thickness deformation.

Fig. 4

Convergence to the full model solution with the increase of coupled modes taken into account in the reduced model. a Frequency response curve of the beam at center, in the vicinity of the first bending mode eigenfrequency; case of the thick beam with 1287 dofs. Red: full model, and convergence curves with increasing number of modes retained in the truncation: 1 mode (yellow), 4 modes (green), 9 modes (blue), all modes statically condensed (dashed light blue). Solution with one NNM in black. b Convergence of the evaluated corrected cubic stiffness coefficient $\varGamma ^1_{111}$, defined in Eq. (16), with increasing number of linear modes kept in the truncation

The frequency-response curves of the thick beam in the vicinity of its first eigenfrequency is investigated in order to illustrate how the static condensation and the NNM approach are able to retrieve the correct nonlinear behaviour. Figure 4a shows the comparison of the solutions obtained by continuation, for different reduced-order models and the full model solution. The latter has been obtained by solving all the degrees of freedom of the system with a parallel implementation of harmonic balance method and pseudo-arc length continuation [2]; the computation of the full forced response with 3 harmonics lasted approximately 36 hours. The convergence of the solution using static condensation with an increasing number of modes to compute the correction is also shown. Despite only few modes have a very high correction factor $\mathcal {C}^{s1}_{111}$, i.e. play a major role in the decrease of $\beta ^1_{111}$ (see Eq. (16)), it is the sum of the contributions from all the coupled modes that makes the reduced model converge to the solution of the full one. In Fig. 4b, the strong stiffening effect coming from not having included enough coupled modes, is slowly reduced by their inclusion in the basis; however, only the response obtained by static condensation of all coupled modes (cyan dashed) approximates the solution correctly (almost overlapped with the full model solution in red). On the other hand, the NNM solution with all the modes taken into account show also a direct convergence to the frequency-response curve.

Figure 4b shows the convergence of the corrected cubic coefficient $\varGamma ^1_{111}$ defined in Eq. (16) with the number of modes retained, i.e. the first mode plus the number of coupled modes condensed. When only one coupled mode is taken into account, the cubic coefficient $\beta ^1_{111}$ overestimates largely $\varGamma ^1_{111}$ (5.5 times): it is first explained by the classical bending-membrane coupling effect, and secondly to Poisson effect relating to the results given in Table 1. With 9 coupled modes the error on $\varGamma ^1_{111}$ is still significant (more than 60%). This strong overestimation of the cubic stiffness value results in the unrealistic stiffening effect observed in the forced response. The number of coupled modes that must be taken into account to ensure an acceptable accuracy makes the use of STEP in its first classic formulation (i.e. with the eigenmodeshape functions as projection basis and without condensation) quite impractical: 44 modes give a 1% error and 68 an error of 0.1%. The condensation of these modes onto the excited one becomes then a viable option to drastically reduce the computational burden without affecting the accuracy of the solution.

4 Alternative computational methods

In the previous section, we have shown that in the case of 3D elements, a strong coupling with thickness modes occurred, thus rendering the convergence of the modal ROM particularly stringent. When one is able to compute all the linear modes and associated coefficients, then static condensation and normal form approach can be used to finally produce accurate ROMs. However in most of the cases, the computation of all the linear modes, including the thickness modes appearing at very high frequencies, is out of reach. In this section, we investigate two alternative methods, for which there is no need to compute all the linear modes: static modal derivative, and a modified STEP.

4.1 Static modal derivatives

Sections 2 and 3 were devoted to the derivation of a reduced order model from a modal point of view. In fact, a modal projection of the quadratic nonlinear forces onto each mode $\phi _s$ is required to obtain the coefficients $\alpha ^s_{pp}$. Here we want to introduce the concept of static modal derivatives (see [12]) because its application provides the same results as the static condensation of all non-bending modes, but without requiring the computation of their associated eigenvectors.

The definition of modal derivatives have arisen from the recognition of the fact that in the nonlinear range, mode shapes and frequencies depend on amplitude [12]. Introducing this dependency in the eigenproblem defining the modes, one arrives at a quantity defined as the modal derivative [12, 51]. Following the definition of static modal derivatives ${{{\varvec{\theta }}}}_{pr}$ (SMD) from [13], it reads:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pr}=-{{{\varvec{K}}}}^{-1}\left. \left( \dfrac{\partial }{\partial q_p} \dfrac{\partial {{{\varvec{f}}}}_\text {nl}}{\partial {{{\varvec{x}}}}}(\varvec{\phi }_p q_p) \right) \right| _{q_p=0} \cdot \varvec{\phi }_r. \end{aligned}$$

(19)

When $p=r$ a more convenient expression (in the point of view of its direct computation from a FE code) for the modal derivative ${{{\varvec{\theta }}}}_{pp}$ writes:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pp}=-{{{\varvec{K}}}}^{-1}\left( \dfrac{{{{\varvec{f}}}}_\text {nl}(\lambda {{{\varvec{\phi }}}}_p)+{{{\varvec{f}}}}_\text {nl}(-\lambda {{{\varvec{\phi }}}}_p)}{\lambda ^2}\right) . \end{aligned}$$

(20)

The equivalent general expression for ${{{\varvec{\theta }}}}_{pr}$, with $p\ne r$ is provided in “Appendix B”. Eq. (20) shows how the SMD can be easily computed from a set of applied static displacement, in a manner having analogies with the STEP. The term in parenthesis in Eq. (20) can be seen as the numerical second order derivatives with respect to $\lambda$ of the nonlinear force along mode p, evaluated at the equilibrium position. This can be easily evaluated in any FE software by imposing, on a nonlinear structure, a displacement proportional to the linear pth mode with first positive and then negative sign in order to isolate the quadratic part of the nonlinear forces.

The SMD can thus be seen as an added displacement vector that enriches the basis constituted by the linear mode p, in a way that takes into account the nonlinear deformation of the structure. The use of SMDs as added vectors in the reduced order model basis is extensively documented in [13, 30, 35, 40, 51], and a complete comparison of quadratic manifolds derived from modal derivatives with normal for theory is given in [48]. Here the focus is on the relationship between modal derivatives and non-bending modes and on the equivalence between static condensation of all non-bending modes and static condensation of modal derivatives.

Given a system with geometric nonlinearities up to cubic order, the static modal derivatives related to the pth and rth linear modes can then be expressed in terms of the vector of quadratic coefficients ${{{\varvec{\alpha }}}}_{pr}$ as:

$$\begin{aligned} {{\varvec{\theta }}}_{pr}=-{{{\varvec{V}}}}{{{\varvec{\varOmega }}}}^{-2}\;{{\varvec{\alpha }}}_{pr}\, , \end{aligned}$$

(21)

and the one relative to the pth linear mode as:

$$\begin{aligned} {{\varvec{\theta }}}_{pp}=-2\;{{{\varvec{V}}}}{{{\varvec{\varOmega }}}}^{-2}\;{{\varvec{\alpha }}}_{pp}\, , \end{aligned}$$

(22)

where the full matrix of eigenvectors ${{{\varvec{V}}}}$ has been introduced. The detailed derivation of these two equations is given in “Appendix B” together with the classical orthonormality properties of the matrix of eigenvectors ${{{\varvec{V}}}}$. By expanding Eq. (22) over all the modes and by noticing that, for a flat structure, $\alpha ^s_{pp}$ is non-zero only when s is a non-bending mode, one obtains this important relationship (see “Appendices B and [48])”:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pp}=-\sum ^{N}_{s=N_B+1}2{{{\varvec{\phi }}}}_s \dfrac{\alpha ^s_{pp}}{\omega ^2_s}. \end{aligned}$$

(23)

The SMD thus appears as a linear combination of coupled modes with factor $-2\alpha ^s_{pp}/\omega ^2_s$; therefore, it can be seen as a displacement field that takes into account the contribution of all non-bending modes into one equivalent vector. To show this property, the SMD relative to the first bending mode of the thick beam test case is depicted in Fig. 5a. The projection of the SMD onto the linear modes of the system recovers Eq. (23). The contributions from the non-bending modes that have been identified in previous calculations, the fourth axial mode as well as various thickness modes, appear in the SMD. For each mode s, The modal amplitudes $q_s$ obtained from the projection coincides with those from Eq. (23) i.e. they are equal to $-2{\alpha ^s_{pp}}/{\omega ^2_s}$. These results recover and elaborate on those obtained in [13, 35], where it was shown that axial modes are contained in the SMDs of bending modes. The result is here extended to thickness modes and specified since the exact participation factor of each mode is made explicit.

Fig. 5

a Static modal derivative ${{\varvec{\theta }}}_{11}$ associated to the first bending mode. b Four first non-bending modes contained in the SMD ${{\varvec{\theta }}}_{11}$, found equivalently by projecting ${{\varvec{\theta }}}_{11}$ on the linear mode basis, or by application of Eq. (23). The relative modal participation factors $q_s$ of each of these modes numbered 34, 167, 328 and 330 (see also Table 1) is also numerically given, normalized by the total amplitude of the SMD ${{\varvec{\theta }}}_{11}$ ($q_{SMD}=1$), and exactly recovers the factors exhibited in Eq. (23). (Color figure online)

Once understood that the SMD allows gathering in a single vector the participation of all coupled modes, we want to show how to retrieve directly, from the calculation of the SMD, the correct nonlinear behaviour of the structure, when the motion is restricted to a single master mode. In the specific case of a flat structure, Eq. (14) shows that the amplitudes of the NB modes can be explicitely related to the squared amplitude of the single master mode labeled p, thanks to the static condensation, as:

$$\begin{aligned} p_s = - \dfrac{\alpha ^s_{pp}}{\omega ^2_s}q_p^2. \end{aligned}$$

(24)

The physical displacement ${{\varvec{x}}}(q_p)$ that corresponds to the solution gathering together the master bending mode p and all its coupled NB modes can be written as:

$$\begin{aligned} {{\varvec{x}}}(q_p)=q_p{{\varvec{\phi }}}_p+\sum ^{N}_{s=N_B+1}p_s\;{{{\varvec{\phi }}}}_s. \end{aligned}$$

(25)

In this last equation, replacing $p_s$ by its value obtained from static condensation, Eq. (24), and then using Eq. (23) defining ${{{\varvec{\theta }}}}_{pp}$ as a summation on the NB modes, one arrives easily at the fact that this physical displacement can be expressed as a function of the modal coordinate plus the participation of the SMD as:

$$\begin{aligned} {{\varvec{x}}}(q_p)=q_p{{\varvec{\phi }}}_p+\frac{1}{2}q_p^2{{\varvec{\theta }}}_{pp}\;. \end{aligned}$$

(26)

If one wants now to derive a reduced-order model composed of a single master coordinate (say $q_p$ here) and that contains the correct nonlinear behaviour, then the equation of motion would simply read:

$$\begin{aligned} \ddot{q}_p+ 2\zeta _p\omega _p\dot{q}_p+ \omega _p^2q_p+ \tilde{\varGamma }_{ppp}^p q_p^3 = Q_p, \end{aligned}$$

(27)

with $\tilde{\varGamma }_{ppp}^p$ a corrected cubic coefficient. Thanks to Eq. (11a), one knows that a cubic coefficient can be found from this computation, provided the imposed displacement is selected correctly. If the imposed displacement is along a linear mode, as in Eq. (11a), then one will retrieve the modal nonlinear coupling coefficient, but other choice of imposed displacement can be made. In particular, if one selects the one given by Eq. (26), then the cubic coefficient $\tilde{\varGamma }_{ppp}^p$ will contain the contribution of the master mode plus that of the SMD. Since Eq. (23) shows that in our particular case (flat structure, one master mode), the SMD is completely equivalent to the static condensation of all coupled NB modes, then one will easily understand that $\tilde{\varGamma }_{ppp}^p = \varGamma _{ppp}^p$, the corrected cubic coefficient given in Eq. (16). The complete proof of this result is provided in “Appendix C”.

The main result, from the SMD perspective, is that if one restricts to a single master mode p, then the SMD can be easily computed thanks to Eq. (20). Then the corrected cubic coefficient can be directly computed from:

$$\begin{aligned} \varGamma ^p_{ppp}={{{\varvec{\phi }}}_p^{{\text {T}}}}\left( {{{\varvec{f}}}}_\text {nl}({{\varvec{x}}}(\lambda ))-{{{\varvec{f}}}}_\text {nl}({{\varvec{x}}}(-\lambda ))\right) /2\lambda ^3, \end{aligned}$$

(28)

where the imposed displacement is selected as in Eq. (26). This procedure allows then to find exactly the same corrected cubic coefficient of the master mode, but without resorting to the computation of all eigenvectors, as needed in the static condensation. It is thus a much more computationally efficient to use this methodology. Numerical examples are provided in Sect. 6.1.

4.2 A modified STEP for 3D elements

Fig. 6

Examples of prescribed displacement field of the M-STEP, in the transverse direction and in the middle surface/line of a plate/beam

As observed in the previous sections, the modal ROM associated to 3D FE discretization shows a slow convergence because of the couplings with very high frequency modes involving thickness deformations. Since these thickness modes are a peculiarity of the 3D model, they have no counterpart in plate or beam theory, which concentrate the kinematical description on the middle plane/line. Also, using the STEP with plate/beam elements show a faster convergence since one has to recover only the well-known coupling between bending and in-plane motions. In order to circumvent these difficulties, we propose here to modify the STEP by prescribing the displacements only on the middle line/plane of the structure, and to let free the other degrees of freedom (Fig. 6). The idea is to include automatically the effects of NB modes, by a kind of implicit condensation of their motion, embedded into the prescribed displacement on the middle line/plane. The obtained method is called the M-STEP, for Modified-STEP. Note that a comparable idea has also been introduced in [15, 49], but for 2D flat structures only, where only transverse motions were prescribed, leaving the other degrees of freedom free and thus building directly a condensed model.

4.2.1 Formulation

We show in this section that it is possible to compute directly the cubic coefficients $\varGamma _{ijk}^r$ of Eq. (15) with a modified STEP, without having to compute all the coefficients $\alpha _{ij}^k$ and $\beta _{ijl}^k$ beforehand. Note that in the classical STEP, the three components of the displacement field are prescribed to all the nodes of the FE mesh: the whole vector of unknown ${{{\varvec{x}}}}$ is imposed and the FE code computation is just an evaluation of the internal force vector ${{{\varvec{f}}}}_e={{{\varvec{f}}}}({{{\varvec{x}}}})$.

Here, we choose to apply the STEP by prescribing the displacement field only to selected nodes and for selected components. Precisely, we denote by $\mathcal {S}$ the middle plane/line of the structure and ${{{\varvec{n}}}}$ the bending direction. To compute $\varGamma _{ppp}^k$ for a given $p\in \{1,\ldots N_B\}$, we choose to perform a FE computation by prescribing (i) only the transverse component of the bending mode ${{{\varvec{\phi }}}}_p$ (ii) only on the nodes of the middle plane/line of the structure (Fig. 6). We then prescribe the following time independent displacement to the structure:

$$\begin{aligned} {{{\varvec{x}}}}|_{\mathcal {S},{{{\varvec{n}}}}} = \lambda \, {{{\varvec{\phi }}}}_p|_{\mathcal {S},{{{\varvec{n}}}}},\quad \lambda \in \mathbb {R}. \end{aligned}$$

(29)

where ${{{\varvec{x}}}}|_{\mathcal {S},{{{\varvec{n}}}}}$ corresponds to displacement ${{{\varvec{x}}}}$ restricted to (i) the nodes of the FE mesh belonging to $\mathcal {S}$ and to (ii) its components along the direction defined by ${{{\varvec{n}}}}$. In all the other nodes and directions, a zero forcing is prescribed. We then use the finite element code to solve this problem and obtain ${{{\varvec{x}}}}$ as well as ${{{\varvec{f}}}}({{{\varvec{x}}}})={{{\varvec{f}}}}_e$ everywhere. Since some components of ${{{\varvec{x}}}}$ are not prescribed, a Newton-Raphson procedure is necessary to solve this nonlinear algebraic problem.

To precise the method, we call master dofs the ones for which the displacement is prescribed, and slave dofs the other ones. The full displacement vector ${{{\varvec{x}}}}$ and the internal forcing ${{{\varvec{f}}}}$ can thus be decomposed as:

$$\begin{aligned} {{{\varvec{x}}}}= \begin{bmatrix} {{{\varvec{x}}}}_{\text {M}}\\ {{{\varvec{x}}}}_{\text {S}} \end{bmatrix},\qquad {{{\varvec{f}}}}= \begin{bmatrix} {{{\varvec{f}}}}_{\text {M}}\\ {{{\varvec{f}}}}_{\text {S}} \end{bmatrix}, \end{aligned}$$

(30)

where the index $\text {M}$ and $\text {S}$ are associated to the master and slave dofs, respectively. The M-STEP consists in prescribing:

$$\begin{aligned} {\left\{ \begin{array}{ll} {{{\varvec{x}}}}_{\text {M}}=\lambda {{{\varvec{\phi }}}}_p^\text {M}, \\ {{{\varvec{f}}}}_{\text {S}}={{{\varvec{0}}}}. \end{array}\right. } \end{aligned}$$

(31)

Then, solving the static problem ${{{\varvec{f}}}}({{{\varvec{x}}}})={{{\varvec{f}}}}_e$ with the FE code leads to compute the internal force vector ${{{\varvec{f}}}}_\text {M}({{{\varvec{x}}}})$ on the master nodes and the displacement ${{{\varvec{x}}}}_S$ on the slave nodes. The solution of the problem then reads:

$$\begin{aligned} {{{\varvec{x}}}}= \begin{bmatrix} \lambda {{{\varvec{\phi }}}}_p^\text {M}\\ {{{\varvec{x}}}}_{\text {S}} \end{bmatrix},\qquad {{{\varvec{f}}}}={{{\varvec{f}}}}_e= \begin{bmatrix} {{{\varvec{f}}}}_{\text {M}}\\ {{{\varvec{0}}}} \end{bmatrix}. \end{aligned}$$

(32)

Translated in the modal space, the above computations are close to the following situation. Prescribing via ${{{\varvec{x}}}}$ only the transverse motion in the form of ${{{\varvec{\phi }}}}_p$, and because ${{{\varvec{\phi }}}}_p$ is orthogonal to the other bending modes ${{{\varvec{\phi }}}}_k$, $k\ne p$, the modal coordinate are $q_p\simeq \lambda$ and $q_k\simeq 0$. Considering the orthogonality relations associated to the stiffness matrix, this leads to assume that, for all $s\in \{1,\ldots N_B\}, \; s\ne p$:

$$\begin{aligned} {{{{\varvec{\phi }}}}_p^{{\text {T}}}}{{\varvec{K}}}{{{\varvec{x}}}}\simeq {{{{\varvec{\phi }}}}_p^{{\text {T}}}}{{\varvec{K}}}(\lambda {{{\varvec{\phi }}}}_p)=\lambda \omega _p^2 m_p,\quad {{{{\varvec{\phi }}}}_s^{{\text {T}}}}{{\varvec{K}}}{{{\varvec{x}}}}\simeq {{{{\varvec{\phi }}}}_s^{{\text {T}}}}{{\varvec{K}}}(\lambda {{{\varvec{\phi }}}}_p)=0. \end{aligned}$$

(33)

In other words, it is assumed that the nonzero slave part ${{{\varvec{x}}}}_\text {S}$ of ${{{\varvec{x}}}}$ is not involved in the orthogonality relations of the bending modes. Moreover, since the part of the displacement associated to longitudinal and thickness displacements is not prescribed by ${{{\varvec{x}}}}$, the associated NB modal coordinates are not zero and their value depend on the nonlinear coupling and the geometric nonlinearities. Finally, because any NB eigenmode ${{{\varvec{\phi }}}}_s$ has zero displacements on the middle plane/line in the transverse direction, the modal forcing of the NB modes is exactly zero:

$$\begin{aligned} {{{\varvec{\phi }}}}_s= \begin{bmatrix} {{{\varvec{0}}}} \\ {{{\varvec{\phi }}}}_s^{\text {S}} \end{bmatrix}\quad \Rightarrow \quad Q_s=\frac{{{{{\varvec{\phi }}}}_s^{{\text {T}}}}{{{\varvec{f}}}}_e}{m_s}=0. \end{aligned}$$

(34)

To summarize, one has:

Error estimates of those assumptions will be introduced in Sect. 4.2.2.

Using the assumptions (35) in Eqs. (12,13) leads to:

The above Eq. (38) leads to:

$$\begin{aligned} p_s =-\dfrac{ \alpha ^s_{pp} \lambda ^2}{\omega ^2_s }, \end{aligned}$$

(39)

that can be condensed into (36), (37) to give:

$$\begin{aligned} \left( \beta ^p_{ppp}-\sum _{l=N_B+1}^{N} \frac{\alpha ^p_{ps} \alpha ^s_{pp}}{\omega ^2_s } \right) \lambda ^3&= {{{{\varvec{\phi }}}}_p^{{\text {T}}}}{{{\varvec{f}}}}({{{\varvec{x}}}}) /m_p - \lambda \omega _p^2, \end{aligned}$$

(40)

$$\begin{aligned} \left( \beta ^k_{ppp}-\sum _{l=N_B+1}^{N} \frac{\alpha ^k_{ps} \alpha ^s_{pp}}{\omega ^2_s } \right) \lambda ^3&= {{{{\varvec{\phi }}}}_k^{{\text {T}}}}{{{\varvec{f}}}}({{{\varvec{x}}}}) /m_k,\quad \forall k\ne p. \end{aligned}$$

(41)

One then recognizes the terms in parenthesis as the sought corrected cubic coefficients $\varGamma ^i_{ppp}$, defined by Eqs. (16) and (51), so that:

$$\begin{aligned} \varGamma ^p_{ppp}=\frac{{{{{\varvec{\phi }}}}_p^{{\text {T}}}} {{{\varvec{f}}}}({{{\varvec{x}}}})}{m_p \lambda ^3}-\frac{\omega _p^2}{\lambda ^2},\qquad \varGamma ^k_{ppp}=\frac{{{{{\varvec{\phi }}}}_k^{{\text {T}}}} {{{\varvec{f}}}}({{{\varvec{x}}}})}{m_k \lambda ^3}. \end{aligned}$$

(42)

Consequently, all $\varGamma _{ppp}^i$, $i=1,\ldots N_B$ relies on a single nonlinear static finite elements computation, defined by Eq. (45). Other coefficients $\varGamma _{ijl}^k$ can be obtained in the same manner, by mixing different $x_\text {M}$ on several bending modes, following the classical STEP.

This above described M-STEP method is a way of automatically embed in the computation the effect of all the NB modes nonlinearly coupled to the bending modes associated to $\varGamma _{ijl}^k$. The essence of the method is to select the prescribed displacement ${{{\varvec{x}}}}$ so that (i) it leaves free the degrees of freedom associated to the NB modes, so that the forcing $Q_s$ of the longitudinal modal coordinates in Eq. (38) is exactly zero and (ii) it is as orthogonal as possible to the other bending modes than the pth.

4.2.2 Quality indicator for the convergence of the method

In order to be able to quantify a priori the quality of the computation, a main idea is to check the validity of the assumptions used in the two first equations (35), which are true at first order but might deteriorate in case of an incorrect selection of master dofs. Equivalently, one can verify the orthogonality of the displacement vector ${{{\varvec{x}}}}$ to the bending modes written in Eq. (33). To that purpose, let us define the following errors:

$$\begin{aligned} e_1^{pp} = \frac{{{{{\varvec{\phi }}}}_p^{{\text {T}}}}{{\varvec{K}}}{{{\varvec{x}}}}}{\lambda \omega _p^2 m_p}-1,\qquad e_1^{pk} = \frac{{{{{\varvec{\phi }}}}_k^{{\text {T}}}}{{\varvec{K}}}{{{\varvec{x}}}}}{\lambda \omega _p^2 m_p},\quad \forall k\ne p, \end{aligned}$$

(43)

that should be small as compared to 1 because of Eq. (33). If one wants to compute those errors with a FE in a non intrusive way, ${{\varvec{K}}}{{{\varvec{x}}}}={{{\varvec{f}}}}_1$ can be computed as the reaction force vector ${{{\varvec{f}}}}_1$ of a linear static computation where ${{{\varvec{x}}}}$ is prescribed to all the nodes of the FE mesh.

Fig. 7

Evolution of the ratio of the cubic coefficients $\beta _{iii}^{i}$ compared with the analytic values $\beta _{iii}^{i,AN}$, with regard to the Poisson ratio, for the first four ($i=1,\ldots 4$) bending modes of the thin beam (Fig. 1a) meshed with shell or 3D elements, as specified on the plot. The heuristic dependences on the Poisson’s ratio are plotted with black dashed lines. All symbols are merged

Another check can also be performed by prescribing Eq. (31) into a linear static computation:

$$\begin{aligned} {{{\varvec{K}}}}{{{\varvec{x}}}}_\text {l} = {{{\varvec{f}}}}_e, \end{aligned}$$

(44)

which gives:

$$\begin{aligned} {{{\varvec{x}}}}_\text {l}= \begin{bmatrix} \lambda {{{\varvec{\phi }}}}_p^\text {M}\\ {{{\varvec{x}}}}_{\text {lS}} \end{bmatrix}. \end{aligned}$$

(45)

Since there are no geometrical nonlinearities, imposing ${{{\varvec{\phi }}}}_p$ on the middle line/surface in the transverse direction should result in a vector almost collinear to ${{{\varvec{\phi }}}}_p$, that is ${{{\varvec{x}}}}_\text {l}\simeq \lambda {{{\varvec{\phi }}}}_p$. In particular, the slave part ${{{\varvec{x}}}}_{\text {lS}}$ of ${{{\varvec{x}}}}_{\text {l}}$ should be very close to ${{{\varvec{\phi }}}}_p^\text {S}$, the slave part of ${{{\varvec{\phi }}}}_p$. We then define the following error:

$$\begin{aligned} e(\hat{{{{\varvec{y}}}}},{{{\varvec{y}}}})=\frac{||\hat{{{{\varvec{y}}}}}-{{{\varvec{y}}}}||}{||{{{\varvec{y}}}}||}, \end{aligned}$$

(46)

where $||\cdot ||$ is the norm of vector $\cdot$, and we check that $e_2^p = e({{{\varvec{x}}}}_l,\lambda {{{\varvec{\phi }}}}_p)$ is very small as compared to 1.

5 Physical mechanisms of the nonlinear couplings

5.1 Poisson effect

Fig. 8

Convergence without Poisson effect ($\nu =0$) for the thick beam. a Nondimensional correction factor $\mathcal {C}^{1s}_{111}/\beta ^1_{111}$ of all the modes over the nondimensional mode number s/N. b Convergence of the evaluated corrected cubic stiffness coefficient $\varGamma ^1_{111}$, defined in Eq. (16), with increasing number of linear modes kept in the truncation. c Mode shapes and order of appearance in the basis (#) of four of the most relevant modes coupled with the first bending mode

The previous sections show that in the case of a 3D model, a given bending mode is nonlinearly coupled to numerous high frequency modes, most of them involving thickness deformations. To understand this effect, a numerical study of the sensitivity of the coefficients computed with the STEP on the Poisson ratio is here given. It is found that a precise dependence of the cubic coefficient $\beta ^r_{iii}$ can be established: in Fig. 7, the values obtained by the direct application of the STEP are fitted to an heuristic law related to the Poisson ratio. In particular, the growths of the cubic coefficients in the case of shell and 3D elements match perfectly the ratios:

$$\begin{aligned} \rho _1=\frac{1}{1-\nu ^2},\qquad \rho _2=\frac{2}{ (1+\nu )(1-2\nu )}, \end{aligned}$$

(47)

related to 2D and 3D constitutive laws. Indeed, the 3D constitutive law for an isotropic elastic material writes:

$$\begin{aligned} {{\varvec{\pi }}} = \frac{E}{(1+\nu )(1-2\nu )}\big [\nu \, {\text {tr}}({{\varvec{\varepsilon }}}) {{\varvec{I}}}_3 + (1-2\nu ) {{\varvec{\varepsilon }}}\big ], \end{aligned}$$

(48)

where ${{\varvec{\pi }}}$ denotes the second Piola-Kirchhoff stress tensor, ${{\varvec{\varepsilon }}}$ the Green-Lagrange strain tensor, ${{\varvec{I}}}_3$ the identity operator in 3D and $(E,\nu )$ the Young’s modulus and the Poisson ratio of the material.

Moreover, in the case of usual plate theories, a plane stress state is assumed, for which the transverse component $\pi _{zz}=0$ of the stress tensor is zero (z being the direction normal to the middle plane of the plate). In this case, the in-plane counterpart of (48) reads:

$$\begin{aligned} \tilde{{{\varvec{\pi }}}} = \frac{E}{1-\nu ^2}\big [\nu \, {\text {tr}}(\tilde{{{\varvec{\varepsilon }}}}) {{\varvec{I}}}_2 + (1-\nu ) \tilde{{{\varvec{\varepsilon }}}}\big ], \end{aligned}$$

(49)

where $\tilde{{{\varvec{\pi }}}}$ and $\tilde{{{\varvec{\varepsilon }}}}$ denote the plane parts of ${{\varvec{\pi }}}$ and ${{\varvec{\varepsilon }}}$, respectively, and ${{\varvec{I}}}_2$ the identity operator in 2D.

Equations (48) and (49) makes directly appear the ratios $\rho _1$ and $\rho _2$, which suggests also a direct relationship between the constitutive law and the results presented in Fig. 7. As a side note, the reference values, used to compare nondimensional values of cubic coefficients in Fig. 7, differ from those in Figs. 4b and 8b. Indeed, in Fig. 7, the normalising coefficient has been selected as $\beta _{111}^{1,AN}$, the analytical value obtained from the beam theories (see e.g. [8]). As shown in the previous sections, this coefficient has to be compared with the corrected coefficient used once the convergnce is obtained from 3D models. On the other hand, in Fig. 4b and 8b, the reported values are normalized with respect to $\beta _{111}^{1}$, i.e. the uncorrected cubic coefficient, the one coming from direct application of STEP and shown in Table 1. This explains the different values observed between the two figures (e.g. $\beta _{111}^{1}/\beta _{111}^{1,AN} \approx 3.8$ for $\nu = 0.3$ on Fig. 7, whereas $\beta _{111}^{1}/\varGamma _{111}^{1} = 5.5$ on Fig. 4b).

5.2 Geometrical nonlinearities

In this section, we focus on the physical explanation of the nonlinear couplings with thickness modes, whose origin is the geometrical nonlinearities. Considering first Fig. 7, it can be inferred that the couplings are amplified by the Poisson ratio, but that they are present even without Poisson effect, since there is a factor 2 between the FE value of $\beta _{iii}^{i}$ with respect to its corresponding analytical value in the case $\nu =0$. This leads us to investigate the couplings in this particular case.

Figure 8 is the analog, with $\nu =0$, of Figs. 3a, 4b and Table 2. Comparing those figures shows that the number of coupled modes is much smaller in the case $\nu =0$ than for $\nu =0.3$: the relevant modes correspond to 5% of the modal basis if $\nu =0$, whereas it was 20% for $\nu =0.3$ (see the modes with a correction factor above $10^{-15}$ in Figs. 3a and 8a). Moreover, the deformations of the cross section in the case $\nu =0$ are purely in the bending transverse-y direction in the case $\nu =0$ (y, as defined in Figs. 1 and 14, is the deformation direction of the first bending mode considered in all computations of the present article), whereas they were more complex (in 2D) in the case $\nu =0.3$ (compare the mode shapes of Table 2 and Fig. 8). Finally, taking a close look at the static modal derivative (SMD), shown in Fig. 5a, that gathers all the corrections brought by the NB modes, shows that it has almost the same shape in both cases $\nu =0$ and $\nu =0.3$. In particular, Fig. 9 shows that the deformations of the SMD cross section occur without distorsion: initially a square, it is deformed in a rectangle. In the case of no Poisson’s ratio, the deformation is purely in the bending y-direction, whereas the Poisson’s effects adds a slight deformation in the lateral z-direction.

Those effects are purely geometrical and come from (i) the particular 3D shape of the eigenvectors and (ii) how these particular shapes, resulting from a linear computation, are modified by the geometrical nonlinearities. A closed form solution in a simple bending case is exposed in “Appendix D”. It shows that the 3D shape of the eigenmodes is the combination of three contributions (see Eq. (89)):

• the deformation of the neutral axis/plane of the structure, described by classical beam/plate theories;

• the 3D rotation of the cross section around the z-axis, that produces an axial deformation with a linear dependence in the thickness coordinate y;

• the 3D Poisson effect, that distorts the cross section in its two (y and z) directions.

Then, by computing the Green-Lagrange strain tensor with this particular linear deformation, it is shown that the geometrical nonlinearities adds two contributions to the classical von Kármán beam/plate model:

• 3D effects that are independent of the Poisson effect, that explains a stretching in the transverse y direction, without any deformation in the lateral z direction. Those effects are a direct consequence of the 3D rotation of the cross section created by the bending;

• 3D Poisson effect, that involve stretching in both the transverse y and lateral z directions.

Those two geometrical effects are purely 3D and are additional to the classical membrane/bending coupling.

Fig. 9

Static modal derivative and wiew of the cross-section in the undeformed and deformed configurations, with both $\nu =0$ and $\nu =0.3$

Having in mind those observations, we can now explain those nonlinear thickness coupling. We have first to remark that in both cases of a beam(1D)/plate(2D) von Kármán model and the present 3D model, the modal expansion of Eq. (7) is exact, provided N is the number of degrees of freedom of the model. Moreover, Table 3 and Fig. 4a show that all models converge to the same solution, which proves that the relevant modes of the basis combine themselves in different ways to give, at the end, the same solution. In fact, as a consequence of the above observations, the thickness modes are here to geometrically compensate (i) the nonlinear deformations in the transverse bending direction due to the 3D rotation of the cross sections, shown in blue on Fig. 9 and (ii) the additional deformations of the cross section due to the Poisson effect. Looking again at the deformed cross section shown in red in Fig. 9, one can understand that at the end, the complex 3D distorsions of the cross section due to the Poisson effect (shown in Fig. 14c) must be fully compensated by the NB modes, which then need to be numerous. This explains the bad convergence of the modal expansion, even worse in the case of a non zero Poisson ratio (5% of the modal basis if $\nu =0$, and 20% for $\nu =0.3$).

6 Numerical examples

Fig. 10

First-mode displacements prescribed on the middle line of the clamped beam and the middle surface of the circular plate. For visualisation purpose, the mesh of the plate is less fine than the one use for the computations of the ROM coefficients

In this section, numerical examples on the three different strategies proposed in order to overcome the bias observed when using 3D elements, are given. In each case, the dominant cubic coefficient of a single mode is compared, using either the M-STEP, the static condensation of all the coupled linear modes, or the static modal derivative. Two test cases are used for the comparisons: the clamped-clamped beam of Fig. 1b, and a clamped circular plate.

6.1 Application to a clamped-clamped beam

Computations of the condensed coefficients $\varGamma _{ijk}^r$ are first performed with the three methods. For the M-STEP, the prescribed displacement is depicted in Fig. 10. The comparison made in Table 3 attests that the condensation with all the eigenmodes and the first static modal derivative are quasi equivalent, whereas the M-STEP gives very close values. The relative errors are very small ($<0.5$ %) in each case. The analytical reference values present slightly larger errors, between 1 % and 5 %, probably due to unavoidable differences between the analytical beam theory and the numerical 3D computation. This could be explained by the aspect ratio $h/L=0.03$ of the beam, which is not so small to fully verify Euler-Bernoulli assumptions.

Table 3 Corrected cubic coefficients with the three condensation methods and comparison with anaytical values. Due to the symmetry of the mode shapes, the nonlinear coefficients which couple the first and the third modes have all nonzero values, and are therefore chosen for these computations

Fig. 11

a Dependence of the cubic coefficient $\varGamma _{iii}^i$ with regard to the prescribed displacement amplitude, for different modes i, with $i=1,2,3,4$. Black dashed line : reference analytical value ($\varGamma _{iii}^{i,FE} = \varGamma _{iii}^{i,AN}$), red dashed lines: limits of the range of validity. b Dependance of the cubic coefficient $\varGamma _{111}^1$ with regard to the prescribed displacement amplitude for different lines where the displacements are prescribed. Black dashed line : reference value ($\varGamma _{111}^{1,FE} = \varGamma _{111}^{1,AN}$). The inset shows the location, in the cross section of the beam, of the lines in which the master displacement is prescribed. (Color figure online)

In the case of the M-STEP with the displacement field prescribed on the neutral fiber, Fig. 11a gives the sensitivity to the prescribed displacement amplitudes of the corrected cubic coefficients $\varGamma _{iii}^i$ for different modes i. It is shown that a range of validity for the displacements amplitude centered around $\text {max}(\lambda \phi _p) \simeq h/2$ can be defined. As it could be expected, the length of this validity range, defined on Fig. 11a by a relative error smaller than $3\%$, decreases with the mode order.

Figure 11b shows what happens if the M-STEP is applied with different selections of the master degrees of freedom. We tried to prescribe the displacement field on three other lines of the beam, parallel to but different than the neutral fiber, as shown in the inset of Fig. 11a. We can conclude that the neutral fiber seems the best choice and that the upper line gives very close results. On the other hand, the values obtained when the displacements are prescribed on one of the lateral lines are far from the reference values. This feature will be analyzed in the following considering error estimators.

Fig. 12

Dependence of the criterion $e_1^{pi}$ with regard to the prescribed displacement amplitude when displacements are prescribed on the neutral fiber. The values are presented on a for $p=1$ and different modes $i=1,2,3,4,5$. On b, the values $e_1^{ii}$ are presented, also for $i=1,2,3,4,5$

Table 4 Values of the relative error $e_1^{pi}$ for different p, i, and master dofs. Nonlinear computations are performed with

The error estimate $e_1^{pi}$, introduced in Sect. 4.2.1, is now analyzed. The criterion $e^{pi}_1$ is first computed with the master transverse dofs prescribed in the neutral fiber, a single master mode $p=1$ and different transverse modes $k\in \{1,\ldots ,5\}$. The values presented in Fig. 12a show that the orthogonality is well verified, until the upper limit of validity range observed on Fig. 11a, after which the values $e_1^{13}$ and $e_1^{15}$ deviate from 0. For $p>1$, the coefficients $e_1^{pp}$ evolve in a similar way as $e_1^{11}$, as depicted on Fig. 12b. Table 4 give the numerical values of $e_1^{1i}$ in the plateau of Fig. 12a, proving that the error is the order of $10^{-4}$. Consequently, the orthogonality of the prescribed displacement field ${{{\varvec{x}}}}$ to the transverse modes ${{{\varvec{\phi }}}}_i$ is well verified, thus validating assumptions (35a,b).

Then, the same error estimate $e_1^{pi}$ is computed in the cases of a master displacement prescribed in the other lines of the inset of Fig. 11b. The obtained values are given in Table 4. In this cases, we quantitatively confirm the observation linked to Fig. 11b: the orthogonality of the displacement are not verified when the master dofs are placed on a line of the lateral surfaces of the beam. In particular, the values of the criterion $e_1^{pi}$ presented with $p=1$ and $i=1,2,3,5$ in Table 4 highlight a loss of orthogonality between the first and the odd modes $i=1,3,5$, in the case of master dofs on a lateral line: indeed, the values of $e_1^{11}$, $e_1^{13}$ and $e_1^{15}$ deviate from 0.

The second error estimate $e_2^1$, also introduced in Sect. 4.2.1 and linked to an estimate of the collinearity of the prescribed displacement ${{{\varvec{x}}}}_l$ to the master transverse modes ${{{\varvec{\phi }}}}_1$, in the case of a linear computation (see Eq. (44)). This estimate confirms the above results, in particular that the displacements must preferentially be prescribed on the neutral line.

A physical explanation of those effects can be deduced from the 3D displacement field of the modes. Because of the Poisson effect and the rotation of the sections , the displacement field on the nodes at other locations from the neutral line is not purely transverse for a bending eigenvector ${{{\varvec{\phi }}}}_p$ and not zero for a NB eigenvector ${{{\varvec{\phi }}}}_s$. This explains the losses of orthogonality observed above, as well as the loss of condition (35c), since the master part of ${{{\varvec{\phi }}}}_s$ is not zero: ${{{\varvec{\phi }}}}_s^\text {s}\ne {{{\varvec{0}}}}$.

6.2 Application to a clamped circular plate

In order to extend the results obtained on the beam test examples, the case of a clamped circular plate is here investigated. The selected plate has a radius $R=0.3~\text {m}$, a thickness $h=0.005\,\text {m}$, and the material properties are: density $\rho =7800~\text {kg/m}^3$, Young’s modulus $E=210~\text {GPa}$ and Poisson ratio $v=0.3$. As for the beam cases, a coarse mesh is chosen so as to compute all the modes and apply the different proposed methodologies. Consequently, 540 HEX20 elements on the face and 2 HEX20 elements in the thickness were used, with a total of 1931 nodes and 4928 degrees of freedom.

The convergence study and appearance of thickness modes are investigated for the fundamental axisymmetric bending mode of the clamped plate, as well as the first asymmetric (1,0) mode, having one nodal line and no nodal circle. The case of the first axisymmetric mode is awaited to share the same complexity as the beam case for symmetry reasons, but the asymmetric mode might be more difficult to achieve convergence.

Fig. 13

a Normalised modal correction factor for the clamped circular plate, as a function of the normalised mode number (normalization by the number of dofs), for the first asymmetric (1,0) mode of the plate. b The correction factors are now sorted by decreasing values, and two cases are shown : the case of the first asymmetic mode, corresponding to sorting (a), and the case of the first axisymmetric mode, showing a faster convergence. Grey points are negligible modes in terms of coupling, magenta points are the important in-plane coupled modes while blue points are the important thickness modes

Figure 13a shows the behaviour of the normalised modal correction factor $2(\alpha _{pp}^n/\omega _n)^2/\beta _{ppp}^p$ used in the previous sections, where p refers to the master mode (either $p=1$ for the first axisymmetric mode, or $p=2$ for the first asymmetric) and $n\,\in \, \{1,N\}$ with N the number of dofs. In Fig. 13a only the case of the first asymmetric mode is shown for the sake of brevity (thus $p=2$), but for $p=1$ the trend was very similar. As for the beam, a strong coupling with very high frequency modes is also observed. Investigating more precisely which modes are involved in the couplings, it is found that the ones having the most important correction factor are once again thickness modes.

Table 5 shows the deformed shape of the first nine modes, sorted according to their correction factor, which are thus the most important in the coupling with the bending (1,0) mode. Two purely in-plane modes are found, in position 5 and 8, and all others are thickness modes. The deformed shapes can be compared to that obtained for the beam and shown in Table 2. Indeed, the first thickness mode having the most important correction factor shows a similar geometry for both structures. Strong similarities are also observed between the second mode of the beam and mode (c) in Table 5, and the ninth mode in each case.

Figure 13b shows the normalized correction factor now sorted by order of decreasing values, and for the two cases of the axisymmetric fundamental mode and first asymmetric mode. It shows in particular that the convergence on the correct cubic coefficient is more rapidly achieved for the axisymmetric mode, where less than 10% of the modes are needed. On the other hand, the convergence is more difficult for the first asymmetric mode. Concerning the coupling with high-frequency modes and thickness modes for these two first bending modes, it is interesting to note that the subset of coupled modes is almost exactly the same in the two cases, showing in particular that the coupling with the thickness modes is not very dependent on the selected bending mode. Indeed, more than 90% of the coupled modes are the same for the two cases investigated.

Table 5 Mode shapes of the 9 most relevant modes coupled with the first flexural asymmetric (1,0) mode. Only two of them are in-plane modes: (e) and (h), while all the others are thickness modes. (a2) is a side view of the top view (a1) of the first thickness mode, in order to show the strong dependence on thickness deformation

Table 6 Corrected cubic coefficient $\varGamma _{ijk}^p$, for two flexural modes, i.e. ${i,j,k,p}\in [2,4]$, where 2 refers to the first (1,0) asymmetric mode while 4 refers to the second (2,0) asymmetric mode; and for three different methods : the modified StEP, the static condensation where all the coupled modes are statically condensed, and the Modal derivative where the added modal derivative is then statically condensed to the master mode

Table 6 gathers the numerical results for the corrected cubic coefficient $\varGamma _{ijk}^p$ defined in Eq. (51), with $p=2$ for the first asymmetric (1,0) mode and $p=4$ for the (2,0) asymmetric mode (the first bending modes being sorted by increasing frequencies, $p=1$ is the fundamental axisymmetric, $p=2,3$ for the two configurations of the (1,0) asymmetric mode and $p=4,5$ for the two configurations of the (2,0) mode with two nodal lines). This choice has been guided by the fact that these two asymmetric modes are coupled and thus shows important cubic coefficients that are needed if one wants to derive a reduced-order model. The three methods presented in the previous sections: M-STEP, static condensation of all the linear modes and of the modal derivative, give the same results, showing the convergence of the methods also in this case. Only $\varGamma ^2_{444}$ shows a different result for the M-STEP, however the value is very small as compared to the other ones so that this coefficient can be compared as negligible.

7 Conclusion

In this article, a nonlinear coupling of bending modes with thickness modes of very high frequency has been exhibited, due to geometrical nonlinearities in thin flat structures. This effect adds itself to the classical longitudinal/bending coupling and is the cause of a very slow convergence of a reduced order model (ROM) blindly built on a modal expansion of the nonlinear problem. It has been shown that if all eigenmodes are computed, it is possible to embed the effect of the non-bending modes into a master bending one, thus obtaining a reduced order model composed of only one nonlinear Duffing oscillator. This procedure can be done either by static condensation or by a normal form reduction, equivalent to the reduction on a single nonlinear mode. Finally, two alternative methods have been proposed to overcome the problem: the use of a static modal derivative or the direct computation of the cubic coefficients by an original method, the M-STEP, inspired by the standard STEP. Those methods have been successfully verified on dedicated examples, showing equivalent results.

Most of the results presented in this paper are restricted to the case of flat structures. Indeed, the specific shape of the equations of motion (see Eqs. (12), (13)) has been used to obtain exact equivalences between different methods. One can await that the obtained results should extend to shallow curved structures. However, for more generic shells with all the nonlinear couplings, most of the equivalences found here won’t probably hold anymore.

We focused on the case of a 3D model discretized by finite elements. In Sect. 3.1, we have shown on an example that some convergence problems might also appear for thin structures meshed with plate or shell elements, and having at least one long edge free. Our experience on thin ribbon have shown that the same kind of phenomenon appears when blindly using the STEP with the modal basis, and are again due to the loss of invariance of the modal eigenspaces. Indeed, high-frequency modes involving lateral deformations of the two free edges appeared. We also made computation on a circular plate with a free edge, and found circumferential modes appearing. Consequently, our finding is not restricted to 3D elements, and is completely linked to the use of the modal basis. Note that STEP calculations can also be realized with other input functions, and this has be done in this paper e.g. in Eq. (28). The complete investigation of the analogy between this contribution, focused on 3D elements, and the problems related to plate and shell finite elements, is postponed to a future work.

Appendices

Definition of cubic nonlinear terms

A particular feature of the manipulations of nonlinear coupling coefficients lies in the fact that they are related to monoms having symmetry relationships. For example, a cubic coefficient $\beta ^p_{ijj}$ is related to the monom $q_i q_j^2$, which is also the case of $\beta ^p_{jij}$ and $\beta ^p_{jji}$. Consequently many formulations use upper triangular forms for quadratic and cubic coefficients $\alpha ^p_{ij}$ (assuming $j \ge i$) and $\beta ^p_{ijk}$ (with $k \ge j \ge i$). However, direct calculations produce coefficients that have not this ordering property built-in. The formula given here allows one to get from to another formulation.

The corrected cubic coefficients in Eq. (15) can be derived by replacing the value of $p_s$ from Eq. (14) into the quadratic term of Eq. (12):

$$\begin{aligned} \sum _{i=1}^{N_B}\sum _{s=N_B+1}^{N} \alpha _{is}^r q_ip_s= \sum _{i=1}^{N_B}\sum _{j=1}^{N_B}\sum _{k=j}^{N_B}\sum _{s=N_B+1}^{N} -\dfrac{\alpha _{is}^r\alpha _{jk}^s}{\omega _s^2}q_iq_jq_k \end{aligned}$$

Manipulating the right hand side to make the sums over i, j, and k consistent with the ones of the cubic term of Eq. (12), i.e. having the sum over j from i to $N_B$, reads:

$$\begin{aligned} \begin{aligned}&\sum _{i=1}^{N_B}\sum _{j=1}^{N_B}\sum _{k=j}^{N_B} \sum _{s=N_B+1}^{N} -\dfrac{\alpha _{is}^r\alpha _{jk}^s}{\omega _s^2}q_iq_jq_k\\&\quad = \sum _{i=1}^{N_B}\sum _{j=i+1}^{N_B}\sum _{k=j+1}^{N_B}\sum _{s=N_B+1}^{N} \left( -\dfrac{{\alpha }^r_{is}\alpha ^s_{jk}}{\omega _s^2} -\dfrac{{\alpha }^r_{js}\alpha ^s_{ik}}{\omega _s^2} -\dfrac{{\alpha }^r_{ks}\alpha ^s_{ij}}{\omega _s^2}\right) q_iq_jq_k\\&\qquad +\sum _{i=1}^{N_B}\sum _{k=i+1}^{N_B} \sum _{s=N_B+1}^{N} \left( -\dfrac{{\alpha }^r_{is}\alpha ^s_{ik}}{\omega _s^2} -\dfrac{{\alpha }^r_{ks}\alpha ^s_{ii}}{\omega _s^2}\right) q_i^2q_k\\&\qquad +\sum _{i=1}^{N_B}\sum _{k=i+1}^{N_B} \sum _{s=N_B+1}^{N} \left( -\dfrac{{\alpha }^r_{is}\alpha ^s_{kk}}{\omega _s^2} -\dfrac{{\alpha }^r_{ks}\alpha ^s_{ik}}{\omega _s^2}\right) q_iq_k^2\\&\qquad +\sum _{i=1}^{N_B} \sum _{s=N_B+1}^{N} \left( -\dfrac{{\alpha }^r_{is}\alpha ^s_{ii}}{\omega _s^2} \right) q_i^3. \end{aligned} \end{aligned}$$

(50)

This term can be rewritten in a more compact form by defining the correction factor:

$$\begin{aligned} \mathcal {C}^{rs}_{ijk}= {\left\{ \begin{array}{ll} \dfrac{{\alpha }^r_{is}\alpha ^s_{jk}}{\omega _s^2}+ \dfrac{{\alpha }^r_{js}\alpha ^s_{ik}}{\omega _s^2}+ \dfrac{{\alpha }^r_{ks}\alpha ^s_{ij}}{\omega _s^2}\qquad i<j<k\\ \dfrac{{\alpha }^r_{is}\alpha ^s_{ik}}{\omega _s^2}+ \dfrac{{\alpha }^r_{ks}\alpha ^s_{ii}}{\omega _s^2}i=j<k\\ \dfrac{{\alpha }^r_{is}\alpha ^s_{kk}}{\omega _s^2}+ \dfrac{{\alpha }^r_{ks}\alpha ^s_{ik}}{\omega _s^2}i<j=k\\ \dfrac{{\alpha }^r_{is}\alpha ^s_{ii}}{\omega _s^2}i=j=k \end{array}\right. } \end{aligned}$$

leading to:

$$\begin{aligned}&\sum _{i=1}^{N_B}\sum _{j=1}^{N_B}\sum _{k=j}^{N_B} \sum _{s=N_B+1}^{N} -\dfrac{\alpha _{is}^r\alpha _{jk}^s}{\omega _s^2} q_iq_jq_k\\&\quad = \sum _{i=1}^{N_B}\sum _{j=i}^{N_B}\sum _{k=j}^{N_B}\sum _{s=N_B+1}^{N} -\mathcal {C}^{rs}_{ijk}q_iq_jq_k \end{aligned}$$

now with summation indexes consistent with the ones of Eq. (12).

Finally, the quadratic and cubic nonlinear terms in Eq. (12) read:

$$\begin{aligned}&\sum _{i=1}^{N_B}\sum _{s=N_B+1}^{N} \alpha _{is}^r q_ip_s+ \sum _{i=1}^{N_B}\sum _{j=i}^{N_B}\sum _{k=j}^{N_B} \beta ^r_{ijk} q_iq_jq_k\\&\quad = \sum _{i=1}^{N_B}\sum _{j=i}^{N_B}\sum _{k=j}^{N_B} \left( \beta ^r_{ijk}-\sum _{s=N_B+1}^{N} \mathcal {C}^{rs}_{ijk} \right) q_iq_jq_k. \end{aligned}$$

and the corrected cubic coefficient in Eq. (15):

$$\begin{aligned} \varGamma ^r_{ijk}=\beta ^r_{ijk}-\sum _{s=N_B+1}^{N} \mathcal {C}^{rs}_{ijk} \end{aligned}$$

(51)

Expression of static modal derivatives in terms of quadratic coupling coefficients

Given the general equation of a system with quadratic and cubic nonlinearities in physical coordinates:

$$\begin{aligned} {{{\varvec{M}}}}\ddot{{{\varvec{x}}}}+{{{\varvec{C}}}}\dot{{{\varvec{x}}}}+{{{\varvec{K}}}}{{\varvec{x}}}+{{{\varvec{f}}}}_\text {nl}({{\varvec{x}}})={{{\varvec{f}}}}_\text {e}, \end{aligned}$$

(52)

the nonlinear force can be written in terms of nonlinear tensors as:

$$\begin{aligned} {{{\varvec{f}}}}_\text {nl}({{\varvec{x}}})={{{\varvec{A}}}}{{\varvec{xx}}}+{{{\varvec{B}}}}{{\varvec{xxx}}}, \end{aligned}$$

(53)

where we used the compact tensor notation also employed in [13, 30]. In order to explicit the notation, the products ${{{\varvec{A}}}}{{\varvec{xx}}}$ and ${{{\varvec{B}}}}{{\varvec{xxx}}}$ are here given with explicit indicial notation:

$$\begin{aligned} {{{\varvec{A}}}}{{\varvec{xx}}} =&\sum _{i=1}^N \sum _{j=1}^N {{{\varvec{A}}}}_{ij}{x}_{j}{x}_{i} \, ,\\ {{{\varvec{B}}}}{{\varvec{xxx}}}=&\sum _{i=1}^N \sum _{j=1}^N \sum _{k=1}^N {{{\varvec{B}}}}_{ijk}{x}_{i}{x}_{j}{x}_{k}. \end{aligned}$$

The most inner product defined above coincides with a matrix product performed on the last index of the tensors.

In order to derive Eq. (23), we first express the ith column of nonlinear stiffness matrix as:

$$\begin{aligned}&\left( \dfrac{\partial {{{\varvec{f}}}}_\text {nl}}{\partial {{{\varvec{x}}}}}\right) _i=\sum _{j=1}^N {{{\varvec{A}}}}_{ij}{x}_{j}+\sum _{j=1}^N {{{\varvec{A}}}}_{ji}{x}_{j}+ \sum _{j=1}^N \sum _{k=1}^N {{{\varvec{B}}}}_{ijk}{x}_{j}{x}_{k}\nonumber \\&\quad + \sum _{j=1}^N \sum _{k=1}^N {{{\varvec{B}}}}_{jik}{x}_{j}{x}_{k} + \sum _{j=1}^N \sum _{k=1}^N {{{\varvec{B}}}}_{jki}{x}_{j}{x}_{k}. \end{aligned}$$

(55)

Exploiting the symmetry of the tensors ${{{\varvec{A}}}}$ and ${{{\varvec{B}}}}$ which implies that ${{{\varvec{A}}}}_{ij}={{{\varvec{A}}}}_{ji}$ and similarly ${{{\varvec{B}}}}_{ijk}={{{\varvec{B}}}}_{jik}={{{\varvec{B}}}}_{jki}$, that stems from the fact that geometric nonlinear forces can be derived from a potential (see [23]), the nonlinear stiffness matrix can be written in compact form as:

$$\begin{aligned} \dfrac{\partial {{{\varvec{f}}}}_\text {nl}}{\partial {{{\varvec{x}}}}}= 2{{{\varvec{A}}}}{{\varvec{x}}}+3{{{\varvec{B}}}}{{\varvec{xx}}}. \end{aligned}$$

(56)

The nonlinear stiffness matrix evaluated along mode p reads:

$$\begin{aligned} \dfrac{\partial {{{\varvec{f}}}}_\text {nl}}{\partial {{{{\varvec{x}}}}}}(\varvec{\varPhi }_p q_p)= 2q_p {{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_p +3q_p^2 {{{\varvec{B}}}} \varvec{\varPhi }_p\varvec{\varPhi }_p \, , \end{aligned}$$

(57)

and its derivatives with respect to $q_p$ evaluated at $q_p=0$ reads:

$$\begin{aligned} \left. \left( \dfrac{\partial }{\partial q_p} \dfrac{\partial {{{\varvec{f}}}}_\text {nl}}{\partial {{{{\varvec{x}}}}}}(\varvec{\varPhi }_p q_p) \right) \right| _{q_p=0}=2{{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_p. \end{aligned}$$

(58)

Hence the static modal derivatives defined in Eq. (19) leads to:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pr}=-2{{{\varvec{K}}}}^{-1}{{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_p{{{\varvec{\varPhi }}}}_r\, , \end{aligned}$$

(59)

and, when $p=r$, to:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pp}=-2{{{\varvec{K}}}}^{-1}{{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_p{{{\varvec{\varPhi }}}}_p. \end{aligned}$$

(60)

We now want to relate these expressions to the quadratic modal coupling coefficients $\alpha ^s_{pr}$ obtained from the STEP method, as well as expliciting directly how to compute the modal derivative from specific static evaluation of the internal force vector. When $p=r$, Eq. (11a) reads:

$$\begin{aligned} \alpha ^s_{pp}={{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\left( \dfrac{ {{\varvec{f}}}_\text {nl}(q_p{{\varvec{\varPhi }}}_p)+ {{\varvec{f}}}_\text {nl}(-q_p{{\varvec{\varPhi }}}_p)}{2q_p^2}\right) /m_s. \end{aligned}$$

(61)

When $p \ne r$, the STEP method needs to call two static evaluations with $q_p{{\varvec{\varPhi }}}_p$ and $q_s{{\varvec{\varPhi }}}_s$ so that the expression writes:

$$\begin{aligned} \alpha ^s_{pr}={{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\left( \dfrac{ {{\varvec{f}}}_\text {nl}(q_p{{\varvec{\varPhi }}}_p{+}q_r{{\varvec{\varPhi }}}_r){+} {{\varvec{f}}}_\text {nl}(-q_p{{\varvec{\varPhi }}}_p{-}q_r{{\varvec{\varPhi }}}_r){-}{{\varvec{f}}}_\text {nl}(q_p{{\varvec{\varPhi }}}_p){-}{{\varvec{f}}}_\text {nl}(-q_p{{\varvec{\varPhi }}}_p){-}{{\varvec{f}}}_\text {nl}(q_r{{\varvec{\varPhi }}}_r)-{{\varvec{f}}}_\text {nl}(-q_r{{\varvec{\varPhi }}}_r)}{2q_p q_r}\right) /m_s. \end{aligned}$$

(62)

Let us now assume that the eigenvectors are mass normalised so that $m_s=1,\;\forall s$, using Eq. (53) we can find the relation between $\alpha ^s_{pp}$ and ${{{\varvec{A}}}}$ as:

$$\begin{aligned} \alpha ^s_{pp}={{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\; {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \, . \end{aligned}$$

(63)

Similarly for the coefficients $\alpha ^s_{pr}$ with $p<r$:

$$\begin{aligned} \alpha ^s_{pr}=2{{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\; {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_r \, , \end{aligned}$$

(64)

where the factor 2 appears for symmetry reasons and is related to the usual problem of representing the polynomial monoms by counting them separately or not. Indeed, in the usual polynomial representation, the coefficients $\alpha ^s_{pr}$ with $p<r$ contains both values of ${{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\;{{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_r$ and ${{{{\varvec{\varPhi }}}}_s^{{\text {T}}}}\;{{{\varvec{A}}}}{{\varvec{\varPhi }}}_r{{\varvec{\varPhi }}}_p$ whereas the coefficients $\alpha ^s_{rp}$ are set to zero. We can now observe that each coefficient $\alpha ^s_{pr}$ can be seen as the row of a vector ${{\varvec{\alpha }}}_{pr}$ whose compact expression would be:

$$\begin{aligned} {{{\varvec{\alpha }}}}_{pr}=2{{{{\varvec{V}}}}^{{\text {T}}}}\; {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_r \, , \end{aligned}$$

(65)

where the full matrix of eigenvector ${{{\varvec{V}}}}$ has been introduced. In the same line, the expression for the vector ${{\varvec{\alpha }}}_{pp}$ when $p=r$ has the same shape but without the factor 2:

$$\begin{aligned} {{{\varvec{\alpha }}}}_{pp}={{{{\varvec{V}}}}^{{\text {T}}}}\; {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \, . \end{aligned}$$

(66)

Introducing Eqs. (65)–(66) respectively in Eqs. (59)–(60), and using the relationships of the quadratic coupling coefficients from the STEP method, Eqs. (62)–(61), one obtains the important formulas allowing one to compute directly the modal derivative from static FE calculations. The expression for ${{{\varvec{\theta }}}}_{pp}$ is given in the main text as Eq. (20) while the formula for ${{{\varvec{\theta }}}}_{pr}$ reads:

$$\begin{aligned} {{{\varvec{\theta }}}}_{pr}=-{{{\varvec{K}}}}^{-1}\left( \dfrac{ {{\varvec{f}}}_\text {nl}(\lambda ({{\varvec{\varPhi }}}_p+{{\varvec{\varPhi }}}_r))+ {{\varvec{f}}}_\text {nl}(-\lambda ({{\varvec{\varPhi }}}_p+{{\varvec{\varPhi }}}_r))-{{\varvec{f}}}_\text {nl}(\lambda {{\varvec{\varPhi }}}_p)-{{\varvec{f}}}_\text {nl}(-\lambda {{\varvec{\varPhi }}}_p)-{{\varvec{f}}}_\text {nl}(\lambda {{\varvec{\varPhi }}}_r)-{{\varvec{f}}}_\text {nl}(-\lambda {{\varvec{\varPhi }}}_r)}{2\lambda ^2}\right) \end{aligned}$$

(67)

We now want to relate more closely the modal derivative to the modal representation and derive an explicit expression showing that the SMD gathers the contributions of all coupled modes. For that purpose, one needs to express the usual orthonormality conditions shared by the matrix of eigenvectors ${{{\varvec{V}}}}$, assumed to be mass normalized:

$$\begin{aligned} {{{{\varvec{V}}}}^{{\text {T}}}}{{{\varvec{M}}}}{{{\varvec{V}}}}={{{\varvec{I}}}},\qquad \qquad {{{{\varvec{V}}}}^{{\text {T}}}}{{{\varvec{K}}}}{{{\varvec{V}}}}={{{\varvec{\varOmega }}}}^2, \end{aligned}$$

where ${{{\varvec{\varOmega }}}}^2$ is the diagonal matrix containing the squared eigenfrequencies.

Recalling eq. (59) we can now express the static modal derivatives in terms of ${{\varvec{\alpha }}}_{pr}$:

$$\begin{aligned} {{\varvec{\theta }}}_{pr}=-{{{\varvec{K}}}}^{-1}{{{\varvec{V}}}}^{-\text {T}}\;{{\varvec{\alpha }}}_{pr}, \end{aligned}$$

(68)

and for $p=r$:

$$\begin{aligned} {{\varvec{\theta }}}_{pp}=-2{{{\varvec{K}}}}^{-1}{{{\varvec{V}}}}^{-\text {T}}\;{{\varvec{\alpha }}}_{pp}. \end{aligned}$$

(69)

Moreover, using the orthogonality conditions:

$$\begin{aligned} {{{\varvec{K}}}}^{-1}{{{\varvec{V}}}}^{-\text {T}}={{{\varvec{V}}}}{{{\varvec{\varOmega }}}}^{-2}, \end{aligned}$$

(70)

we can express the static modal derivatives in the general case as:

$$\begin{aligned}&{{\varvec{\theta }}}_{pr}= -{{{\varvec{V}}}}{{{\varvec{\varOmega }}}}^{-2}\;{{\varvec{\alpha }}}_{pr},\end{aligned}$$

(71)

$$\begin{aligned}&{{\varvec{\theta }}}_{pp}= -2{{{\varvec{V}}}}{{{\varvec{\varOmega }}}}^{-2}\,{{\varvec{\alpha }}}_{pp}. \end{aligned}$$

(72)

In the particular case of a flat structure, the equations of motions have a simple structure recalled in Eqs. (12)–(13), so that $\alpha ^s_{pr}$ is nonzero only for the non-bending modes. Thanks to this simplification, one can express the static modal derivatives as a linear combination of the non-bending modes only as:

$$\begin{aligned}&{{\varvec{\theta }}}_{pr}= -\sum ^{N}_{s=N_B+1}{{\varvec{\varPhi }}}_s \dfrac{\alpha ^s_{pr}}{\omega ^2_s}, \end{aligned}$$

(73)

$$\begin{aligned}&{{\varvec{\theta }}}_{pp}= -\sum ^{N}_{s=N_B+1}2\,{{\varvec{\varPhi }}}_s \dfrac{\alpha ^s_{pp}}{\omega ^2_s}. \end{aligned}$$

(74)

Fig. 14

Analytical solution of a beam in pure bending. a Stress state $\sigma _{xx}$ in a cross section $\mathcal {S}$; b, c reference (in grey) and deformed (in blue) configurations of the beam. The Poisson effect has been magnified with a factor 20 in the cross section view regard to the prescribed displacement amplitude for. (Color figure online)

Corrected cubic coefficient obtained from static modal derivatives

This appendix aims at demonstrating that coefficient $\tilde{\varGamma }_{ppp}^p$ introduced in Eq. (27), i.e. by using an imposed displacement composed of the master mode plus a quadratic part containing the SMD, is exactly equal to that given in Eq. (16), and obtained thanks to the static condensation. For that purpose, let us recall that the imposed displacement in the first case reads:

$$\begin{aligned} {{\varvec{x}}}(q_p)=q_p\;{{\varvec{\varPhi }}}_p+\frac{1}{2}\;q_p^2\;{{\varvec{\theta }}}_{pp}. \end{aligned}$$

(75)

The cubic coefficient $\tilde{\varGamma }_{ppp}^p$ can be computed by using the general formula from the STEP method, Eq. (11a), by replacing the imposed displacement by the one given in Eq. (75). Consequently, one arrives at:

$$\begin{aligned} \tilde{\varGamma }_{ppp}^p={{{\varvec{\varPhi }}}_p^{{\text {T}}}}\left( {{{\varvec{f}}}}_\text {nl}(\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})-{{{\varvec{f}}}}_\text {nl}(-\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})\right) /2\lambda ^3. \end{aligned}$$

(76)

Using the explicit expression of the nonlinear forces from Eq. (53):

$$\begin{aligned}&{{{\varvec{f}}}}_\text {nl}(\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})= \lambda ^2{{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p\nonumber \\&\quad + \lambda ^3\left( \frac{1}{2} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\theta }}}_{pp}+\frac{1}{2} {{{\varvec{A}}}}{{\varvec{\theta }}}_{pp}{{\varvec{\varPhi }}}_p+ {{{\varvec{B}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \right) \nonumber \\&\quad +\mathcal {O}(\lambda ^4) \end{aligned}$$

(77)

$$\begin{aligned}&{{{\varvec{f}}}}_\text {nl}(-\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})= \lambda ^2{{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p\nonumber \\&\quad - \lambda ^3\left( \frac{1}{2} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\theta }}}_{pp}+\frac{1}{2} {{{\varvec{A}}}}{{\varvec{\theta }}}_{pp}{{\varvec{\varPhi }}}_p+ {{{\varvec{B}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \right) \nonumber \\&\quad +\mathcal {O}(\lambda ^4). \end{aligned}$$

(78)

Thus the difference between the above nonlinear forces reads:

$$\begin{aligned}&{{{\varvec{f}}}}_\text {nl}(\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})-{{{\varvec{f}}}}_\text {nl}(-\lambda {{\varvec{\varPhi }}}_p+\frac{1}{2}\lambda ^2{{\varvec{\theta }}}_{pp})\nonumber \\&\quad =2\lambda ^3\left( \frac{1}{2} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\theta }}}_{pp}+\frac{1}{2} {{{\varvec{A}}}}{{\varvec{\theta }}}_{pp}{{\varvec{\varPhi }}}_p+ {{{\varvec{B}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \right) \;+\;\mathcal {O}(\lambda ^5) \end{aligned}$$

(79)

and the coefficient $\tilde{\varGamma }_{ppp}^p$, neglecting high order terms, reads:

$$\begin{aligned} \tilde{\varGamma }_{ppp}^p=\frac{1}{2}{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\theta }}}_{pp}+ \frac{1}{2}{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\theta }}}_{pp}{{\varvec{\varPhi }}}_p+ {{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{B}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \end{aligned}$$

(80)

Following a similar argument as the one of Eq. (61) it is possible to show that the last term on the right hand side of Eq. (80) is the uncorrected cubic coefficient:

$$\begin{aligned} \beta ^p_{ppp}={{{\varvec{\varPhi }}}_p^{{\text {T}}}}{{{\varvec{B}}}}{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p{{\varvec{\varPhi }}}_p \end{aligned}$$

(81)

As regards to the first two terms on the right hand side of Eq. (80), using the link between SMD and non-bending modes of Eq. (23) valid in case of a flat structure, they can be written as:

$$\begin{aligned}&\frac{1}{2}{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{\varvec{\theta }}}_{pp} = -\frac{1}{2}\sum ^{N}_{s=N_B+1} {{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{{\varvec{\varPhi }}}}_s \;\dfrac{2\alpha ^s_{pp}}{\omega ^2_s} \end{aligned}$$

(82)

$$\begin{aligned}&\frac{1}{2}{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\theta }}}_{pp}{{\varvec{\varPhi }}}_p = -\frac{1}{2}\sum ^{N}_{s=N_B+1} {{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_s{{\varvec{\varPhi }}}_p \;\dfrac{2\alpha ^s_{pp}}{\omega ^2_s} \end{aligned}$$

(83)

where the index s spans over all the non-bending modes. Using the symmetry of the tensor ${{{\varvec{A}}}}$ once again:

$$\begin{aligned} {{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{{\varvec{\varPhi }}}}_s{{\varvec{\varPhi }}}_p+{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{{\varvec{\varPhi }}}}_s=2\;{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{{\varvec{\varPhi }}}}_s \end{aligned}$$

(84)

and in light of Eq. (64):

$$\begin{aligned} 2\;{{{\varvec{\varPhi }}}_p^{{\text {T}}}} {{{\varvec{A}}}}{{\varvec{\varPhi }}}_p{{{\varvec{\varPhi }}}}_s = \alpha ^p_{ps} \end{aligned}$$

(85)

because $p>s$ for each non-bending mode.

It is now possible to express the corrected cubic coefficient as:

$$\begin{aligned} \tilde{\varGamma }_{ppp}^p=\beta ^p_{ppp}-\sum ^{N}_{s=N_B+1} \alpha ^p_{ps} \;\dfrac{\alpha ^s_{pp}}{\omega ^2_s}=\varGamma _{ppp}^p \end{aligned}$$

(86)

and recover the same expression of the corrected coefficient given in Eq. (16) obtained with static condensation of all non-bending modes.

It is worth mentioning that the demonstration could have been done with another pathway. Indeed, the displacement introduced in Eq. (75) follows the general strategy proposed in [13, 30], consisting of using SMD to build a quadratic manifold approach in order to define a reduced-order models thanks to a nonlinear mapping between physical and reduced coordinates. However, the SMD can also be used more simply, by considering ${{\varvec{\theta }}}_{pp}$ as an enrichment of the modal basis, composed here of the single master eigenvector ${{\varvec{\varPhi }}}_p$. One could then derive a two-dofs reduced-order model, by projecting the general equations of motion onto these two vectors. Then, since the motions associated to the SMD are linked to NB modes having high frequencies, one can neglect their inertia and proceed to the static condensation of the part coming from the static modal derivative. By doing so, one would show again that the corrected cubic coefficient is still equal to $\varGamma _{ppp}^p$. Hence the static condensation of the SMD (a single vector) is thus strictly equivalent, in our simplified case of a flat structure (and considering only one master mode), to the static condensation of all coupled NB modes (including thickness modes).

Analytical solution for the pure bending of a beam

We consider the linear elastic solution of a beam of rectangular cross section under pure bending (see Fig. 14). The beam is thus subjected to a uniform bending moment ${{{\varvec{M}}}}_b=M{{{\varvec{e}}}}_z$. The material of the beam is assumed linear elastic with Young’s modulus E and Poisson’s ratio $\nu$. The orthonormal frame $({{{\varvec{e}}}}_x,{{{\varvec{e}}}}_y,{{{\varvec{e}}}}_z)$ is used, with ${{{\varvec{e}}}}_x$ colinear to the middle axis of the beam, ${{{\varvec{e}}}}_y$ the direction of bending and ${{{\varvec{e}}}}_z={{{\varvec{e}}}}_x\wedge {{{\varvec{e}}}}_y$ with $\wedge$ the vector product. The local equilibrium of the beam is exactly verified by the following axial stress state, linear through the thickness of the beam:

$$\begin{aligned} {\text {div}}\sigma ={{{\varvec{0}}}}\quad \Rightarrow \quad \sigma =\begin{pmatrix} \sigma _{xx} &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{pmatrix} ,\quad \sigma _{xx}=-\alpha y, \end{aligned}$$

(87)

where ${\text {div}}$ is the divergence operator, $\sigma$ is the stress tensor and $\alpha$ is a constant. The bending moment writes $M=\alpha I$ where I is the second moment of inertia of the beam. The linear strain tensor $\varepsilon$ verifies:

$$\begin{aligned}&{\varepsilon }=\frac{1}{2}\left( {{\varvec{\nabla }}}{{{\varvec{U}}}} + {{{\varvec{\nabla }}}^{{\text {T}}}}{{{\varvec{U}}}}\right) = \frac{1+\nu }{E}{\sigma }-\frac{\nu }{E}{\text {tr}}\sigma {\textsf {I}}_3\nonumber \\&\quad \Rightarrow \quad {\epsilon }=-\alpha \begin{pmatrix} y /E &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad -\nu y/E &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad -\nu y/E \end{pmatrix} \end{aligned}$$

(88)

As exposed in [31], the following displacement field verifies exactly the above equations:

$$\begin{aligned} {{{\varvec{U}}}}=\alpha xy\,{{{\varvec{e}}}}_x-\frac{\alpha }{2}\left[ x^2+\nu (y^2-z^2)\right] \,{{{\varvec{e}}}}_y-\nu \alpha y z \,{{{\varvec{e}}}}_z. \end{aligned}$$

(89)

This displacement is composed of three parts:

• the transverse (${{{\varvec{e}}}}_y$) component, proportional to $x^2$, which is the standard transverse displacement of the neutral line due to the bending. This term is the one directly computed in a beam theory;

• the axial (${{{\varvec{e}}}}_x$) component, which is the 3D linearized rotation of the cross section around vector ${{{\varvec{e}}}}_z$, also due to bending. In a beam theory, it is the consequence of the Euler-Bernoulli kinematics;

• two additional terms in the transverse (${{{\varvec{e}}}}_y$) and lateral (${{{\varvec{e}}}}_z$) directions, proportional to $\nu$ and thus directly linked to the Poisson effect. These terms are responsible of the distortion of the cross section, as seen in Fig. 14c.

To see the effect of the geometrical nonlinearities, we use the displacement field (89) to compute the nonlinear Green-Lagrange strain tensor. We obtain:

$$\begin{aligned} {\gamma }&=\frac{1}{2}\left( {{\varvec{\nabla }}}{{{\varvec{U}}}} + {{{\varvec{\nabla }}}^{{\text {T}}}}{{{\varvec{U}}}} + {{{\varvec{\nabla }}}^{{\text {T}}}}{{{\varvec{U}}}}{{\varvec{\nabla }}}{{{\varvec{U}}}}\right) \end{aligned}$$

(90)

$$\begin{aligned}&= {\varepsilon } + \underbrace{\frac{\alpha ^2}{2} \begin{pmatrix} x^2 &{}\quad 0 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{pmatrix} }_{{\gamma }_1}+\underbrace{\frac{\alpha ^2}{2} \begin{pmatrix} y^2 &{}\quad xy &{}\quad 0 \\ xy &{}\quad x^2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{pmatrix} }_{{\gamma }_2}\nonumber \\&\quad +\underbrace{\frac{\nu \alpha ^2}{2} \begin{pmatrix} 0 &{}\quad xy &{}\quad -xz \\ xy &{}\quad \nu (y^2+z^2) &{}\quad 0 \\ - xz &{}\quad 0 &{}\quad \nu (y^2+z^2) \end{pmatrix} }_{{\gamma }_3} \end{aligned}$$

(91)

The above equation shows that, in addition to the linear part ${\varepsilon }$, the nonlinear part has three type of components:

• the first nonlinear term ${\gamma }_1$ is purely axial, with the $\gamma _{xx}$ term $\alpha ^2x^2/2$ being the leading term responsible of the standard axial/bending coupling due to the geometrical nonlinearities. It is the term predicted by the von Kármán theory: if $v(x)=\alpha /2x^2$ denotes the transverse displacement of the neutral fiber, the nonlinear terms added by the von Kárman theory in $\gamma _{xx}$ is $v'(x)^2/2=\alpha ^2x^2/2$ (see [8]);

• the second nonlinear term ${\gamma }_2$ comes from purely 3D effects independent of the Poisson effect, that add (i) a stretch in the axial ${{{\varvec{e}}}}_x$ direction, proportional to $y^2$; (ii) a positive and homogeneous stretch in the transverse ${{{\varvec{e}}}}_y$ direction (proportional to ${{{\varvec{x}}}}^2$); (iii) a transverse shear;

• the third nonlinear term ${\gamma }_3$ gathers the effects of the 3D Poisson effect. It involves stretching in both the transverse ${{{\varvec{e}}}}_y$ and lateral ${{{\varvec{e}}}}_z$ direction as well as shear.

Even if the present results are strictly valid for a beam in pure bending, they can be extended and applied to understand qualitatively any bending state.