The construction of non-linear normal modes for systems with internal resonance

Abstract A numerical method, based on the invariant manifold approach, is presented for constructing non-linear normal modes for systems with internal resonances. In order to parameterize the non-linear normal modes of interest, multiple pairs of system state variables involved in the internal resonance are kept as ‘seeds’ for the construction of the multi-mode invariant manifold. All the remaining degrees of freedom are then constrained to these ‘seed’, or master, variables, resulting in a system of non-linear partial differential equations that govern the constraint relationships, and these are solved numerically. The computationally-intensive solution procedure uses a combination of finite difference schemes and Galerkin-based expansion approaches. It is illustrated using two examples, both of which focus on the construction of two-mode models. The first example is based on the analysis of a simple three-degree-of-freedom example system, and is used to demonstrate the approach. An invariant manifold that captures two non-linear normal modes is constructed, resulting in a reduced order model that accurately captures the system dynamics. The methodology is then applied to a larger order system, specifically, an 18-degree-of-freedom rotating beam model that features a three-to-one internal resonance between the first two flapping modes. The accuracy of the non-linear two-mode reduced order model is verified by comparing time-domain simulations of the two DOF model and the full system equations of motion.


Introduction
In order to obtain accurate reduced order models for non-linear systems, "non-linear modal analysis" has been proposed as an analogy to its linear counterpart. The concept of non-linear normal modes (NNM) was initiated by Rosenberg [1] and subsequently considered by a number of other investigators and their co-workers, including Rand [2,3], Nayfeh [4][5][6], and Vakakais [7][8][9], who have studied the existence, construction, stability, and bifurcations of NNMs. The definition of NNMs was generalized by Shaw and Pierre [10,11] through the introduction of invariant 1 manifolds. The works by Nayfeh [12], Vakakis et al. [13], and Camillacci [14] offer overviews of the various methods of constructing NNMs.
A NNM invariant manifold is a two-dimensional surface in the system phase space that is tangent to the corresponding linear modal eigenspace at the equilibrium point [10]. In order to parameterize these manifolds for vibratory systems, a single pair of state variables in linear modal coordinates (typically a modal displacement-velocity pair, or a modal amplitude and phase) are chosen as master coordinates for an individual NNM. Then, all the remaining degrees of freedom (DOF), the so-called slave coordinates, are constrained to these master coordinates in a particular manner, dictated by the equations of motion. The nonlinear partial differential equations (PDEs) describing the geometry of the manifold are produced using an approach that follows center manifold construction. Based on this methodology, a numerical framework for constructing NNMs, namely a Galerkin projection method [15], has been proposed and effectively applied to a variety of non-linear systems, including systems with non-smooth restoring forces [16], as well as systems with non-proportional damping forces, non-symmetric non-linearities, and gyroscopic terms [17]. Once the NNM invariant manifold is obtained, motions on it are governed by the dynamics of the corresponding master coordinates, which are described by two first-order ordinary differential equations. Note that these equations of motion are valid only for initial conditions on the invariant manifold.
Unlike linear modes, NNMs will interact during a general motion that is initiated by general initial conditions. Moreover, an invariant manifold approach that is based on a single mode reduction will break down in the presence of internal resonances between the master and any slave coordinates [18]. Hence, a non-linear normal multi-mode methodology is required if one is interested in the multi-mode responses of non-linear systems. Previous studies of NNMs with internal resonances were primarily based on perturbation methods [6,8] or polynomial series expansions with the invariant manifold approach [18,19]. These are applicable only in the weakly non-linear regime. In order to obtain accurate reduced-order models for non-linear systems with internal resonances in strongly non-linear amplitude regimes, a new method for constructing invariant manifolds is proposed in this paper, as follows.
For an n-DOF non-linear system with M modes involved in an internal resonance (or, more generally, with M modes to be retained for any reason), the multi-mode invariant manifold can be defined and obtained numerically in terms of 2M master coordinates (M displacement-velocity pairs or M amplitudes and phases). The procedure is outlined as follows. First, a transformation to polar coordinates is applied, as in the single-mode expansion case proposed by Pesheck et al. [15], to each pair of master coordinates. Then, the constraint functions for the slave coordinates are expressed in terms of these master coordinates, which are M amplitudesā and M phases. Since the master coordinates are represented by more than one amplitudephase pair, the computation of the manifolds is much more complicated than in the single-mode case. To address this problem, a computational approach is proposed, which combines Galerkin projections in the phase coordinates and finite difference discretizations in the amplitude coordinates. Using this methodology, a 2M-dimensional invariant manifold can, in principle, be constructed for the system, and motions on this manifold are governed by a set of 2M first-order differential equations in the master coordinates.
It shall be noted that the procedure presented in this paper is utilized to examine the non-linear modal dynamics for the M modes involved in the internal resonance. In order to obtain a complete description of the system's non-linear modal structure, the remaining n − M NNMs will have to be constructed.
The multi-NNM approach is applied herein to a simple three-DOF system, as well as to a rotating blade model in which transverse motions are non-linearly coupled with axial extensions of the blade. For the latter system, an 18-DOF discretized model derived from linear modal analysis is examined, which features an internal resonance between the first and second flapping modes. Using the multi-NNM procedure, the internally resonant four-dimensional invariant manifold is constructed and the resonant dynamics are shown to be accurately captured by the two-DOF (four state) reduced-order model.
The paper is organized as follows. The class of nonlinear systems under consideration and the formulation of the multi-mode invariant manifold equations are described in Section 2. In Section 3, the solution procedure for the invariant manifold is demonstrated on a 3-DOF example system. In Section 4, the method-ology is applied to an 18-DOF rotating beam system featuring an internal resonance between the first two flapping modes. Finally, some conclusions are drawn in Section 5.

