Nonlinear reduced-order modelling for limit-cycle oscillation analysis

A nonlinear model reduction based on eigenmode decomposition and projection for the prediction of sub- and supercritical limit-cycle oscillation is presented herein. The paper focuses on the derivation of the reduced-order model formulation to include expansion terms up to fifth order such that higher-order nonlinear behaviour of a physical system can be captured. The method is applied to a two degree-of-freedom pitch–plunge aerofoil structural model in unsteady incompressible flow. Structural stiffness nonlinearity is introduced as a fifth-order polynomial, while the aerodynamics follow linear theory. It is demonstrated that the reduced-order model is capable of accurately capturing sub- and supercritical limit-cycle oscillations arising both from initial disturbances and gust excitation. Furthermore, an analysis of the computational cost associated with constructing such reduced-order model and its applicability to more complex aeroelastic problems is given.


Introduction
Limit-cycle oscillation (LCO) is an inherently nonlinear phenomenon and a critical aspect of stability analysis in aircraft aeroelasticity. Limit-cycle oscillations can arise from concentrated structural nonlinearities which are inherent in aeroelastic systems such as freeplay at interconnection between different components including the hinge of a trailing-edge control surface or the connection between wing, pylons and engines [18]. Another type of concentrated structural nonlinearity is cubic stiffness which can typically represent a wing section undergoing torsional deformation. In fact, from an experimental perspective, concentrated structural nonlinearity in the form of polynomial stiffness (not limited to cubic) has been identified in experimental models through data curve fitting [24,26]. Sources of LCO have also been attributed to distributed structural nonlinearity, namely the geometric nonlinear stiffening of plate-like structures such as low-aspectratio wings [9]. Limit-cycle oscillation is also induced by aerodynamic nonlinearities arising from flow phenomena such as the separation at high incidence angles and unsteady shock oscillations in transonic flow [18].
In terms of LCO associated with concentrated structural nonlinearity, if the nonlinearity is purely hardening, supercritical LCO can occur which is characterised by a smooth growth in oscillation amplitude with respect to increasing flight speed above the flutter point. The system is stable below the flutter point. If the nonlinearity is softening however, the system can become unstable even below the flutter point. This is shown in [30] as one of the earliest studies of LCO arising from concentrated structural nonlinearity. The nonlinearities examined were namely cubic stiffness, freeplay and hysteresis. In the case of the cubic hardening spring, LCO was observed at flow velocities above the linear flutter speed, while softening springs in general had a destabilising effect on flutter. If, in addition, higher-order hardening nonlinearities exist, then subcritical LCO can arise which is characterised by stable large amplitude oscillations below the flutter point growing with increasing flight speed above the flutter point [9]. Subcritical LCO is highly undesirable due to the sudden large amplitude nature and the potential to manifest at velocities below the flutter point.
Limit-cycle oscillations have been studied extensively in the time domain. The amplitudes of LCO arising from cubic hardening stiffness nonlinearity in relation to the flow velocity have been studied explicitly using both numerical integration and asymptotic theory [17]. Monte Carlo simulations have been applied in the study of aerofoil LCO resulting from a subcritical bifurcation subject to parametric uncertainty in the cubic softening and quintic hardening stiffness coefficients of the nonlinear torsional spring and initial pitch deflection [22]. The stability of LCO with respect to control surface freeplay nonlinearity has also been studied both numerically and experimentally [6,25,28]. In the presence of freeplay, the LCO amplitude is strongly dependent on initial conditions. The harmonic balance method, which evaluates steady-state oscillations directly, has been used in studies of LCO to avoid the time-consuming computation of transient behaviour, in particular close to the linear flutter point. A general conclusion is that the accuracy of LCO predictions increases with the number of harmonics retained. Harmonic balance method has been applied for predictions in the context of control surface freeplay [20], where excellent agreement to direct time integration has been demonstrated, and also in the context of LCO arising from both cubic hardening and softening stiffness nonlinearity [12]. In the latter study, when discussing cubic softening stiffness no additional higher-order hardening nonlinearity was considered, and consequently, LCO was not observed at velocities below or above the linear flutter point.
Prediction of nonlinear aeroelastic phenomena such as LCO, which can provide insight into physical mechanisms, requires reliable analysis tools. Time-domain simulation is a direct and general methodology for this task. However, computational cost is then often prohibitive. This is emphasised in, for example, sensitivity studies where large numbers of simulations are required. In an industrial scenario, there will be many parameters under consideration such as multiple payload/mass cases, various gust lengths and shapes, and control laws. A complete evaluation of all parameters can easily result in millions of cases to be simulated. Furthermore, computational cost increases significantly when nonlinear aerodynamics are involved where computational fluid dynamics (CFD) is necessary to capture the relevant flow physics. While being able to provide good predictions of steady-state LCO amplitudes at reduced cost in comparison with timedomain simulation, the harmonic balance method cannot provide information on the relationship between LCO and the initial conditions which give rise to them. This is particularly relevant to subcritical LCO, as stable large amplitude oscillations can occur below the bifurcation point and these are dependent on initial conditions.
Consequently, reduced-order modelling should be exploited. One approach to construct reduced-order models (ROM) is by using data obtained from time-or frequency-domain solutions of the full-order system. Proper orthogonal decomposition [21] uses samples of the full-order dynamical response to form a set of basis functions that can represent the system unknowns. Typically, the snapshots are taken during a dynamic history of interest where critical information relating to the dominant system dynamics can be obtained. The benefit of this method is that the models can be specifically tuned to capture relevant physical behaviour. However, the need to generate samples incurs the computational cost of running full-order simulations.
An alternative method is to construct reduced models based on eigenmode decomposition. The general idea of this method involves using a small number of eigenmodes of the system's Jacobian matrix as a basis to project a Taylor series expansion of the nonlinear full-order residual function. The determination of eigenmodes for aerodynamic models based on linear theory can be achieved using standard eigensolution extraction techniques. For nonlinear aerodynamic models such as the Euler or Navier-Stokes equations on the other hand, the dimension of the system is large, easily exceeding millions of unknowns for realistic problems. Thus, it has been argued that direct extraction of an eigensolution from such system is prohibitive, making proper orthogonal decomposition the preferable method [10]. However, solving the eigenvalue problem for large CFD systems has been demonstrated successfully for test cases ranging from simple pitchplunge aerofoils to three-dimensional industrial problems of realistic aircraft [2,27]. The approach is based on the Schur complement rearrangement of the complete Jacobian matrix of the fluid and structural equations and the use of preconditioned sparse iterative linear solvers [1].
An application of eigenmode decomposition based on the centre manifold theory was investigated to compute LCO response on a rectangular wing in transonic flow [29] and wing rock motion on a delta wing [3]. In this approach, a single critical eigenmode of the Jacobian matrix is used as the basis for the projection. The physical space is thus spanned by the critical eigenmode only, while the influence of the non-critical space is accounted for by the centre manifold theorem. The application of eigenmode decomposition using multiple eigenmodes as basis for the projection of the fullorder system is later discussed [7]. A unified approach for CFD-based flutter and gust response analysis of realistic aircraft configurations is then presented in [27] using this multiple-mode model reduction framework. This demonstrated the applicability of the method when applied to a CFD context. In general, the approach using the centre manifold theory takes a single critical mode focusing on the prediction of steady-state LCO phenomena as these are often dominated by the critical mode. In situations when transient behaviour is of interest, such as transition into LCO due to an initial or gust disturbance, the multiple-mode framework has greater applicability. Here, the non-critical modes become important and are included in the basis.
In the above works [3,7,29], the Taylor expansion of the nonlinear residual includes terms up to third order. Consequently, the resulting nonlinear ROM is capable of capturing up to third-order behaviour of the nonlinear system, while any higher-order dynamics will be lost due to this truncation of the Taylor series. If the nonlinear residual contains higher-order terms with respect to the system states, then higher-order terms in the Taylor expansion are necessary. In other words, a higher-order ROM is required to capture higher-order nonlinear dynamics. This is the key point that must be addressed for a ROM to predict subcritical LCO behaviour.
Subcritical LCO prediction by model reduction remains largely unexplored. Two studies have been identified in the literature. Prediction of a subcritical LCO has been investigated using a reduced-order cyclic method for a pitch-plunge aerofoil in incompressible flow [4,13]. The method combines the techniques of proper orthogonal decomposition and computation of LCO in the modal subspace using an orbital approach and allows for the direct computation of LCO. The aerofoil model here exhibits a subcritical Hopf bifurcation through its nonlinear cubic softening and quintic hardening pitch spring. Prediction of subcritical LCO has also been achieved using proper orthogonal decomposition for the AGARD 445.6 wing [5]. The aerodynamics is modelled using inviscid CFD.
This paper builds on previous work by the authors and presents an assessment of the nonlinear model reduction for the prediction of LCO arising from concentrated structural nonlinearity in the form of polynomial stiffness. First, the existing framework is extended to include up to fifth-order terms in the Taylor expansion allowing, in principle, for the prediction of higherorder nonlinear dynamics. The aim is to determine whether subcritical LCO prediction can be achieved using the resulting ROM. Then, assessments are made comparing the framework based on a single critical mode and multiple modes in their range of applicability and computational cost. A nonlinear two degreeof-freedom pitch-plunge aerofoil structural model coupled with linear aerodynamic theory is formulated here as the test case for the discussion.

