Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

Detecting internal resonances during model reduction

Evangelia Nicolaidou

Evangelia Nicolaidou

Department of Mechanical Engineering, University of Bristol, Bristol BS8 1TR, UK

[email protected]

Google Scholar

Find this author on PubMed

,
Thomas L. Hill

Thomas L. Hill

Department of Mechanical Engineering, University of Bristol, Bristol BS8 1TR, UK

Google Scholar

Find this author on PubMed

and
Simon A. Neild

Simon A. Neild

Department of Mechanical Engineering, University of Bristol, Bristol BS8 1TR, UK

Google Scholar

Find this author on PubMed

    Abstract

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

    1. Introduction

    Engineering structures vibrating at large amplitudes experience geometric nonlinearity, which couples the underlying linear normal modes of the model via nonlinear stiffness functions in the equations of motion. These couplings foster energy exchange between modes and can give rise to complex dynamic behaviours, which linear systems cannot exhibit and hence cannot be studied using methods from the linear vibration theory.

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

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

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

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

    In the second family of indirect methods, the dataset used for computing the nonlinear stiffness coefficients is obtained by applying a set of static forces to the structure, which are proportional to the reduced modeshapes, and extracting the corresponding displacement—this technique is referred to as the applied loads procedure or implicit condensation and expansion (ICE) [2932]. This approach relies on static condensation to account for the effect of membrane stretching and thus removes the need for including any high-frequency in-plane modes in the reduction basis as independent DOFs. Instead, these are assumed to be quasi-statically coupled to the reduced modes, and their effect can be captured implicitly. The ICE method has been successfully applied for the reduction of a range of flat structures, see, for example, refs. [20,33]; however, it has traditionally suffered from issues related to the lack of invariance of the reduced subspace. A main drawback of the ICE method is the fact that the computed ROM parameters vary with respect to the scaling of the forces used to obtain the static solution dataset, resulting in ROMs that lack robustness [29]. In addition, the ICE method assumes that the inertia of the condensed modes is negligible; this assumption breaks down when considering structures that can undergo large in-plane displacements, such as cantilever beams, and can lead to significantly inaccurate results. Recent works have addressed these limitations. Specifically, it has been shown that the quasi-static coupling between the reduced and condensed modes generates higher-order nonlinear terms in the reduced dynamics. When these are taken into account, the ROM parameter estimation procedure becomes robust with respect to the scaling of the static solution dataset [33,34]. In addition, it has been shown how the effect of the inertia of the condensed modes can be accounted for in the reduced dynamics; this additional treatment is referred to as inertial compensation (ICE(-IC)) [35]. This introduces some additional acceleration- and velocity-dependent terms in the reduced equations of motion, whose coefficients can be easily computed using the existing static solution dataset. Similar expressions in the reduced dynamics are seen in the direct model order reduction techniques such as the quadratic manifold with modal derivatives [36,37], as well as the more general concept of an invariant manifold based on the theory of normal forms [3841].

    Nevertheless, a major challenge associated with the ICE method, as well as reduced-order modelling techniques in general, remains: selecting the reduction basis without a priori knowledge of the full-order dynamics. Retaining redundant modes will lead to computationally suboptimal ROMs, while omitting dynamically significant modes will lead to inaccurate results, and important features such as internal resonances may not be captured. In other words, the ICE method relies on a slow/fast decomposition and is unable to capture any internal resonances between the reduced and condensed modes [42]. In this study, we aim to address this limitation and propose a method that can be used to predict the existence of internal resonances in conservative systems and thus guide the reduction basis selection process, without the need for full-order simulations. Specifically, we represent each condensed coordinate as the superposition of two components: one that is statically coupled to the reduced coordinates and another that is dynamically independent of them: the latter may be considered as the error associated with the static condensation of the mode in question. By using this framework, we show how these errors may be approximated during model reduction in a computationally efficient manner. This may serve as a tool for predicting internal resonances between the reduced and condensed coordinates, or, equivalently, for verifying the accuracy of ROMs by ensuring that static condensation approximation is sufficiently accurate for all operating conditions of interest.

    This article is structured as follows. In §2, we explore the nature of quasi-static and dynamic modal coupling in geometrically nonlinear systems using a simple oscillator as a motivating example. In §3, we show how the error associated with static condensation may be approximated during model order reduction, which can be used to predict the existence of internal resonances. In §§4 and 5, we demonstrate the proposed method using the simple oscillator and an FE model of a clamped–clamped beam, respectively. Finally, conclusions are presented in §6.

    2. Motivation

    In this section, we explore the nature of modal coupling in general conservative, geometrically nonlinear systems. Specifically, we investigate the quasi-static coupling approximation, which is often used in reduced-order modelling frameworks, and its applicability in different scenarios. To this end, we consider a discrete 4-DOF system composed of two point masses, m, as a motivating example. The masses are free to move in the x-y plane and are connected to a fixed frame and to each other through a set of linearly elastic springs, with stiffness kii{1,2,3,4,5} and unstretched length ℓ0, as shown in figure 1. At the equilibrium position, all springs are undeformed and oriented either horizontally or vertically. This system may be considered an extension to the single-mass, 2-DOF oscillator previously studied in refs. [34,43,44].

    Figure 1.

    Figure 1. Schematic diagram of the 4-DOF, two-mass oscillator used as a motivating example, shown at equilibrium with the springs unstretched and of length ℓ0. (Online version in colour.)

    Its equations of motion can be expressed in the form

    Mx¨+Kx+fx(x)=Fx, 2.1
    where x is the vector of physical displacements, M and K are the linear mass and stiffness matrices, respectively, and Fx and fx are the vectors of external and nonlinear restoring forces, respectively. These are given by
    x=(y1y2x1x2),M=[m0000m0000m0000m],K=[k40000k50000k1+k2k200k2k2+k3], 2.2a
    fx(x)=(k40(0+y1)d4k1y1(0d1)d1k2(y1y2)(0d2)d2+k40k50(0+y2)d5k3y2(0d3)d3k2(y2y1)(0d2)d2+k50k4x1(0d4)d4k10(0+x1)d1k20(x1x20)d2+(k1k2)0k5x2(0d5)d5k30(x20)d3k20(0+x2x1)d2+(k2k3)0), 2.2b
    d1=y12+(0+x1)2,d2=(y1y2)2+(0x1+x2)2,d3=y22+(0x2)2,d4=(0+y1)2+x12,d5=(0+y2)2+x22, 2.2c
    where di(x) i{1,2,3,4,5} are the lengths of the springs in the deformed configuration. The derivation of these can be found in appendix A. Here, we consider the dynamics of the system in its linear modal space, where the equations of motion are linearly uncoupled. The mode shape (ϕn) and natural frequency (ωn) of the nth mode of the underlying linear system are computed by solving the eigenproblem (Kωn2M)ϕn=0. After applying the transform x=Φq and pre-multiplying by ΦT, equation (2.1) can be rewritten as follows:
    q¨+Λq+f(q)=F, 2.3
    where Φ is the matrix containing the mass-normalized mode shapes1 in its columns, such that ΦTMΦ=I, Λ=ΦTKΦ is the diagonal matrix containing the corresponding squares of the natural frequencies along its leading diagonal, and F=ΦTFx and f(q)=ΦTfx(Φq) are, respectively, the vectors of external and nonlinear restoring forces in the modal space.

    For this motivating example, the physical parameters of the system are set to the following values: m = 1 kg, ℓ0 = 1 m, k1 = 1000 N m−1, k2 = 10 N m−1, k3 = 1000 N m−1, k4 = 1 N m−1 and k5 = 21.2 N m−1. Note that the horizontal grounding springs (k1 and k3) are significantly stiffer than the vertical grounding ones, which creates a dichotomy between the natural frequencies of the first two modes (ω1 = 1.0 rad s−1 and ω2 = 4.6 rad s−1), and those of the other two modes (ω3 = 31.6 rad s−1 and ω4 = 31.9 rad s−1). This system aims to emulate, in a highly simplified manner, the slow/fast dynamics that characterize the low-frequency transverse modes and high-frequency in-plane modes in plate- or beam-like structures.

    First, we consider the quasi-static behaviour of the system when only the first mode is forced directly. We numerically solve the static equations Λq+f(q)=F for q for a series of load cases where F1 ∈ [ − 3, + 3] and Fn=0n{2,3,4}, where Fn denotes the nth element in F. Figures 2a,b show the quasi-static modal response of the system against the static force applied to the first mode and against the corresponding response of the first mode, respectively. In reduced-order modelling methods such as the ICE(-IC) [30,35], this dataset, or more commonly a subset thereof, is used to approximate the functions describing the nonlinear stiffness of the reduced mode(s) and sometimes the quasi-static relationship between the condensed modes and the reduced mode(s). Specifically, the latter task is carried out only for a small set of high-frequency in-plane (or membrane) modes, for which the inertial forces are assumed to be small relative to the internal restoring forces. Here, we generalize this approach by treating all condensed modes equally, irrespective of their natural frequency or characteristics of their mode shapes. We denote the quasi-static relationship between the nth mode and the reduced mode using the function gn and approximate it as the Kth-order polynomial, i.e.

    gn(q1)=k=2KBk(n)q1k,n{2,3,4}, 2.4
    where the coefficients Bk(n) are identified via least-squares regression based on the static solution dataset, i.e. fitting to the curves shown in figure 2b.
    Figure 2.

    Figure 2. Quasi-static modal response of the 4-DOF oscillator, plotted against (a) the static force applied in the first mode and (b) the resulting response of the first mode. (Online version in colour.)

    Figure 3a shows the first backbone curve of the 4-DOF oscillator, computed according to equations (2.1) and (2.2) (with Fx = 0) using the MATLAB-based numerical continuation toolbox Continuation Core (COCO) [13]. The branch emerging near Ω=1.55rad s1 corresponds to a 1:3 internal resonance between the first and second modes. Figure 3b shows the time history of the modal response of the system (solid black lines), for three different nonlinear normal modes (NNMs) associated with the first backbone curve, which are represented by red dots in figure 3a and correspond to fundamental response frequencies of 1.01, 1.39 and 1.53 rad s−1, respectively. The quasi-static response of modes 2–4 is computed by evaluating the functions gn(q1) during the NNM motion: this is represented by dash-dotted lines in figure 3b. The difference between the dynamic and the quasi-static response of each mode, i.e. qn − gn(q1), is represented by dashed blue lines. This may be considered as the error arising from the quasi-static approximation/implicit condensation.

    Figure 3.

    Figure 3. (a) First backbone curve of the 4-DOF oscillator, shown in the projections of the maximum amplitudes of the first two modes, Q1 and Q2, against the response frequency, Ω. (b,c,d) Three nonlinear normal modes of the system, corresponding to fundamental response frequencies of 1.01, 1.39 and 1.53 rad s−1 and represented by red dots on the backbone curves. These are plotted as modal displacement, qn, against time, t (solid black lines). The dash-dotted red, purple and green lines show the quasi-static component of modes 2, 3 and 4, respectively, and the dashed blue lines show the corresponding error, qngn(q1),n={2,3,4}.(Online version in colour.)

    It can be seen that, as expected, the quasi-static approximation is sufficiently accurate when applied to the high-frequency modes, as qngn(q1),n{3,4},t,Ω. However, in the case of the second, low-frequency mode, the quasi-static approximation is initially moderately accurate near the first linear natural frequency, but becomes increasingly inaccurate as the system approaches internal resonance. Interestingly, for all NNM solutions, the error arising from this approximation appears to be a single-harmonic signal of frequency 3Ω. This suggests that the response of each condensed mode may be naturally decomposed into two parts: a component that is quasi-statically coupled to the reduced mode(s) and a component that is dynamically independent of the reduced mode(s). In the next section, we exploit this idea to show how the existence of dynamic interactions can be predicted and discuss how this is fundamentally relevant to reduced-order modelling.

    3. Predicting dynamic interactions

    (a) Derivation

    We start by considering the semi-discretized equations of motion of a conservative, linearly elastic, geometrically nonlinear structure, which are typically obtained through a finite element procedure. As mentioned earlier, these can be written as a series of N second-order ordinary differential equations in the form

    Mx¨+Kx+fx(x)=Fx 3.1a
    and
    q¨+Λq+f(q)=F 3.1b
    in the physical and mass-normalized linear modal spaces, respectively, where N is the number of DOFs of the FE model.

    Following ref. [35], but included here in the condensed form for completeness, we now separate the modal coordinates into three distinct classes. The first class, denoted by the subscript r (reduced), consists of a small set of modes, which, for a given set of operating conditions, contain the majority of the total energy in the system and are dynamically independent: these modes must form the basis of the ROM. The number of modes in this class, R ≪ N, dictates the lower limit to the number of DOFs that an accurate ROM must have. For a single backbone curve in the absence of any internal resonance, R = 1. The second class, denoted by the subscript s (static), is composed of S ≪ N modes, which may contain a substantial fraction of the energy of the full-order system, yet their response can be approximated as being quasi-statically coupled to the reduced modes. These modes need not be included as independent DOFs in the reduction basis, as their effects can be incorporated implicitly in the reduced dynamics [30,31,35]. Finally, the third class, denoted by the subscript u (unmodelled), contains the remaining U = (N − R − S) modes, which are very weakly coupled to the reduced modes, such that they always contain a negligible amount of energy under the operating conditions of interest. As such, it is assumed that these modes can be ignored during the reduction process with negligible loss of accuracy. By using this framework, the modal equations of motion of the FE model, equation (3.1b), can be rewritten as follows:2

    [q¨rq¨sq¨u]+[Λr000Λs000Λu][qrqsqu]+[fr(qr,qs,qu)fs(qr,qs,qu)fu(qr,qs,qu)]=[Fr00]. 3.2
    Ignoring the third group of weakly coupled modes, these can be approximated as follows:
    [r¨s¨]+[Λr00Λs][rs]+[f^r(r,s)f^s(r,s)]=[Fr0], 3.3
    where f^r(r,s):=fr(r,s,0) and f^s(r,s):=fs(r,s,0), such that qr ≈ r, qs ≈ s, qu ≈ u = 0 and xΦrr+Φss. Equivalently, the kinetic energy of the full-order system, T, can be approximated as follows:
    T(r,s)=12(r˙)Tr˙+12(s˙)Ts˙, 3.4
    while the potential energy function, V(r,s), is such that
    Vr=Λrr+f^r(r,s) 3.5a
    and
    Vs=Λss+f^s(r,s). 3.5b

    As discussed in ref. [35], when the modes in the second group, s, can be expressed as functions of the reduced modes, r, i.e. s = g(r), then equation (3.3) can be exactly reduced to

    r¨+(gr)Tgrr¨+(gr)T(2gr2r˙)r˙+Λrr+f~r(r)=Fr, 3.6
    where g(r) is an S × 1 vector of quasi-static coupling functions, ∂g/∂r is its S × R Jacobian matrix, ∂2g/∂r2 is its S × R × R second derivative tensor and f~r(r):=f^r(r,g). Reduced-order models based on equation (3.6) were found to produce remarkably accurate results for systems where a clear slow/fast dynamic behaviour can be observed, e.g. between the low-frequency transverse modes and the highly stiff in-plane modes in thin plates and slender beams.

    Here, we aim to broaden the scope of the implicit condensation approach and seek to quantify the error introduced by the static condensation. We define this error using the S × 1 time-dependent vector h(t), i.e.

    s=g(r)+h. 3.7
    Using equation (3.7), and noting that s˙=(g/r)r˙+h˙, the Lagrangian of the system, L, can be expressed as follows:
    L(r,h)=T(r,g+h)V(r,g+h)=12(r˙)Tr˙+12(r˙)T(gr)Tgrr˙+(r˙)T(gr)Th˙+12(h˙)Th˙V(r,g+h). 3.8
    From this, the equations of motion for r and h can be derived using the Euler–Lagrange equation. This leads to
    r¨+(gr)Tgrr¨+(gr)Th¨+(gr)T(2gr2r˙)r˙+Λrr+f^r(r,g+h)+(gr)T(Λs(g+h)+f^s(r,g+h))=Fr 3.9a
    and
    h¨+grr¨+(2gr2r˙)r˙+(Λs(g+h)+f^s(r,g+h))=0, 3.9b
    as shown in appendix B.

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

    Λrr+f^r(r,g+h)=Λrr+f^r(r,g)+B(r)h+O(h2) 3.10a
    and
    Λs(g+h)+f^s(r,g+h)=Λsg+f^s(r,g)+C(r)h+O(h2)=C(r)h+O(h2), 3.10b
    where
    B(r)=f^r(r,s)ss=g 3.11a
    and
    C(r)=Λs+f^s(r,s)ss=g. 3.11b
    Note that, by definition, Λsg+f^s(r,g)=0, as g(r) is computed based on the static response of the system, with a static force only applied to the reduced modes (i.e. Fs = 0). Substituting equation (3.10) into equation (3.9), and noting that B + (∂g/∂r)TC = 0 as shown in appendix C, leads to
    r¨+(gr)Tgrr¨+(gr)T(2gr2r˙)r˙+Λrr+f~r(r)=Fr(gr)Th¨+O(h2) 3.12a
    and
    h¨+C(r)h+O(h2)=p(r,r˙,r¨), 3.12b
    where
    p(r,r˙,r¨)=grr¨(2gr2r˙)r˙. 3.13

    As expected, it can be seen that, when h = 0, equation (3.12a) is equivalent to the reduced dynamics obtained for the perfectly statically coupled case (equation (3.6)). The additional h- and h¨-dependent terms that appear on the right-hand side of equation (3.12a) may be considered as a forcing arising due to the dynamic coupling between r and s. In addition, when the kinetic energy of the statically coupled modes is negligibly small (i.e. s˙Ts˙0), the g-dependent terms in equation (3.6) are negligible, and the reduced dynamics can be simply expressed as follows:

    r¨+Λrr+f~r(r)=Fr. 3.14
    As discussed in ref. [35], this more traditional form of the ROM equations of motion is suitable for structures in which in-plane displacements are limited.

    The remainder of this article demonstrates how the equations of motion for h (equation (3.12b)) may be used to efficiently predict the presence of dynamic coupling between the reduced modes (r) and the condensed modes (s). This allows features such as internal resonances to be predicted and ROMs to be validated without the need for full-order FE simulations.

    (b) Use in reduced-order modelling frameworks

    The linear properties of the ROM, Λr, Φr and Φs, can be computed directly using the linear mass and stiffness matrices of the FE model. The reduced nonlinear stiffness functions, f~r(r), and the quasi-static coupling functions, g(r), are approximated indirectly in a least-squares manner using a force–displacement dataset for a series of nonlinear static solutions extracted from the FE model, as detailed in refs. [29,35]. These are approximated as the Kth-order polynomials of the reduced coordinates. Once the ROM parameters are identified, the reduced backbone curves can be computed, e.g. using numerical continuation, based on either equation (3.6) or equation (3.14).

    Using equation (3.12b), additional insight can now be gained by simulating the error dynamics for each NNM of the ROM. Given that the static modal coupling is well captured through the functions g(r), then, in the absence or near the onset of a dynamic interaction, the error h is expected to be relatively small such that any nonlinear monomials of h become negligible, i.e. O(h2)0. As such, the linearized version of equation (3.12b) may be considered, which can be solved efficiently using, for example, the harmonic balance method [1]. To this end, the S × S matrix containing the linear coefficients of h must be approximated as functions of r. As with f~r(r) and g(r), these can be computed based on least-squares polynomial regression. In this case, the tangent stiffness matrix, Ktan, must be extracted for each nonlinear static load case, in addition to the vectors of applied forces and resulting displacement, which are needed for the standard ICE(-IC) method. The coefficients of each function in C(r) are then computed using the {r,ΦsTKtanΦs} dataset evaluated at each load case. The practical implications of this are considered later in §5, where the proposed techniques are demonstrated using the commercial FE software Abaqus.

    Once these functions are approximated, then the matrix C, as well as the vector on the right-hand side of equation (3.12b), p, can be evaluated during a periodic solution of the ROM. These can then be expressed as summations of sinusoidal components, i.e.

    Cij=a(ij)0+n=1Nh[a(ij)ncos(nωt)+b(ij)nsin(nωt)] 3.15a
    and
    pi=α0(i)+n=1Nh[αn(i)cos(nωt)+βn(i)sin(nωt)], 3.15b
    where Nh is the number of harmonics, and a(ij)n, αn(i)n[0,Nh]Z, and b(ij)n, βn(i) n[1,Nh]Z, are coefficients, which can be identified via a discrete Fourier transform for each element in C and p, during an NNM motion with fundamental response frequency ω. Similarly, each element in h can be expressed as a sum of its Fourier components, i.e.
    hi=A0(i)+n=1Nh[An(i)cos(nωt)+Bn(i)sin(nωt)] 3.16a
    and
    h¨i=n=1Nh(nω)2[An(i)cos(nωt)+Bn(i)sin(nωt)]. 3.16b
    Using equations (3.15) and (3.16), the terms in equation (3.12b) can be expanded, the coefficients of like harmonic terms on either side of the equation equated, and the harmonic amplitudes of h computed by solving the resulting set of simultaneous linear equations. This can be expressed as follows:
    ch=Υ1cp, 3.17
    where ch and cp are vectors containing the harmonic coefficients of h and p, respectively, and Υ is a matrix that can be algorithmically populated with the harmonic coefficients of C. Further details are given in appendix D. Finally, the computed amplitudes in ch are used to estimate the time history of the error corresponding to each condensed mode during a periodic motion of the ROM (equation (3.16a)).

    As discussed in §1, the main limitation of the ICE(-IC) method remains the fact that it relies on a slow/fast decomposition between the reduced and condensed modes. With the method proposed herein, this limitation is overcome, as the condensed basis can now be formed using any modes, irrespective of their natural frequency. This method can be used to monitor the error associated with the static condensation and thus identify if/when a condensed mode becomes resonant. This would suggest the quasi-static coupling assumption is no longer appropriate for the mode in question, and instead the mode must be included as an independent DOF in the reduction basis (Φr). Equivalently, this method can be used as a tool for validating ROMs, by ensuring that the component of the condensed modes that is dynamically independent of the reduced modes (i.e. h) remains sufficiently small for all operating conditions of interest.

    4. Application to the 4-DOF oscillator

    (a) Single-mode ROM results

    We now revisit the 4-DOF oscillator considered in §2 to demonstrate our proposed method. We first compute a seventh-order (K = 7), single-DOF ROM of the first mode using the ICE(-IC) method (equation (3.6)), while the remaining three modes are included in the statically condensed basis, i.e. Φr=[ϕ1], Φs=[ϕ2ϕ3ϕ4], and Φu is unpopulated as no modes are unmodelled. The dataset used to calibrate the ROM consists of a series of static solutions where the the static force applied to the first mode is equally spaced between −3 and +3. An example of the fitting procedure for the reduced nonlinear stiffness function (f~r(r)) and the quasi-static coupling functions (g(r)) is shown in figures 4a,b, respectively.

    Figure 4.

    Figure 4. Examples of the calibration procedure for the single-DOF ROM of the oscillator: (a) reduced nonlinear stiffness, (b) quasi-static coupling, and (c,d) tangent stiffness matrix components, all approximated as seventh-order polynomial functions of the reduced coordinate. (Online version in colour.)

    Figure 5 shows the backbone curve of the computed single-DOF ROM, in the projection of modal amplitudes against fundamental response frequency (blue lines). Note that modes 2–4 are not modelled directly, but their response is approximated using the quasi-static coupling functions. It can be seen that the primary response of the reduced mode, and that of the high-frequency condensed modes, can be accurately predicted by the ROM. As expected, however, the internal resonance near Ω=1.6rad s1 cannot be captured. This is due to the dynamic energy transfer between modes that takes place during an internal resonance, an effect that the static condensation approach is unable to capture.

    Figure 5.

    Figure 5. Backbone curve of the seventh-order, single-DOF ICE(-IC) ROM of the oscillator (blue lines). The backbone curve of the ROM when the static condensation error is taken into account, using our proposed method, is represented by red lines. The first backbone curve of the full-order system is also shown for reference (dashed grey lines).(Online version in colour.)

    (b) Internal resonance prediction

    By using equation (3.12b) and the harmonic balance method, as described in §3b and appendix D, the error associated with the static condensation of each mode can now be estimated. To achieve this, the elements of the tangent stiffness matrix corresponding to the condensed modes, i.e. Λtan,s=ΦsTKtanΦs, must be approximated as functions of the reduced coordinates.3 Similar to the approximation of the reduced nonlinear stiffness functions and the quasi-static coupling functions, these are computed in a least-squares manner, based on the same set of full-order nonlinear static solutions. Examples of this are shown in figures 4c,d.

    The improved prediction of the reduced backbone curves, which takes into account the approximated error arising from the static condensation (with Nh = 7), is now represented by red lines in figure 5. From this, it can be seen that, as the response frequency increases, the magnitude of h2 relative to g2 becomes increasingly large until Ω1.6 rad s−1, before it rapidly decreases again.4 This singularity-type behaviour is caused by the third harmonic component of h2, and it indicates that a dynamic interaction between the first and second modes exists, without the need for simulating the dynamics of both modes simultaneously. This suggests that, in this region, a single-mode ROM is no longer appropriate, and the second mode must be included in the reduction basis. Note that, when considering FE models using commercial packages, the backbone curves of the full-order model cannot be readily computed. As such, significant features of the system, such as internal resonances, can be often overlooked during model reduction. With our proposed method, the existence of such dynamic interactions can be predicted, without the need for expensive full-order simulations.

    (c) Two-mode ROM results

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

    Figure 6.

    Figure 6. First backbone curve of the seventh-order, two-DOF ICE(-IC) ROM of the oscillator (blue lines). The backbone curve of the ROM when the static condensation error is taken into account, using our proposed method, is represented by red lines. The first backbone curve of the full-order system is also shown for reference (dashed grey lines).(Online version in colour.)

    5. Application to a finite element model

    We now consider a geometrically nonlinear model of a clamped–clamped beam, constructed using the commercial FE software Abaqus, with the Abaqus2Matlab [45] toolbox used for pre- and post-processing. The beam model is identical to the one studied in ref. [35]. It has a length of 300 mm and a rectangular cross-section of 25 mm × 1 mm and is made of steel with a Young’s modulus of 205 GPa, mass density of 7800 kg m−3 and Poisson’s ratio 0.3. The mesh consists of 120 shear deformable, three-node quadratic beam elements (B32), resulting in 1434 DOFs.

    We compute a quintic5 single-DOF ROM of the first (bending) mode using the ICE method and include modes 3, 6, 72 and 129 in the condensed basis: these correspond to the second and third symmetric bending modes, and the first and second symmetric axial modes, respectively.6 The static solution dataset used to calibrate the ROM consists of four load cases where the static force applied to the first mode is F1 = { − 45, − 22.5, + 22.5, + 45}, resulting in a maximum static transverse deflection of 1.13 mm at the beam midspan. For this model, the computation time for each load case is about 25 seconds on a standard computer. As mentioned earlier, the fitting procedure for the reduced nonlinear stiffness functions, the quasi-static coupling functions and the tangent stiffness functions for the condensed modes7 is carried out in a least-squares manner.

    Figure 7a shows the backbone curve of the single-DOF ROM (top), as well as the corresponding normalized amplitude of the error in the condensed modes (bottom). As the backbone curves of the full-order FE model cannot readily be computed, and thus cannot be directly compared with those of the ROM, we assess the accuracy of the ROM by comparing a set of reduced NNMs with the dynamic response of the FE model when the corresponding initial conditions are applied for a wide range of response frequencies. The initial conditions are enforced in the form of initial applied modal forces, as proposed in ref. [47]. The periodic orbits of the ROM are compared with the FE response in the time domain, as shown in figure 7b. From this, it can be seen that the ICE method can provide accurate results for the reduction of the clamped–clamped beam model, even with a single-mode ROM: this agrees with observations seen in the literature [20,29,30].

    Figure 7.

    Figure 7. (a) Backbone curve of the quintic single-DOF ROM of the beam (top), and the corresponding error associated with the static condensation of each mode (bottom). (b) Comparison between the periodic responses predicted by the ROM and the responses of the FE model, for 10 different sets of initial conditions, which correspond to the black dots on the backbone curve. (Online version in colour.)

    Interestingly, while the standard ROM results, as shown in figure 7b, indicate good agreement with the full-order model, the results obtained from the error-monitoring treatment suggest that there is a strong dynamic interaction between the first and third modes of the beam near Ω=410rad s1, which the single-mode ROM is unable to capture. Note that, while the h components of modes 6, 72 and 129 are also large, the component of mode 3 is the largest and begins its rapid growth at a lower frequency than other modes. As such, mode 3 is treated as the candidate for a dynamic interaction. Specifically, this is associated with the amplitude of the fifth harmonic of the third mode. This result points to the existence of a 1:5 internal resonance between the first and third modes, a feature of flat clamped–clamped beams that has been widely observed in the literature, e.g. [19,20].

    To verify this observation, we compute a quintic two-DOF ICE ROM of the beam, with Φr=[ϕ1ϕ3] and Φs=[ϕ6ϕ72ϕ129]. The calibration dataset consists of 24 load cases, where the maximum force applied in each reduced mode is F1 = 45 and F3 = 360. The resulting backbone curve is shown in figure 8a. It can be seen that the ROM now exhibits a 1:5 internal resonance, while the error predicted using our proposed method remains relatively small. As mentioned earlier, the accuracy of the reduced internally resonant NNMs is verified by comparing them with the corresponding set of responses of the full-order model (figure 8b).

    Figure 8.

    Figure 8. (a) First backbone curve of the quintic two-DOF ROM of the beam (top) and the corresponding error associated with the static condensation of each mode (bottom). (b) Comparison between the periodic responses predicted by the ROM and the responses of the FE model, for three different sets of initial conditions that correspond to the black dots on the backbone curve. (Online version in colour.)

    It should be noted that, in addition to enabling the internal resonance to be modelled, the two-DOF ROM also leads to more accurate response predictions on the primary backbone curve. We quantify the accuracy of each reduced NNM using the periodicity metric ϵ, as defined in ref. [48], i.e.

    ϵ=xTx0x0, 5.1
    where x0 is the displacement vector of the FE model at t = 0, which is imposed as an initial condition based on the reduced NNM, and xT is the displacement vector of the FE model after one period. A smaller ϵ value indicates a response, which is closer to being periodic and thus a more accurate ROM. The computed periodicity values for different NNMs on the primary backbone curves, both for the single- and two-DOF ROMs, are shown in figure 9. The results suggest that the third mode of the FE model becomes increasingly important for response frequencies higher than ∼420 rad s−1. This is in qualitative agreement with the results shown in figure 7a and further confirms the validity of the error-approximation procedure.
    Figure 9.

    Figure 9. Plot of the periodicity error of the FE model for different sets of initial conditions given by the single-DOF (blue crosses) and two-DOF (red circles) ICE ROMs of the clamped–clamped beam.(Online version in colour.)

    6. Conclusion

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

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

    Data accessibility

    FE model submitted as electronic supplementary material.

    Authors' contributions

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

    Competing interests

    We declare we have no competing interests.

    Funding

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

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

    The potential energy of spring i is given by

    Vi=12ki(di0)2, A 1
    where di is its length in the deformed configuration, as defined in equation (2.2c). The Lagrangian of the system can be expressed as follows:
    L=T(x˙)V(x)=12m(y˙12+y˙22+x˙12+x˙22)12i=15ki(di0)2, A 2
    where T denotes the kinetic energy of the system. Applying the Euler–Lagrange equation,
    ddt(Lx˙)Lx=Fx, A 3
    and after some algebraic manipulation, leads to equations (2.1) and (2.2).

    Appendix B. Derivation of the reduced and error dynamics

    Using equation (3.8), the partial derivatives of L with respect to r˙, h˙, r and h can be written, respectively, as follows:

    L^r˙=r˙+(gr)Tgrr˙+(gr)Th˙ B 1a
    L^h˙=grr˙+h˙ B 1b
    L^r=(gr)T(2gr2r˙)r˙+(2gr2r˙)Th˙(Λrr+f^r(r,g+h))(gr)T(Λs(g+h)+f^s(r,g+h)) B 1c
    L^h=(Λs(g+h)+f^s(r,g+h)). B 1d
    Differentiating equations (B 1a) and (B 1b) with respect to time gives
    ddt(L^r˙)=r¨+(gr)Tgrr¨+2(gr)T(2gr2r˙)r˙+(2gr2r˙)Th˙+(gr)Th¨ B 2a
    and
    ddt(L^h˙)=grr¨+(2gr2r˙)r˙+h¨. B 2b
    The equations of motion for r and h can be derived using the Euler–Lagrange equation, i.e.
    ddt(Lr˙)Lr=Fr B 3a
    and
    ddt(Lh˙)Lh=0. B 3b
    Substituting equations (B 1) and (B 2) into equation (B 3) leads to equation (3.9).

    Appendix C. Relationship between B and C

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

    (Λss+f^s(r,s))s=g=0. C 1
    Taking the partial derivative of equation (C 1) with respect to r leads to
    f^s(r,s)rs=g+(Λs+f^s(r,s)s)s=ggr=0. C 2
    Substituting equation (3.5b), i.e. f^s(r,s)=V/sΛss, and equation (3.11b), i.e. C(r)=Λs+f^s/ss=g, into equation (C 2) leads to
    2Vrss=g+Cgr=0and(2Vsr)Ts=g+Cgr=0.} C 3
    Finally, substituting equation (3.5a), i.e. (V/r)=Λrr+f^r(r,s), and equation (3.11a), i.e. B(r)=f^r/ss=g, into equation (C 3), gives
    (f^r(r,s)s)Ts=g+Cgr=0andBT+Cgr=0,} C 4
    This is equivalent to B + (∂g/∂r)TC = 0, since C = CT.

    Appendix D. Estimation of the harmonic coefficients of h

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

    h¨i+j=1SCijhj=j=1Sc0(ij)+k=1Nhcos(kωt)[(kω)2Ak(i)+j=1Sck(ij)]+k=1Nhsin(kωt)[(kω)2Bk(i)+j=1Sdk(ij)], D 1
    where
    c0(ij)=a(ij)0A(j)0+12n=1Nh(a(ij)nA(j)n+a(ij)nB(j)n) D 2a
    ck(ij)=a(ij)0A(j)k+a(ij)kA(j)0+12n=1k1(a(ij)knA(j)nb(ij)knB(j)n)+12n=1Nhk(a(ij)n+kA(j)n+b(ij)n+kB(j)n)+12n=k+1Nh(a(ij)nkA(j)n+b(ij)nkB(j)n),k[1,Nh]Z D 2b
    dk(ij)=a(ij)0B(j)k+b(ij)kA(j)0+12n=1k1(a(ij)knB(j)n+b(ij)knA(j)n)+12n=1Nhk(a(ij)n+kB(j)n+b(ij)n+kA(j)n)+12n=k+1Nh(a(ij)nkB(j)nb(ij)nkA(j)n),k[1,Nh]Z. D 2c
    By using equations (D 1) and (3.15b), the expressions on either side of equation (3.12b) can be directly compared. The coefficients of like harmonic terms are equated and expressed in the form
    Υch=cp, D 3
    where the S(2Nh + 1) × 1 vectors ch and cp contain the harmonic amplitudes of h and p, respectively, i.e.
    ch(i)=[A0(i)A1(i)B1(i)A2(i)B2(i)ANh(i)BNh(i)] D 4a
    ch=[ch(1)ch(2)ch(S)]T D 4b
    cp(i)=[α0(i)α1(i)β1(i)α2(i)β2(i)αNh(i)βNh(i)] D 4c
    cp=[cp(1)cp(2)cp(S)]T. D 4d
    The S(2Nh + 1) × S(2Nh + 1) matrix Υ(a(ij)n,b(ij)n) is populated according to equation (12). Finally, the harmonic coefficients of h are computed by inverting equation (13), i.e. ch=Υ1cp.

    Footnotes

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

    2 Note that, by definition, only the reduced modes have static forces applied during the model reduction procedure, i.e. Fs = 0 and Fu = 0.

    3 Note that Λ:=Λtan(q=0).

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

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

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

    7 For this additional step, the tangent stiffness matrix of the full-order model must be evaluated at each static solution: this can be readily extracted from the FE package.

    Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5462613.

    Published by the Royal Society. All rights reserved.

    References