Multi-mode invariant manifolds
The vibratory system considered is an n-DOF autonomous non-linear system, which can be obtained directly from Newton's laws, Lagrange's equations, finite element methods (FEM), or through any assumedmode method. In order to simplify the construction procedure for the invariant manifolds, systems without damping or gyroscopic terms are considered (the method can be generalized to include these effects, but it is much more computationally intensive). For these systems, one can transform to linear modal coordinates, and the original non-linear system can be expressed in the following standard form: where the i s are the linear modal coordinates, the i s are the corresponding linear modal frequencies, and f i ( j ) is the ith non-linear force, which generally depends on all the linear modal coordinates, expressed in terms of the linear modal coordinates. It should be noted that the non-linear forces, the f i ( j )s, are assumed to be independent of the linear modal velocities,˙ j . With this assumption, the computational cost for the construction of the invariant manifolds can be significantly reduced. In order to obtain accurate reduced-order models for non-linear systems with internal resonances, multimode invariant manifolds must be constructed [18]. Here, an invariant manifold is defined as a multidimensional surface spanned by all the linear modal displacements and velocities involved in the internal resonance, such that any motion initiated on the manifold will remain on it for all times. Following this definition, the invariant manifold can be constructed in the following manner.
Let us assume that there are a total of M modes involved in an internal resonance for the non-linear system defined in Eq. (1). These modes are described by a set of indices, denoted as S M . According to the definition of invariant manifolds, each pair of state variables involved in the internal resonance, ( k ,˙ k ) k ∈ S M , are chosen as master coordinates. Then all the remaining DOFs, namely the slave coordinates, are constrained such that they are dependent on the 2M master coordinates. A straightforward expression for the slave constraint functions is This form has been utilized by Pesheck et al. [18] to carry out an asymptotic procedure for the construction of the invariant manifolds. In this approach, the constraint relationships, Eq. (2), are approximated by polynomial expansions in the master coordinates. Using this method, the domain of validity of the resulting solution is limited to some neighborhood of the system's equilibrium position. In addition, this approach does not allow one to systematically control the accuracy of the approximate solution. Hence, a new form for the master coordinates and the slave constraints is used here in order to overcome some of the inherent deficiencies of the polynomial expansion method. For each pair of master coordinates, ( k ,˙ k ) k ∈ S M , a polar coordinate transformation is applied, where k is the kth linear modal frequency, and (a k , k ) are the master coordinates in amplitudephase form. Then, all the slave coordinates are expressed as functions of these amplitude-phase master coordinates, where the slave coordinates are restricted to the domain defined by the M pairs of amplitude-phase variables. This domain is easily bounded, since the phase coordinate is periodic and the positive amplitude region can be dictated to a range of interest during the construction of the invariant manifold, as described below.
With the above polar form for the master coordinates and the slave constraint relationships, the partial differential equations governing the invariant manifold can be obtained as follows. For each pair of slave coordinates, ( i ,˙ i ), we have where i and f i are the ith linear modal frequency and the non-linear force defined in Eq. (1), respectively. In order to eliminate the explicit time dependence,ȧ k and k in Eq. (5) are replaced using the governing equation of motion for each pair of master coordinates, which arė Substituting Eq. (6) into Eq. (5), the governing PDEs for the invariant manifold are found to be In Eq. (7), there are a total (n − M) pairs of equations governing the invariant manifold, in terms of the constraint equations. These equations are non-linear and have to be solved in some approximate manner; here this is done numerically. Once they are solved, the results for all of the slave constraints, P i and Q i for i / ∈ S M , can be substituted into the M pairs of ordinary differential equations governing the dynamics of the master coordinates, Eq. (6). As a result, the response of the original system restricted to the invariant manifold is captured by this 2M-DOF reduced-order model. The procedure for obtaining the numerical solution for the invariant manifolds is described in the context of a simple example, and then applied to a more substantial problem. Fig. 1 depicts a three-DOF mass-spring system with two cubic non-linear springs of coefficients 1 and 2 , attached to masses M 1 and M 3 , respectively. The system parameters are tuned so that the second linear modal frequency is approximately three times the first one, i.e. 2 ≈ 3 1 . Consequently, a three-to-one internal resonance occurs between the first and second linear modes.