Nonlinear full-order model
The fully coupled nonlinear model describing the dynamics of an aeroelastic system can be represented in semi-discrete state-space form. Denote by W the ndimensional state-space vector partitioned into structural states W s and aerodynamic states W f . Written as a set of first-order ordinary differential equations, the state-space equations arė where R is the nonlinear residual vector corresponding to the unknowns W , while Θ is a vector of independent system parameters. The system has a reference equilibrium point W 0 for given constants Θ 0 . At this equilibrium point, the residual vanishes with R (W 0 , Θ 0 ) = 0. Define w = W − W 0 as the increment in the statespace vector with respect to an equilibrium solution. Define also θ = Θ − Θ 0 as the increment in the system's independent parameters with respect to the equilibrium values. The nonlinear residual in Eq. (1) can then be expanded in a multivariate Taylor series about the reference equilibrium point with respect to the system states W and parameters Θ as where A = ∂ R/∂ W is the system Jacobian matrix and F represents all higher-order derivatives in W . Subscript Θ denotes differentiation with respect to it. While only first-order derivatives in Θ are retained here for clarity, higher-order derivatives are included as well in the subsequent discussion. Note that while the residual at the equilibrium is zero, its derivative with respect to the parameters is not in general.
The function F is explicitly written as where B to E are multilinear vector functions of higherorder derivatives. Here functions up to fifth order with respect to the arguments w are retained. More specifically, evaluated about the equilibrium point (indicated by subscript 0), these are Note that B to E are symmetric multilinear functions with respect to their arguments. The multivariate Taylor expansion of Eq. (2) is the starting point for the model reduction formulation described hereafter.

Multiple-mode nonlinear model reduction
In the model reduction approach using multiple eigenmodes, the full-order system is projected onto a small basis of m eigenvectors of the Jacobian matrix A evaluated at the equilibrium point. The eigensolutions of the Jacobian matrix are complex-valued in general. Such eigensolutions exist, for example, as modes of structural vibration. If the aerodynamics is modelled using CFD, complex-valued fluid modes can be found as well. Otherwise, for linear aerodynamics such as discussed in this paper, eigensolutions associated with the fluid unknowns are purely real-valued. Often, a suitable choice is then to retain lower frequency, weakly damped modes which are associated with large amplitudes and hence dominate the system dynamics.
The set of right eigenvectors φ i is obtained by while the corresponding adjoint problem gives the set of left eigenvectors ψ i . The superscript H denotes the conjugate transpose (i.e. Hermitian). The right and left eigenvectors are used to form the corresponding modal matrices, denoted by Φ and Ψ , It is convenient to scale the eigenvectors to satisfy the biorthonormality conditions, where matrices I and O are the identity matrix and a zero matrix, respectively. The biorthonormality conditions also provide the following results where Λ ∈ C m×m is a diagonal matrix containing the eigenvalues. 1 1 Note in the case of a real-valued eigensolution, the biorthonormality conditions can no longer be satisfied as φ =φ. These eigenvectors are then scaled such that ψ H φ = 1 2 giving ψ H Aφ = 1 2 λ which is convenient in order to use consistent notation when dealing both with real-and complex-valued modes.
The vector w is represented by a small set of m eigenvectors using the following coordinate transformation where z ∈ C m is the state-space vector governing the dynamics of the reduced-order nonlinear system. Essentially, w is represented as a linear combination of right eigenvectors with z as the time-dependent amplitude. The nonlinear ROM is then formed by substitution. Premultiplying each term by the Hermitian of the left modal matrix, the nonlinear reduced formulation takes the forṁ Contained in function F, the functions B to E depending on the state variables w have to be transformed as well to be written in terms of z. Observe while the biorthonormality conditions of Eq. (9) hold true for the Jacobian matrix A, they do not apply to A Θ . The evaluation of F Θ simply requires the terms B Θ to E Θ which means differentiating all the functions in Eq. (3) with respect to Θ. Once the equilibrium point and eigensolution are determined, the terms of the reduced formulation involving full-order operations only need to be calculated once and stored initially. Operations involved, when solving Eq. (11) for arbitrary parameter changes, scale with m rather than n. A special case arises when a single critical mode is used as the basis with m = 1. The critical mode corresponds to a Hopf bifurcation where a single pair of complex conjugate eigenvalues of the Jacobian matrix crosses the imaginary axis. Then, Eqs. (5) and (6) reduce to a single eigenvalue problem, and the evaluation of the vector functions as shown in Eqs. (49)-(52) becomes straightforward as all summations vanish. It is expected that steady-state LCO is dominated by the critical mode as the contribution from all other modes are small eventually. Hence, the ROM constructed based on this single mode is able to predict the steady-state LCO amplitude. However, this ROM would fail to predict any transient behaviour where the contributions from non-critical modes have significant influence. Interestingly, for the current test cases with a trivial equilibrium solution, it can be shown that the centre manifold approach [15] becomes equivalent to such critical mode ROM when expanding the Taylor series up to third order only.