A three-DOF example system
Using linear modal coordinates, the system can be transformed to the standard form shown in Eq. (1) where the modal coordinates i are defined by the linear modal transformation: In the presence of an internal resonance between the first two modes, the master coordinates are chosen as the state variable pairs ( 1 ,˙ 1 ) and ( 2 ,˙ 2 ). The corresponding master coordinate index set is The polar coordinate transformation is applied to these state variable pairs according to the definition in Eq.
(3). As a result, the transformed master coordinates are expressed in terms of the amplitude-phase pairs, (a 1 , 1 ) and (a 2 , 2 ). As defined in Eq. (4), the constraint relationships for the slave coordinates are 3 N. X 1 , X 2 , and X 3 denote the displacements of masses M 1 , M 2 , and M 3 , respectively.
The governing PDEs for the invariant manifold are given in Eq. (7) and are rewritten here as A numerical solution scheme is given here to approximate the unknown constraint relationships, P 3 (a 1 , a 2 , 1 , 2 ) and Q 3 (a 1 , a 2 , 1 , 2 ), which define the geometry of the four-dimensional invariant manifold. In the two-dimensional phase region, defined by the constraint relationships P 3 and Q 3 , are periodic in both 1 and 2 . Hence, they can be efficiently approximated by two-dimensional Fourier series. In the two-dimensional amplitude domain, defined by where the upper limits a 1 max and a 1 max are set during the numerical construction procedure, finite difference discretization methods can be used to approximate the unknown constraint equations by a sequence of overlapping polynomials that interpolate P 3 and Q 3 at a set of grid points. It should be noted that the region of the two lines, {(a 1 , a 2 ) | a 1 = 0 or a 2 = 0}, is excluded from the two-dimensional amplitude domain in the finite difference scheme, since this region is not defined in the governing partial differential equations for the invariant manifold. With the combination of finite difference methods and two-dimensional Fourier series expansions, the unknown constraint equations P 3 and Q 3 , can be approximated at each grid point as where (a i 1 , a j 2 ) is the grid point determined by the finite difference scheme in the amplitude region; indices i and j denote the location of the grid point along the a 1 and a 2 directions, respectively; the Fourier terms, F l ( 1 ) and F m ( 2 ), are defined as F l () = cos l−1 2 , l is odd, sin l 2 , l is even (14) 5 and N 1 and N 2 are the number of terms of the Fourier expansions in 1 and 2 , respectively. As can be seen in Eqs. (12)- (14), the total number of the unknown quantities are determined by N 1 and N 2 , as well as the number of grid points. Once the unknown coefficients, the Cs and Ds in expressions (12) and (13), have been obtained at all grid points, the invariant manifold is completely determined in this approximate manner. Given the expression of P 3 and Q 3 at each grid point, the derivatives of the local interpolant are used to approximate the derivatives of P 3 and Q 3 with respect to a 1 or a 2 . Simple two-point backward interpolation gives where are the distances between adjacent grid points along the a 1 and a 2 directions, respectively. The function O(·) denotes that the errors in these approximations are ordersof-magnitude h 1 and h 2 , respectively. If more grid points are used in the approximation for a given amplitude range, higher accuracy is obtained. The approximation of the derivatives, jQ 3 /ja 1 and jQ 3 /ja 2 , is determined using a similar scheme. Along the 1 and 2 directions, the derivatives of the unknown functions P 3 and Q 3 can be easily obtained using the two-dimensional Fourier series expansions given in expressions (12) and (13).
The computational time associated with the construction of the invariant manifold depends on the number of unknown coefficients, which depends on the number of grid points selected for the amplitude variables and the number of harmonics employed in the Fourier series. For this example system, the total number of unknown coefficients, Cs and Ds, can be re-duced to one-fourth its original number by exploiting the inherent relationship between Q 3 and P 3 , along with the symmetric nature of the non-linear forces, which are cubic and depend only on the displacement variables. The details of these simplifications are now described.
For a given set of values for the Cs, the expression of the velocity constraint Q 3 , can be explicitly determined from the following relationship: which is the algebraic form of Eq. (10). Note that the non-linear forces f 1 and f 2 , in Eq. (16) are only dependent on the Cs, since all the non-linear forces are defined in terms of the displacement field only, Eq. (9). Otherwise, the relationship for Q 3 would be implicit. The velocity constraint Q 3 , would then have to be expanded as given in expression (13). The unknown coefficients, Ds, would need to be solved for simultaneously with the Cs, in an iterative manner. From expression (16), the velocity constraint Q 3 at each grid point (a i 1 , a j 2 ) is evaluated numerically at the following set of phase angles: Then, a two-dimensional fast Fourier transform (FFT) can be applied to these 2N 1 × 2N 2 discrete grid point values in order to obtain the two-dimensional Fourier coefficients corresponding to function Q 3 where theDs are the complex version of the D coefficients, as defined in Eq. (13), and N 1 , N 2 are set to be even. Note that 2N 1 × 2N 2 grid points in the phase domain are used to evaluate the N 1 × N 2 complex Fourier coefficients in Eq. (17), in order to reduce aliasing errors in the Fourier transform. Once the Fourier coefficients are obtained, the first-order derivatives, jQ 3 /j 1 and jQ 3 /j 2 , in Eq. (11) can be efficiently calculated by the twodimensional Inverse FFT (IFFT). The derivatives of Q 3 with respect to a 1 and a 2 can also be obtained by the backward finite difference scheme defined in Eq. (15). By this method, the total number of unknown coefficients in Eqs. (12) and (13) can be cut in half. Specifically, only the C coefficients need to be obtained. For this example system, the non-linear forces in Eq. (9) are only cubic. As a result, half of the Fourier series terms in Eq. (12) can be eliminated. In particular, in the two-dimensional Fourier series expansion of the constraint function P 3 , only the trigonometric basis functions whose combination orders are odd need to be included. With respect to each individual basis function, F l ( 1 )F m ( 2 ), in Eq. (12), the subscripts of the corresponding C coefficient can be used as a guide to determine whether or not the associated basis function shall be retained. If int[l/2] + int[m/2] is odd, the corresponding basis function is retained in the expansion, otherwise it is removed (where the operator int[a] denotes the maximum integer which is not larger than a). As a result, the total number unknown coefficients in Eq. (12) is halved again.
At this stage, the constraint function P 3 at a grid point is approximated by the reduced two-dimensional Fourier series and can be expressed in the following simplified form: where (a i 1 , a j 2 ) is a grid point in the amplitude domain of interest, N is the total number of expansion functions, T ( 1 , 2 ) is a simplified notation for the individual basis functions defined in Eq. (12), with odd harmonic combination order, and C (i,j ) is the corresponding unknown coefficient. Given an initial guess for the C s at all grid points in the amplitude domain, the complex Fourier coefficients,Ds, for the velocity constraint Q 3 are obtained from Eqs. (16) and (17). Then, the value of Q 3 and the corresponding partial derivatives are substituted into Eq. (11) along with the value of P 3 . The corresponding residual function R 3 , is defined at each grid point (a i 1 , a j 2 ) as follows: In order to minimize the residual function R 3 , a "weighted residuals" Galerkin method is used, in which the projection of the residual function onto the where T is the basis function in Eq. (18). Eq. (20) yields a set of non-linear algebraic equations in the unknown coefficients C . These are solved using a nonlinear solver found in the subroutine package NAG, which is based on Powell's hybrid method [20]. The user must, (i) provide an initial guess for the unknown coefficients in Eq. (18), and (ii) evaluate the "weighted residuals", the left-hand side of Eq. (20), by numerical integration.
The resulting four-dimensional invariant manifold cannot be visualized in three-dimensional space. However, we can show specific cross-sections of the manifold. In Fig. 2, the slave constraint relationship P 3 is depicted at the phase angles ( 1 , 2 )=(0, 0). The twodimensional amplitude domain (a 1 , a 2 ) over which the invariant manifold is numerically constructed is taken to be a 1 ∈ [0.01, 0.35] and a 2 ∈ [0.01, 0.35]. Using a convergence study, an 11-by-11 grid finite difference scheme is chosen to discretize this domain. The mesh sizes for h 1 and h 2 are 0.034 in the a 1 and a 2 directions. In Fig. 2, the invariant manifold looks smooth with this mesh scheme, which indicates that the chosen discretization is sufficient to capture the manifold geometry.
The domain defined by ( 1 , 2 ), where 1 ∈ [0, 2 ] and 2 ∈ [0, 2 ], is a two-dimensional torus. The invariant manifold at any grid point (a i 1 , a j 2 ) can be visualized in this torus domain. In Fig. 3(a), the phase angle 1 is defined by the angle from vector − − → Or 0 to vector − → Or. At phase angle 1 , a cross-section of the invariant manifold along the plane rOz is shown in Fig. 3(b). In this cross-section, the shape of the invariant manifold is shown by the vector − → , which indicates the magnitude of the manifold at phase angle 2 . For comparison, the nominal torus domain is also shown here as a circle with a dashed line, on which the magnitude of the invariant manifold is zero. For the solution shown in Fig. 2, the number of terms for the two-dimensional Fourier series in expression (12) is taken to be N 1 = N 2 = 12 for each grid point. Since the non-linear forces are cubic for this example system, the total number of basis functions in the two-dimensional phase space is reduced to N = 72 in expression (18). As a result, the total number of coefficients required for the construction of the constraint relationship P 3 , i.e. the C s in expression (18) for all of the 11 × 11 grid points in the two-dimensional am-plitude domain, is equal to 8712. For each unknown coefficient C , zero is used as the initial guess value to start Powell's hybrid solution method.
Once the constraint relationship P 3 is obtained for this internal resonance case, system motions on the invariant manifold can be captured by the reducedorder model, which involves the master coordinates only. As shown in Eq. (6), numerical time simulations can be carried out for these four first-order differential equations for given initial conditions, a 1 (0), 1 (0), a 2 (0), and 2 (0). In Fig. 4, time simulations for the master coordinates are shown using the reducedorder model with two DOFs and the original three-DOF model restricted to the invariant manifold. Based on the reduced-order model, Eq. (6), the responses of the amplitude-phase pairs are simulated. Then, the responses of the modal coordinates are obtained using the definition of the polar coordinate transformation. For comparison, the time responses for the master coordinates can also be acquired by the direct time simulation of the original system, Eq. (8), since the initial conditions for the slave coordinates 3 (0) and˙ 3 (0) can be obtained using the slave constraint function P 3 and Q 3 . According to the definition of the invariant manifold, any motion initiated on the manifold will remain on it for all the time, indicating that the simula-tion obtained from the original model should match the response from the reduced-order model if the invariant manifold has been constructed accurately enough to capture the actual geometry. In Fig. 4, the two categories of time responses are undistinguishable, indicating that the manifold geometry is accurate.
With the reduced-order model, time responses for the slave coordinates ( 3 ,˙ 3 ) can be obtained from the constraint relationships P 3 and Q 3 . In Figs. 5(a) and (b), these responses are compared to simulations based on the original system model. An excellent match between these two results is observed, which is further evidence to the accuracy of the invariant manifold. The time response of any physical coordinate, i.e. the displacement or velocity of any mass in Fig. 1, can be determined from the simulation of the reduced-order model, since the responses of all modal coordinates, both masters and slaves, are calculated. The displacement and velocity of mass M 3 , X 3 (t) andẊ 3 (t), are shown in Figs. 5(c) and (d). Again, excellent agreement is found between the results for the reducedorder and original system models.
This example was used to demonstrate the details of the process. A more substantial example is now considered.