Discussion of computational cost
Now the cost of forming and integrating the reduced formulation is discussed in more detail. As explicitly shown in "Appendix 2", each function of increasing order of derivative B to E involves increasing numbers of terms for its evaluation when written as function of z. This relation is presented in Fig. 1. The symmetry properties of the multilinear functions, e.g. B (x, y) = B ( y, x), C (x, y, y) = C ( y, x, y) = C ( y, y, x) and so on, can be exploited to gain two advantages. First, fewer terms are required to construct the ROM initially, and secondly, fewer operations are performed for the Consider the test cases discussed hereafter. All the terms of the ROM are evaluated analytically as a simple pitch-plunge aerofoil structural model with linear aerodynamics is presented. However, if the physical model is more complex, such as a finite element model for the structure and CFD for nonlinear aerodynamics, then analytical evaluation is usually not possible. Thus, a matrix-free approach using finite differences is necessary to construct the functions of higher-order derivatives requiring only the evaluation of the residual function. The number of residual evaluations per finite difference approximation depends on the order of derivative to be approximated and the order of the finite difference. As an example in "Appendix 2", the secondorder central difference schemes in Eq. (53) involve between two to six residual evaluations for the functions B to E with a single argument.
In addition, the finite difference approximations of the functions B to E are defined for a single argument. For a function with mixed arguments, i.e. different eigenvectors, a set of identities must be used to place the evaluation in a suitable form with a single argument, e.g.
and so on. Thus, several finite differences are required for a term with mixed arguments. Additionally, the arguments of the terms are, in general, complex-valued. In situations where complex arithmetic is not possible, such as industrial finite element or CFD codes, these identities must be used again to isolate the real and imaginary parts of the arguments which results in additional finite difference operations [15].
It can be concluded that the total number of residual evaluations required to approximate a single function of a higher-order derivative is significantly higher than the number of terms shown in Fig. 1. Using many modes in the ROM basis as well as higher-order Taylor expansion, millions of full-order residual evaluations could be required making the construction of the ROM prohibitive. Few, well-selected, dominant modes must therefore be chosen.
Time integrating the ROM demands the evaluation of the dynamics of the reduced state z in Eq. (11). The computational cost associated with running the ROM is directly related to the results presented in Fig. 1. Once the ROM terms are evaluated either analytically or by finite differences, projection with Ψ H gives complexvalued terms, the number of which scales with m. Integration methods for ordinary differential equations, such as schemes of the Runge-Kutta family, are readily available working with complex arithmetic as well. Using industrial finite element or CFD codes, the cost of integrating the governing equations of the reducedorder model will remain small in comparison. It is important to note that the construction of the ROM from the full-order system is as expensive as running the full-order simulation. This is the case with most ROM formulations of any flavour. However, the current formulation allows the construction of the ROM once and for all, while subsequent parameter changes require minimal additional cost.

Two degree-of-freedom aerofoil model
As shown in Fig. 2, a standard two degree-of-freedom aerofoil model elastically supported in plunge and pitch is considered. The plunge deformation is indicated by h positive downward, and the pitch angle α is positive nose-up. The terms m, b and a h are the aerofoil mass, the semi-chord length and the nondimensional distance from the mid-chord to the elastic axis, respectively. The term I α represents the moment of inertia about the elastic axis given in dimensionless form as radius of gyration r α = I α /mb 2 . The static moment about the elastic axis is S α = mx α b where x α b denotes the offset of the centre of gravity from the elastic axis. The terms K h and K α denote the spring stiffness coefficients in bending and torsion, respectively. While the Fig. 2 Two degree-of-freedom pitch-plunge aerofoil bending stiffness is linear, the torsional spring is modelled by a fifth-order polynomial. The natural uncoupled frequency in plunge is denoted ω h = √ K h /m and similarly for pitch ω α = √ K α /I α . The frequency ratio is given byω = ω h /ω α . The structural damping coefficients (not explicitly given in the figure) are denoted by ζ ξ and ζ α , respectively.
Using the nondimensionalisation whereū, as a nondimensional representation of the freestream velocity U , is referred to as the reduced velocity, the equations of motion are derived as with( ) indicating derivatives with respect to nondimensional time τ and μ = m/ πρb 2 describing the aerofoil-to-fluid mass ratio (where ρ is the freestream density). The nonlinear torsional spring polynomial f α is given by where β α3 and β α5 are the specific stiffness coefficients corresponding to the cubic and quintic nonlinearity. Figure 3 presents the torsional spring stiffness f α for the linear case as well as the two baseline configurations discussed below. In this study, the aerodynamic coefficients of lift C L and moment C M about the elastic axis originate from  Fig. 3 Torsional spring polynomial f α two physical sources, i.e. contributions due to the wing motion and due to gust disturbance. The coefficients can then be written as where the superscripts m and g refer to contributions from aerofoil motion and gust, respectively.
To proceed with the aerodynamic formulation, the nondimensional downwash at the three-quarter chord point due to the structural motion in pitch and plunge is introduced as Then the lift is obtained following [11] as where the Wagner function Φ w (τ ) is approximated by with constants Ψ 1 = 0.165, Ψ 2 = 0.335, ε 1 = 0.0455 and ε 2 = 0.3 following [14]. The Wagner function describes the ratio of transient to steady-state lift (with circulatory origin) for a general aerofoil motion. Note that Φ w (0) = 0.5. Thus, while terms involving the Wagner function denote circulatory contribution, the remaining terms are of non-circulatory origin. Accordingly, the pitching moment coefficient is defined as The lift and moment coefficients are modified here by additional terms to account for an arbitrary gust excitation W g (τ ) as discussed in [8]. These terms are explicitly given by where the Küssner function Ψ k (τ ) is approximated in the form with constants Ψ 3 = 0.5792, Ψ 4 = 0.4208, ε 3 = 0.1393 and ε 4 = 1.802 as given in [19]. The Küssner function describes the ratio of transient to steady-state lift for an aerofoil penetrating a sharp-edged gust. Note that Ψ k (0) = 0, indicating that the initial lift build up due to the aerofoil entering a sharp-edged gust is zero. The above expressions for aerodynamic lift and moment coefficients are integro-differential equations. The terms involving the convolution integral can be eliminated by introducing auxiliary variables [16] the dynamics of which are evaluated aṡ using the Leibniz integral rule. Both the structural equations and the aerodynamic force and moment expressions depend on the same shared system states, which are X s for the structural degrees of freedom and W f for the augmented aerodynamic states, Combining the equations and collecting the coefficients of common terms, a set of governing ordinary differential equations describing the dynamics of the structural system are obtained. This is expressed in matrix-vector form as (28) where the terms M, C and K are the effective mass, damping and stiffness matrices containing structural and aerodynamic contributions, k N is a nonlinear vector arising from the polynomial stiffness, and f a arises from the influence of initial conditions upon the unsteady aerodynamic forces. The matrix D f relates the structural equations to the augmented aerodynamic states. Similarly, the aerodynamics system in Eq. (25) can be formulated aṡ with matrix A f f relating the fluid unknowns to their first temporal derivatives and matrix A f x coupling the fluid equations to the structural degrees of freedom. The vector f g describes the influence of the external gust disturbance on the aerodynamic equations. The explicit form of these matrices and vectors is given in "Appendix 1". In the final step, Eqs. (28) and (29) are recast in a coupled first-order ordinary differential equation of the general form as given in Eq. (1) where the unknowns are partitioned into structural and aerodynamic contributions and the system parameters Θ consist of the reduced velocityū and the gust velocity W g . The corresponding residual vector R is given by where the matrix A L contains the linear contributions only and is defined as The matrix blocks Aẋ f = −M −1 D f and A f x couple the structural and aerodynamic equations. The vec-tors b N , b a and b e denote the nonlinear terms, aerodynamic contribution due to initial conditions and external inputs, respectively, Note that the equilibrium point for the aerofoil model presented herein is the trivial solution.
The set of first-order ordinary differential equations in time is the starting point for the application of the model reduction presented in Sect. 2.2. As the aerodynamics are linear following the theories of Wagner and Küssner, the higher-order terms in the model reduction arise only from the nonlinear terms contained in vector b N . The functions B to E are evaluated analytically herein, which is straightforward for the pitchplunge aerofoil structural model with linear aerodynamics, and have been tested by comparison with finite difference evaluations. Since the nonlinear residual, the Jacobian matrix and the higher-order derivatives are all known analytically, their derivatives with respect to the reduced velocityū follow immediately as well.
To account for the gust input in the reduced model, the derivative of the external input vector b e with respect to the gust velocity W g is required adding a contribution to R Θ in Eq. (11). As is explicitly shown in "Appendix 1", the only nonzero entries in vector b e take the values of the gust velocity and the derivative with respect to gust velocity becomes trivial.
Operations with higher-order derivatives, such as vector multiplication, become prohibitive using dense array storage schemes. Since these higher-dimensional arrays are very sparse due to the discrete polynomial form of the structural nonlinearity, a higher-order extension of the well-known compressed sparse row format [23] has been devised for the work presented herein. All sparse array operations have been compared to equivalent dense operations for validation, and very significant savings in computing time are achieved.

Results
The structural model used for the analysis is described by the parametersω = 0.2, μ = 100, a h = −0.5, x α = 0.25 and r α = 0.5. While the plunge spring is linear, two different nonlinear springs in the pitch degree of freedom are considered. The first configuration, following [17], is a cubic hardening spring with β α3 = 3, while the second configuration, following [22], is cubic softening and quintic hardening with β α3 = −3 and β α5 = 20. No structural damping is assumed in the analysis herein. From here on, direct integration of Eq. (1) with the residual defined in Eq. (31) is referred to as reference solution.

Linear stability analysis
While the model reduction approach based on multiple modes can be constructed at any value of the bifurcation parameter (i.e. the reduced velocityū), the reduction approach based on a single mode requires the critical eigensolution. Locating the system's bifurcation point is achieved by varying the reduced velocityū and solving for the eigenvalues of the Jacobian matrix A. The system has two complex conjugate pairs of eigenvalues corresponding to the two structural degrees of freedom and four purely real-valued eigenvalues corresponding to the aerodynamic states. There are additionally two purely real-valued eigenvalues corresponding to the gust states w 5 and w 6 in Eq. (27). The point where the real part of a pair of complex conjugate eigenvalues vanishes indicates the Hopf bifurcation.
The tracing of the eigenvalues, originating in the structural part of the coupled system, with respect to the reduced velocity is shown in Fig. 4. The results are confirmed with the results obtained in [17] with a bifurcation point found at aboutū L = 6.285. This is the point where the reduced-order models are constructed in the subsequent discussion. The construction is done once and for all, while arbitrary parameter changes (i.e. reduced velocity and gust disturbances) can then be included at negligible extra cost.