The rotating beam system
Here the methodology is applied to a vibratory system of more practical interest, and with more degrees of freedom. A uniform rotating Euler-Bernoulli beam, shown in Fig. 6, is considered. This system has been studied by Pesheck et al. [18], who approximated the invariant manifold in the case of an internal resonance using asymptotic methods. Apiwattanalunggarn et al. [21] also studied this system and obtained the single-mode non-linear invariant manifold for large amplitude motions using a non-linear finite element approach. This rotating beam system can be considered as a highly simplified model of a helicopter rotor blade, in which the following effects are neglected: lead-lag motion, torsional motion, aerodynamic loading, and the weight of the blade. Even with this simplified model, it has been shown that typical discretization procedures for this system suffer from very slow modal convergence, since a comparatively large number of axial modes must be included in order to accurately describe the transverse bending motion of the beam [18,21]. Due to the non-linear axial/bending coupling effects, the resulting discrete models are computationally cumbersome, even for direct time simulations. Hence, a practical model order reduction technique would be very useful for the analysis of this system.
A detailed derivation of the PDE's governing the transverse bending, w(x, t), and the axial elongation, u(x, t), of this beam can be found in Ref. [18]. The derivation procedure is briefly described here. The potential energy, U, and kinetic energy, T, may be expressed as follows: where w(x, t) and u(x, t) are the transverse and axial displacements, respectively, () ,x denotes a partial derivative with respect to the spatial variable x, and() represents a time derivative. It should be noted that the standard linear curvature assumption is made in the energy expressions; this is done since the model still provides the non-linear coupling of interest, resulting in the slow modal convergence described above [18].
Hamilton's principle is used to develop the weak formulation for the equation of motion.
where () denotes the variation of a quantity. In order to obtain the discretized version of the equations of motion in the standard form (1), the linearized partial differential equations off the non-linear rotating beam system are given here where u d is the dynamic component of the axial elongation, u(x, t), defined as The static part, u s (x), is the static elongation of the beam due to rotation when the transverse deflection is zero. Linear mode shapes corresponding to Eqs. (24) and (25) can be obtained using a Rayleigh-Ritz procedure. Once these modes are determined, the solutions to the non-linear system (23) are sought in the form of an expansion as whereŪ i andW i are the linear modes corresponding to Eqs. (24) and (25), and the integers N and N denote the number of axial and transverse linear modes used, respectively. These expansions are substituted into the weak formulation, Eq. (23), and the discretized non-linear equations of motion are obtained as follows: where ,i and ,i are the linear modal frequencies associated with the ith modes in the axial and transverse directions, respectively. These equations of motion are non-linearly coupled because the quadratic non-linear forces corresponding to the axial motion f ,i , depend on the transverse motion j , while the quadratic and cubic non-linear forces in the transverse direction f ,i , are dependent on both axial and transverse motions, j and k .
The convergence of this model has been thoroughly investigated in Ref. [18]. It has been found that at least an 18-DOF model, with N = N = 9, must be used to accurately capture the periodic response in the vicinity of the first non-linear mode for an energy level corresponding to a transverse deflection amplitude of about 0.1 m at the beam tip, for a 9.0 m beam. This 18-DOF discretized model is used here as the reference model. Based on this model, the invariant manifold is constructed. Consequently, a reduced-order model can be obtained for the representation of the dynamics on the invariant manifold.
As in the case studied in Ref. [18], the parameters of the uniform rotating beam are set as follows: L = 9 m, m=10 kg/m, EI=3.99×10 5 N m 2 , EA=2.23×10 8 N, = 23.85 rad/s, and h = 0.5 m. Under these conditions, a three-to-one internal resonance occurs between the first two transverse modes, ,2 ≈ 3 ,1 . The master coordinates are chosen as the state variable pairs, ( 1 ,˙ 1 ) and ( 2 ,˙ 2 ). Polar coordinate transformations, defined in Eq. (3), are applied to these two pairs of state variables, resulting in two amplitudephase pairs as the transformed master coordinates, (a 1 , 1 ) and (a 2 , 2 ). All the remaining DOFs, including 7 transverse deflection modes and 9 axial modes, form the slave coordinates, which are constrained as follows: Thus, there are a total of 16 pairs of constraint relationships, Eq. (30), that need to be solved. The governing PDEs for these constraint functions are given in Eq. (7). The invariant manifold is solved for numerically in the following four-dimensional domain: 1 , a 2 , 1 , 2  where the amplitude range is carefully chosen so that the non-linear effect in the system is sufficiently strong and the invariant manifold obtained from the asymptotic expansion method in Ref. [18] becomes inaccurate over this domain.
In the solution procedure, the constraint relationships for the velocities Q i s, in Eq. (30) are not solved for, due to the fact that the velocity constraint is the time derivative of the corresponding displacement constraint. The details of the reduction have been given in Section 3, Eqs. (16) and (17), for the construction of the invariant manifold of the 3-DOF example system. Hence, only the displacement constraint relationships, P i s, need to be solved.
For the two-dimensional phase domain, ( 1 , 2 ), the two-dimensional Fourier series, defined in Eq. (12), is utilized for the expansion functions for the displacement constraints. Because both quadratic and cubic non-linear terms exist in this system, the expansion cannot be further simplified. For the twodimensional amplitude domain, (a 1 , a 2 ), the finite difference discretization scheme, which was used in the 3-DOF example system, cannot be utilized here due to limitations in computational capacity. The numerical difficulty is clearly shown by the following case: Let us divide the two-dimensional amplitude domain into 8/8 grid points, and set the number of terms in the Fourier expansion as N 1 = N 2 = 8 at each grid point. Then, the total number of unknown coefficients is 4096 for each displacement constraint relationship P i . With 16 slave constraints in Eq. (30), the final number of unknowns is 65,536. It is inefficient to solve for the invariant manifold with such a large number of unknowns.
A strategy to overcome this numerical difficulty is to discretize the two-dimensional amplitude domain into small elements, and then utilize low-order polynomials as expansion functions in the discretized elements. For this example system, the amplitude domain, {(a 1 , a 2 ) | a 1 ∈ [0.01, 0.75], a 2 ∈ [0.01, 0.4]}, is evenly divided into 7/7 equal-sized patches. The width of each patch is 0.1057 along the a 1 direction, and 0.05571 along the a 2 direction. The displacement constraint relationships, the P i s in Eq. (30), are then expanded in each discretized four-dimensional element where T j (a 1 ) are piecewise linear functions defined in the amplitude segments, a 1 ∈ [a low 1 , a up 1 ], as follows: The definition of the piecewise linear functions, T k (a 2 ), is the same as for T j (a 1 ), where the lower and upper limits of the amplitude segments are set as a 2 ∈ [a low 2 , a up 2 ]. The Fourier terms, F l ( 1 ) and F m ( 2 ), are defined in Eq. (14).
In each element, given by Eq. (31), the deduction of the velocity constraint (Q i ) from the corresponding displacement constraint (P i ), and the evaluation of the residue function, R i , are again given by Eqs. (16), (17), and (19). It should be noted that in the present example, the numerical values of the velocity constraints (Q i ) and the residue functions (R i ) are now evaluated at the Gaussian quadrature points for polynomials in the discretized two-dimensional amplitude domain. A three-by-three-point Gaussian quadrature formula is sufficient in the region, {(a 1 , a 2 ) | a 1 ∈ [a low 1 , a We set the number of the Fourier terms in expansion (32) to be N 1 = N 2 = 8. As a result, the total number of the unknown coefficients, the Cs in expansion (32), is equal to 4096 for each element, for each of the 16 slave constraint relationships. Note that the total number of unknown quantities resulting from the finite difference discretization scheme in the whole amplitude domain is equal to 65, 536. Thus, it is seen that the computational cost is tremendously reduced for the non-linear solver, since the invariant manifold, defined by Eq. (30), is now solved for independently in each four-dimensional discretized element.
The initial values used in the numerical solution of the Cs for each discretized element are determined as follows. For the first element, which has a twodimensional amplitude domain given by   {(a 1 , a 2 ) | a 1 ∈ [0.01, 0.1157], a 2 ∈ [0.01, 0.06571]}, zeros are good initial values, due to the fact that the non-linearities are weak near the origin. Then, for subsequent elements, which have incremental values in the a 1 or a 2 directions, the expansion coefficients obtained from the preceding element are used as the initial values. Once the results for all discretized elements are obtained, the expansion coefficients from contiguous elements are averaged at their interface. The resulting solution for the invariant manifold is stitched together to cover the entire domain of interest. With the obtained invariant manifold, the reference model with 36 states can be reduced to a 4-state model for this internally resonant case.
A cross-section of the invariant manifold is shown in Fig. 7. The slave constraint relationship for the third transverse deflection mode, P 3 in Eq. (30), is depicted at the phase angles ( 1 , 2 ) = (0, 0). The amplitude domain, a 1 ∈ [0.01, 0.75] and a 2 ∈ [0.01, 0.4], is evenly divided into 7/7 patches, and the invariant manifold appears smooth with this mesh. Note that the invariant manifold defined in Eq. (30) is the ensemble of the displacement and velocity constraint relationships for all 16 slave coordinates, and Fig. 7 represents simply the cross-section of one slave coordinate among the 16.
Time responses for the displacements of the master and slave coordinates are shown and compared in Figs. 8 and 9 using three different simulation approaches: (i) direct time simulations based on the 36-state reference model, with initial conditions that satisfy the constraint relationships; (ii) time simulations for the master coordinates using the 4-state reduced-order model, along with the reconstruction of the slave coordinate responses using the constraint functions; and (iii) simulations based on the reducedorder model obtained by the asymptotic expansion method described in Pesheck et al. [18], wherein the invariant manifold and the corresponding reduced- order model were generated using asymptotic series expansions. It is seen that simulations obtained from the reduced-order model match the reference model results precisely, while the results from the asymptotic method depart from the reference response rather quickly as time progresses. This is not surprising, because the combination orders of the multi-dimensional polynomials used in the asymptotic expansion method are limited to three [18]. Here, the combination order of the trigonometric functions in Eq. (32) can be as large as eight. Consequently, the invariant manifolds constructed here are more accurate than the manifolds obtained in [18], and the simulations will match more closely over a longer period of time, especially at larger amplitudes. The better accuracy of the reduced-order model can also be verified by observing simulations of the transverse displacement of the tip of the beam, as shown in Fig. 10. In Figs. 8-10, the accuracy of the reduced-order model has been verified by comparisons of the simulated time responses. The 4-state reduced-order model can then be utilized to investigate the dynamic behavior of this system, which arises from the existence of the internal resonance. Amplitude modulation of the responses for the two master coordinates, a 1 (t) and a 2 (t), is demonstrated in Fig. 11. Note that there exists a continuous exchange of energy between the two modes. Within the first second, the time period of the energy exchange can be approximately determined as 0.12 s. Similar properties can also be found in Fig. 12, where motions are simulated over a long time period. The energy exchange shown in Fig. 12 occurs at a much slower time scale, with a period of about 9 s.

Conclusions
The following conclusions can be drawn from this study: (i) Multi-NNMs can be effectively generated by the invariant manifold approach. A systematic solution methodology for the invariant manifold has been proposed, which uses the polar form of the master coordinates. Four-dimensional invariant manifolds have been successfully constructed for the 3-DOF exam- ple system and for the rotating beam system, using a combination of finite difference or finite element discretization schemes in the amplitude domain and twodimensional Fourier series expansions in the phase domain. (ii) A reduced-order model can be generated once the multi-NNM is obtained, and motions on the invariant manifold can be accurately captured by this model. The precision of the reduced-order model is controlled by the numerical parameters used in the solution procedure. (iii) Although only quadratic and cubic order non-linear forces were considered in the systems considered, the construction method can be extended to systems with more complicated non-linear forces, and this is relatively straightforward if the nonlinear forces do not depend on velocities. Otherwise, the numerical solution algorithm will need to simultaneously solve for the displacement (P i ) and velocity (Q i ) constraint functions, and will involve many more unknown coefficients. (iv) For complicated dynamic systems, such as more realistic rotating blade models that include lead-lag and torsional motions, gyroscopic effects and damping forces must be considered in the linear order model. The multi-mode invariant manifold approach can be extended to such systems. However, complex linear modal analysis must be used to obtain a revised form of the master and slave coordinates [17].