Supercritical bifurcation
The first scenario considered here is LCO arising from a supercritical bifurcation corresponding to the case of a cubic hardening spring with β α3 = 3. The nonlinear reduced-order models based both on multiple modes and the critical mode are constructed at the bifurcation point. Since the nonlinearity is cubic only, the terms up to C must be evaluated in the ROM construction, Limit-cycle oscillation induced by an initial disturbance is studied first. The reduced velocity is set to five per cent above the bifurcation point, and the system is excited by an initial disturbance in pitch of 5 • . Due to the nonlinear cubic hardening restoring moment in pitch, the system exhibits a stable LCO for reduced velocities above the flutter point. Figure 5 shows the time history of the LCO as predicted by both critical and multiple-mode ROM compared with the full-order reference solution. The two modes used in the construction of the multiple-mode ROM originate from the structural vibration problem. Excellent agreement with the full-order solution is observed throughout. The crit- ical mode ROM predicts the amplitude of the LCO well, although discrepancies are found during the transient response up to approximately 100 time units. Since the critical mode ROM formulation relies on information exclusively from the single critical mode, it is clear why this type of behaviour is observed. In the transition period to a steady-state limit-cycle response, the dynamics of the system are influenced by contributions from the non-critical structural mode as well until it is damped out in time. Once the steady-state limit-cycle is reached, the dynamics are dominated by the single critical mode. A point should also be noted about the discrepancy in the initial conditions between the solutions in Fig. 5. Specifically, in the full-order reference solution the initial pitch is set to 5 • , while after transformation to and from the reduced space the initial pitch obtained is about 3 • for the ROM solution. Similar behaviour is found in the plunge degree of freedom, which is not reproduced as zero. The critical mode ROM (m = 1) depends exclusively on the critical mode, and the reduced space spanned by this single mode cannot fully represent the corresponding initial condition in physical space. Making the transformation to the reduced space with z (0) = Ψ H w (0) to obtain the initial condition of the reduced state variable z, and back to physical space with w (0) = Φ z (0) +Φz (0), information is lost and the difference in the initial condition arises. Indeed, retaining all modes of a physical system, the expression ΦΨ H +ΦΨ H = I is satisfied otherwise. The multiple-mode ROM (m = 2) solution is already significantly improving the results and thus confirming the choice of keeping two modes.
Limit-cycle oscillation can also be induced by gust disturbance. Figure 6 shows the same system disturbed by a 1-Cosine gust. The gust profile is defined as with the dimensionless parameters being gust intensity W 0 = 0.1 (i.e. 10 % of freestream speed), gust wavelength L g = 20 semi-chord lengths and gust initial time τ 0 = 10. For the construction of the multiple-mode ROM, the effect of including two additional modes has been assessed as only including the structural vibration modes is not sufficient for accuracy during the gust excitation. The first additional mode corresponds to the real-valued eigenvalue of λ = −0.1393 which is the lower time constant used in the approximation of the Küssner function. This mode is demonstrated to be dominant in coupling the structural response to the gust input as discussed in [7]. The multiple modes ROM constructed using these three modes shows good agreement with the full-order solution. However, differences in the transient response up to 100 time units are still observed, particularly in the plunge response. This prediction is significantly improved by including an additional fourth mode of a real-valued eigenvalue λ = −0.03178 corresponding to the Wagner aerodynamics. Its value changes with respect to reduced velocity and can always be identified as it is close to the value of the lower time constant used in the approximation of the Wagner function, while the remaining eigenvalues corresponding to the Wagner aerodynamic states appear as repeated exact values of the time constants. The critical mode ROM predicts the steady-state LCO well but fails during the gust disturbance. The consideration here is that the critical mode ROM lacks the gust mode and is thus unable to correctly capture the transition into the LCO. By producing a series of time histories, it is possible to trace the plots in Figs. 7 and 8 showing the evolution of the stable limit-cycle amplitude with respect to the reduced velocity. Results are compared with predictions presented in [17] which were obtained via full-order time-domain simulations. In Fig. 7, the ROM construction retains only first-order derivatives in the reduced velocity ratio u/u L non-dimensional plunge amplitude 1 1.05 1. Taylor expansion with respect to the reduced velocity, while in Fig. 8 every term up to third order is included. Good agreement is observed between the different predictions of reduced models and reference solution. It should be noted here that the current reference solution reproduces the results from [17] accurately. In Fig. 7, the plunge and pitch amplitude predictions by the reduced models show generally good agreement with each other, while larger amplitude discrepancies to the reference solution are observed at higher flow velocities. In contrast, the amplitude predictions at higher flow velocities as shown in Fig. 8 are in much closer agreement with the reference solution. This demon-reduced velocity ratio u/u L non-dimensional plunge amplitude 1 1.05 1. strates the improvement when including higher-order derivatives with respect to the reduced velocity. As is shown by the results presented, a source of error in the reduced-order formulation is the truncation of the Taylor series expansion, while the overall quality of the prediction is improved by including higher-order derivatives with respect to the system parameter. The reduced model is also limited by the set of basis vectors used in its construction, and any information lost with the excluded vectors is an inherent and necessary concession. Thus, identifying the dominant modes is critical for constructing an accurate and representative reduced-order model. It should also be noted that in the high velocity region the amplitudes in pitch and plunge are large, exceeding for instance 20 • in pitch angle. This is unrealistic in the context of the linear aerodynamic theory applied here due to the appearance of nonlinear aerodynamic phenomena such as massive boundary layer separation. However, the results presented herein are nevertheless solutions of the models used to approximate real physics.
Finally, to demonstrate the robustness of the modelling approach, we present results for various values of the nonlinear term β α3 including subcritical flutter instabilities, rather than limit-cycle responses, for negative stiffness constants. Figure 9 compares the fullorder reference solutions with results from the reduced formulation based on the structural modes only and retaining derivatives up to third order with respect to the reduced velocity. Excellent agreement is found throughout. Note that the subcritical flutter boundary (dashed line) is indicated by the critical initial pitching angle α(0) to cause the instability. All remaining initial values are set to zero.

Subcritical bifurcation
Subcritical limit-cycle responses corresponding to the case of a cubic softening combined with quintic hardening spring with β α3 = −3 and β α5 = 20 are considered next. The general mechanism is that the cubic softening nonlinearity contributes to a destabilising effect, while at large amplitudes the quintic hardening nonlinearity acts to constrain the system response from diverging. From Fig. 3, a stable LCO amplitude in pitch deflection of about 20 • -25 • is expected since here the hardening effect exceeds the softening term. In the presence of the softening nonlinearity, limit-cycle responses can occur at velocities below the linear flutter point. As above, the reduced models are constructed once and for all at the bifurcation point. Since the discrete nonlinearity now extends to quintic order, terms up to E must be evaluated. Note that functions B and D are zero as the equilibrium solution is trivial.
Most importantly, the reduced-order formulation is capable of predicting subcritical limit-cycle response resulting from structural softening nonlinearity as a consequence of extending the formulation to include up to fifth-order derivatives in the state variables. A lowerorder formulation, such as of cubic order, used previously cannot predict this response. Figure 10 shows the transition to LCO at a reduced velocity three per cent below the linear flutter point. The initial condition here is a disturbance in pitch of 13 • , the reason of which will be explained below. Similarly to the case of supercritical bifurcation, the modes used for the construction of the multiple-mode ROM are the two complex-valued eigenmodes corresponding to the structural degrees of freedom. Excellent agreement with the full-order solution is found. The critical mode ROM predicts the steady-state limit-cycle response well but shows significant amplitude differences in the transition period. This is in line with the discussion for the case of supercritical bifurcation in Fig. 5 and due to the limited information contained in the critical mode.
As with the supercritical case, subcritical LCO can also be induced by a gust disturbance, now even below the flutter point. Figure 11 shows the system subjected to the 1-Cosine gust profile with gust intensity W 0 = 0.1, wavelength L g = 20 and initial time τ 0 = 10. Corresponding to the supercritical case, the same two additional modes are assessed in the prediction of the transient behaviour. The conclusion remains unchanged. The best prediction is obtained by including the fourth mode corresponding to the purely real-valued eigenvalue of λ = −0.03178. As expected, the critical mode ROM predicts the steadystate limit-cycle amplitude reasonably well but fails to accurately capture the response during the gust disturbance phase. In the case of subcritical bifurcation, the presence of the quintic hardening stiffness constrains the system at high amplitudes, and consequently, stable large amplitude LCO is observed in Figs. 12 and 13 as indicated by the characteristic stable branch. The stable limit-cycle pitch amplitude is confirmed with results presented in [22]. The stable branch extends to reduced velocities below the bifurcation point. In this region, if the system is subjected to large enough disturbances the response will jump into the stable branch. Accordingly, if the disturbances are not large enough to incite the jump, the response will decay to zero. The critical initial pitching angle and plunge deflection causing Transient leading to subcritical LCO due to gust disturbance,ū = 6.097. a Plunge response, b pitch response the transition to stable large amplitude LCO are highlighted by the respective, as referred to herein, unstable branch. Consequently, to obtain stable LCO response below the flutter point, the initial condition required is dictated by the unstable branch. This justifies the choice of the initial pitch disturbance of 13 • stated above for the response obtained in Fig. 10 as it is well above the unstable branch. For completeness, several other values above the unstable branch are tested and result in the same steady-state LCO amplitude. It is important to note that the stable branch is independent of the initial conditions applied.
In the following, the evolution of the stable limitcycle amplitude with respect to the reduced velocity is reduced velocity ratio u/u L non-dimensional plunge amplitude  First, the critical mode ROM does predict the subcritical limit-cycle response even though discrepancies are observed. Also, the stable LCO amplitude as predicted by the multiple-mode ROM, based on the two structural modes only, shows good agreement with the full-order solution except for higher values of reduced velocity. Here the multiple-mode ROM based on two structural modes fails. This behaviour is not improved when including higher-order derivatives with respect to the reduced velocity and is thus not shown in Fig. 13.
Since it is expected that the reduced model forms a better representation of the full-order system with increasing number of modes, excellent agreement with the reference solution is found throughout when including the aerodynamic mode. This emphasises again the importance of identifying the dominant modes.
Secondly, as is found for the supercritical case, including higher-order derivatives in the reduced velocity improves the ROM predictions. This can be seen clearly from the turning point where the unstable branch transitions into the stable branch. In Fig. 12, the turning point as predicted by the reduced models are clearly offset from the reference solution, while in Fig. 13 excellent agreement between the solutions is observed. Note that for Fig. 13 the unstable branch in pitch as predicted by the critical mode ROM better agrees with the reference solution compared to the unstable branch in plunge. As the critical mode originates in the pitch degree of freedom, it lacks information to fully represent the plunge response for the critical mode formulation. The multiple-mode ROM on the other hand includes both the pitch and plunge degrees of freedom giving an overall excellent agreement to the reference solution.
As with the supercritical configuration in Sect. 4.2, we include a parameter study to demonstrate the robustness of the modelling approach. Figure 14 presents two additional scenarios of a fifth-order nonlinearity, besides the baseline configuration. The reduced-order formulation includes the three modes discussed above and expands the residual with respect to the reduced velocity up to third order. Overall good agreement with the reference solution is found for all cases.

Conclusions
In this work, a nonlinear model reduction approach based on eigenmode decomposition and projection is presented focusing on the prediction of sub-and supercritical limit-cycle behaviour. The model reduction is formulated to include up to fifth-order derivatives in the Taylor expansion of the nonlinear full-order residual function. The underlying motivation is such that higherorder nonlinear behaviour of the full-order system can therefore be captured with the extended reduced model formulation. The methods are applied to a two degreesof-freedom pitch-plunge aerofoil structural model coupled with linear aerodynamics. Stiffness nonlinearity is introduced into the pitch degree of freedom in polynomial form up to fifth order such that the system exhibits the desired limit-cycle behaviour. Including multiple modes in the basis, used for projection of the full-order system, is unnecessary in the case of a supercritical limit-cycle oscillation if interest lies in the amplitude prediction. Results using only the critical mode for model reduction are in excellent agreement with the reference solution. If transient behaviour is important however, as is the case when gust disturbance is discussed, multiple modes are required. The situation changes for subcritical limitcycle behaviour. Even if the interest is in the amplitude only, the critical mode alone is not sufficient. While the reduced model based on the critical mode does predict the subcritical limit-cycle response, multiple modes are mandatory to resolve the discrepancy between the reduced model prediction and the reference solution. If the gust-induced transient is important, the same conclusion as for the supercritical case can be reached. Finally, the order of expansion with respect to the parameters is important as well with higher-order formulations giving superior results.
Analysis of computational cost shows that if the number of modes used in the reduced model construction is large, then forming the reduced model becomes increasingly expensive. This is particularly true when the physical problem is complex and analytical expressions of the Jacobian matrix and higher derivatives are not possible such that finite difference approximations must be used. The number of residual evaluations required for the finite differences can easily grow to the order of millions for the reduced model formulation including fifth-order terms. In the case where only the single critical mode is used, the associated com-putational cost is always feasible. The identification of few dominant modes is thus an important topic which needs to be addressed further. These symmetry properties can be exploited to reduce the total number of individual terms required. For convenience, define the following auxiliary terms for the second-and third-order functions b rr = B φ r , φ r z r z r + 2B φ r ,φ r z rzr + B φ r ,φ r z rzr b rs = B φ r , φ s z r z s + B φ r ,φ s z rzs + B φ r , φ s z r z s + B φ r ,φ s z rzs (46) c rrr = C φ r , φ r , φ r z r z r z r + 3C φ r , φ r ,φ r z r z rzr + 3C φ r ,φ r ,φ r z rzrzr + C φ r ,φ r ,φ r z rzrzr c rrs = C φ r , φ r , φ s z r z r z s + C φ r , φ r ,φ s z r z rzs + 2C φ r ,φ r , φ s z rzr z s + C φ r ,φ r , φ s z rzr z s + 2C φ r ,φ r ,φ s z rzrzs + C φ r ,φ r ,φ s z rzrzs c rss = C φ r , φ s , φ s z r z s z s + C φ r , φ s , φ s z r z s z s + 2C φ r , φ s ,φ s z r z szs + C φ r ,φ s ,φ s z rzszs + 2C φ r , φ s ,φ s z r z szs + C φ r ,φ s ,φ s z rzszs c rst = C φ r , φ s , φ t z r z s z t + C φ r , φ s ,φ t z r z szt + C φ r ,φ s , φ t z rzs z t + C φ r ,φ s ,φ t z rzszt + C φ r , φ s , φ t z r z s z t + C φ r , φ s ,φ t z r z szt + C φ r ,φ s , φ t z rzs z t + C φ r ,φ s ,φ t z rzszt The total number of terms based on these new formulations, as given in Table 1, as well as the direct approach is illustrated in Fig. 1.
As discussed in Sect. 2.3, finite difference evaluations are often required to form the terms when analytical expressions of the derivatives are not available. For example, the following second-order central difference schemes in one argument vector x can be used where R l = R (W 0 + lεx).