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

Dongying Jiang, Christophe Pierre, Steven Shaw

► To cite this version:

Dongying Jiang, Christophe Pierre, Steven Shaw. The construction of non-linear normal modes for systems with internal resonance. International Journal of Non-Linear Mechanics, 2005, 40 (5), 10.1016/j.ijnonlinmec.2004.08.010. hal-01350807

HAL Id: hal-01350807

https://hal.science/hal-01350807

Submitted on 1 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The construction of non-linear normal modes for systems with internal resonance

D. Jiangª, C. Pierreᵇ, S.W. Shawᶜ

ªMKP Structural Design Associates, Inc., 3003 Washtenaw Ave Suite 1-E, Ann Arbor, MI 48104, USA

ᵇDepartment of Mechanical Engineering, University of Michigan, 3112 G. G. Brown Bldg., Ann Arbor, MI 48109-2125, USA

ᶜ2328H Engineering, Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824-1226, USA

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.

# 1. 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–6], and Vakakais [7–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
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 non-linear 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 Pescheck 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 $\bar{a}$ and *M* phases $\phi$. Since the master coordinates are represented by more than one amplitude-phase 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 non-linear 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.

## 2. 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 assumed-mode 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:

$$ \ddot{\eta}_i + \omega_i^2 \eta_i = f_i(\eta_j), \quad i, j = 1-\dots n, \quad (1) $$

where the $\eta_i$s are the linear modal coordinates, the $\omega_i$s are the corresponding linear modal frequencies, and $f_i(\eta_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(\eta_j)$s, are assumed to be independent of the linear modal velocities, $\dot{\eta}_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, multi-mode 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, $(\eta_k, \dot{\eta}_k) \in 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

$$ \begin{aligned} \eta_i &= X_i(\eta_k, \dot{\eta}_k) \\ \dot{\eta}_i &= Y_i(\eta_k, \dot{\eta}_k), \quad \text{for } i \notin S_M, k \in S_M. \end{aligned} \quad (2) $$

This form has been utilized by Pescheck 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, $(\eta_k, \dot{\eta}_k) \in S_M$, a polar coordinate transformation is applied,

$$ \begin{aligned} \eta_k &= a_k \cos(\phi_k) \\ \dot{\eta}_k &= -a_k \omega_k \sin(\phi_k), \quad \text{for } k \in S_M, \end{aligned} \quad (3) $$

where $\omega_k$ is the kth linear modal frequency, and $(a_k, \phi_k)$ are the master coordinates in amplitude-phase form. Then, all the slave coordinates are expressed as functions of these amplitude-phase master coordinates,

$$ \begin{aligned} \eta_i &= P_i(a_k, \phi_k) \\ \dot{\eta}_i &= Q_i(a_k, \phi_k), \quad \text{for } i \notin S_M, k \in S_M, \end{aligned} \quad (4) $$

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, $(\eta_i, \dot{\eta}_i)$, we have

$$
\begin{aligned}
Q_i &= \dot{P}_i \\
&= \sum_{k \in S_M} \left( \frac{\partial P_i}{\partial a_k} \dot{a}_k + \frac{\partial P_i}{\partial \phi_k} \dot{\phi}_k \right) \\
\dot{Q}_i &= \sum_{k \in S_M} \left( \frac{\partial Q_i}{\partial a_k} \dot{a}_k + \frac{\partial Q_i}{\partial \phi_k} \dot{\phi}_k \right), \quad \text{for } i \notin S_M, \quad (5) \\
&= -\omega_i^2 P_i + f_i
\end{aligned}
$$

where $\omega_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, $\dot{a}_k$ and $\dot{\phi}_k$ in Eq. (5) are replaced using the governing equation of motion for each pair of master coordinates, which are

$$
\begin{aligned}
\dot{a}_k &= \frac{-f_k(\eta_j)}{\omega_k} \sin \phi_k \\
\dot{\phi}_k &= \omega_k - \frac{f_k(\eta_j)}{\omega_k a_k} \cos \phi_k,
\end{aligned}
\quad \text{for } k \in S_M, \quad j=1-n. \tag{6}
$$

Substituting Eq. (6) into Eq. (5), the governing PDEs
for the invariant manifold are found to be

$$
\begin{aligned}
Q_i &= \sum_{k \in S_M} \left[ \frac{\partial P_i}{\partial a_k} \left( -\frac{f_k}{\omega_k} \sin \phi_k \right) \right. \\
    &\qquad \left. + \frac{\partial P_i}{\partial \phi_k} \left( \omega_k - \frac{f_k}{\omega_k a_k} \cos \phi_k \right) \right] \\
\sum_{k \in S_M} &\left[ \frac{\partial Q_i}{\partial a_k} \left( -\frac{f_k}{\omega_k} \sin \phi_k \right) \right. \\
&\qquad \left. + \frac{\partial Q_i}{\partial \phi_k} \left( \omega_k - \frac{f_k}{\omega_k a_k} \cos \phi_k \right) \right] = -\omega_i^2 P_i + f_i \\
\text{for } i \notin S_M.
\end{aligned}
\tag{7}
$$

In Eq. (7), there are a total (n - M) pairs of equa-
tions 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ᵢ and Qᵢ
for i ≠ Sₘ, can be substituted into the M pairs of or-
dinary differential equations governing the dynamics
of the master coordinates, Eq. (6). As a result, the re-
sponse of the original system restricted to the invariant
manifold is captured by this 2M-DOF reduced-order
model. The procedure for obtaining the numerical so-
lution for the invariant manifolds is described in the
context of a simple example, and then applied to a
more substantial problem.

3. A three-DOF example system

Fig. 1 depicts a three-DOF mass-spring system with two cubic non-linear springs of coefficients $\gamma_1$ and $\gamma_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. $\omega_2 \approx 3\omega_1$. Consequently, a three-to-one internal resonance occurs between the first and second linear modes.

Using linear modal coordinates, the system can be
transformed to the standard form shown in Eq. (1)

$$
\dot{\eta}_i + \omega_i^2 \eta_i = f_i(\eta_j) \quad \text{for } i, j = 1, 2, 3, \quad (8)
$$

where the modal coordinates $\eta_i$ are defined by the
linear modal transformation:

$$
\begin{Bmatrix} x_1 \\ x_2 \\ x_3 \end{Bmatrix} = \begin{bmatrix} 0.5870 & -0.5147 & 0.6249 \\ 0.6039 & -0.2357 & -0.7614 \\ 0.5392 & 0.8243 & 0.1725 \end{bmatrix} \begin{Bmatrix} \eta_1 \\ \eta_2 \\ \eta_3 \end{Bmatrix}.
$$

The linear modal frequencies are $\omega_1 = 0.5972 \text{ rad/s}$,
$\omega_2 = 1.792 \text{ rad/s}$, and $\omega_3 = 3.405 \text{ rad/s}$, and the non-
linear forces are

$$
\begin{align}
f_1 &= 0.587(-0.587\eta_1 + 0.5147\eta_2 - 0.6249\eta_3)^3 \nonumber \\
    &\quad - 0.8088(0.5392\eta_1 + 0.8243\eta_2 + 0.1725\eta_3)^3, \nonumber \\
f_2 &= -0.5147(-0.587\eta_1 + 0.5147\eta_2 - 0.6249\eta_3)^3 \nonumber \\
    &\quad - 1.2365(0.5392\eta_1 + 0.8243\eta_2 + 0.1725\eta_3)^3, \nonumber \\
f_3 &= 0.6249(-0.587\eta_1 + 0.5147\eta_2 - 0.6249\eta_3)^3 \nonumber \\
    &\quad - 0.2588(0.5392\eta_1 + 0.8243\eta_2 \nonumber \\
    &\qquad + 0.1725\eta_3)^3. \tag{9}
\end{align}
$$

In the presence of an internal resonance between the
first two modes, the master coordinates are chosen as
the state variable pairs $(\eta_1, \dot{\eta}_1)$ and $(\eta_2, \dot{\eta}_2)$. The cor-
responding master coordinate index set is $S_M = \{1, 2\}$.
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, \phi_1)$ and $(a_2, \phi_2)$.

As defined in Eq. (4), the constraint relationships
for the slave coordinates are

$$
\eta_3 = P_3(a_1, a_2, \phi_1, \phi_2),
$$

$$
\dot{\eta}_3 = Q_3(a_1, a_2, \phi_1, \phi_2).
$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic diagram of the mass-spring system: masses $M_1 = M_2 = M_3 = 1$ kg; spring stiffnesses $K_1 = 0.5$ N/m, $K_2 = 5$ N/m, $K_3 = 2.03$ N/m, $K_4 = 0.6$ N/m; non-linear spring forces $\gamma_1 = 1.0 \times X_1(t)^3$ N, $\gamma_2 = 1.5 \times X_3(t)^3$ N, $X_1, X_2$, and $X_3$ denote the displacements of masses $M_1, M_2$, and $M_3$, respectively.</figcaption></figure>

The governing PDEs for the invariant manifold are given in Eq. (7) and are rewritten here as

$$Q_3 = \frac{\partial P_3}{\partial a_1} \left( \frac{-f_1 \sin \phi_1}{\omega_1} \right) + \frac{\partial P_3}{\partial \phi_1} \left( \omega_1 - \frac{f_1 \cos \phi_1}{\omega_1 a_1} \right) + \frac{\partial P_3}{\partial a_2} \left( \frac{-f_2 \sin \phi_2}{\omega_2} \right) + \frac{\partial P_3}{\partial \phi_2} \left( \omega_2 - \frac{f_2 \cos \phi_2}{\omega_2 a_2} \right), \quad (10)$$

$$\frac{\partial Q_3}{\partial a_1} \left( \frac{-f_1 \sin \phi_1}{\omega_1} \right) + \frac{\partial Q_3}{\partial \phi_1} \left( \omega_1 - \frac{f_1 \cos \phi_1}{\omega_1 a_1} \right) + \frac{\partial Q_3}{\partial a_2} \left( \frac{-f_2 \sin \phi_2}{\omega_2} \right) + \frac{\partial Q_3}{\partial \phi_2} \left( \omega_2 - \frac{f_2 \cos \phi_2}{\omega_2 a_2} \right) = -\omega_3^2 P_3 + f_3. \quad (11)$$

A numerical solution scheme is given here to approximate the unknown constraint relationships, $P_3(a_1, a_2, \phi_1, \phi_2)$ and $Q_3(a_1, a_2, \phi_1, \phi_2)$, which define the geometry of the four-dimensional invariant manifold. In the two-dimensional phase region, defined by

$$\{(\phi_1, \phi_2) | \phi_1 \in [0, 2\pi], \phi_2 \in [0, 2\pi]\},$$

the constraint relationships $P_3$ and $Q_3$, are periodic in both $\phi_1$ and $\phi_2$. Hence, they can be efficiently approximated by two-dimensional Fourier series. In the two-dimensional amplitude domain, defined by

$$\{(a_1, a_2) | 0 < a_1 < a_1 \text{max}, 0 < a_2 < a_2 \text{max}\},$$

where the upper limits $a_1 \text{max}$ and $a_1 \text{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 \text{ 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

$$P_3(a_1^i, a_2^j, \phi_1, \phi_2) \approx \sum_{l=1}^{N_{\phi_1}} \sum_{m=1}^{N_{\phi_2}} C_{lm}^{(i,j)} F_l(\phi_1) F_m(\phi_2), \quad (12)$$

$$Q_3(a_1^i, a_2^j, \phi_1, \phi_2) \approx \sum_{l=1}^{N_{\phi_1}} \sum_{m=1}^{N_{\phi_2}} D_{lm}^{(i,j)} F_l(\phi_1) F_m(\phi_2), \quad (13)$$

where $(a_1^i, a_2^j)$ 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(\phi_1)$ and $F_m(\phi_2)$, are defined as

$$F_l(\phi) = \begin{cases} \cos\left(\frac{l-1}{2}\phi\right), & l \text{ is odd,} \\ \sin\left(\frac{l}{2}\phi\right), & l \text{ is even} \end{cases} \quad (14)$$
and $N_{\phi_1}$ and $N_{\phi_2}$ are the number of terms of the Fourier expansions in $\phi_1$ and $\phi_2$, respectively. As can be seen in Eqs. (12)–(14), the total number of the unknown quantities are determined by $N_{\phi_1}$ and $N_{\phi_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

$$ \begin{aligned} & \frac{\partial P_3}{\partial a_1}(a_1^i, a_2^j, \phi_1, \phi_2) \\ & \approx \frac{1}{h_1} [P_3(a_1^i, a_2^j, \phi_1, \phi_2) \\ & \qquad - P_3(a_1^{i-1}, a_2^j, \phi_1, \phi_2)] + O(h_1), \\ & \frac{\partial P_3}{\partial a_2}(a_1^i, a_2^j, \phi_1, \phi_2) \\ & \approx \frac{1}{h_2} [P_3(a_1^i, a_2^j, \phi_1, \phi_2) \\ & \qquad - P_3(a_1^i, a_2^{j-1}, \phi_1, \phi_2)] \\ & \qquad + O(h_2), \end{aligned} \quad (15) $$

where $h_1 = a_1^i - a_1^{i-1}$, $h_2 = a_2^j - a_2^{j-1}$ are the distances between adjacent grid points along the $a_1$ and $a_2$ directions, respectively. The function $O(\cdot)$ denotes that the errors in these approximations are orders-of-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, $\partial Q_3/\partial a_1$ and $\partial Q_3/\partial a_2$, is determined using a similar scheme. Along the $\phi_1$ and $\phi_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 reduced

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:

$$ \begin{aligned} & Q_3(a_1^i, a_2^j, \phi_1, \phi_2; C) \\ & \approx \frac{\partial P_3}{\partial a_1}(a_1^i, a_2^j, \phi_1, \phi_2; C) \\ & \qquad \times \left[ -f_1(a_1^i, a_2^j, \phi_1, \phi_2; C) \frac{\sin \phi_1}{\omega_1} \right] \\ & \qquad + \frac{\partial P_3}{\partial \phi_1}(a_1^i, a_2^j, \phi_1, \phi_2; C) \\ & \qquad \times \left[ \omega_1 - f_1(a_1^i, a_2^j, \phi_1, \phi_2; C) \frac{\cos \phi_1}{\omega_1 a_1} \right] \\ & \qquad + \frac{\partial P_3}{\partial a_2}(a_1^i, a_2^j, \phi_1, \phi_2; C) \\ & \qquad \times \left[ -f_2(a_1^i, a_2^j, \phi_1, \phi_2; C) \frac{\sin \phi_2}{\omega_2} \right] \\ & \qquad + \frac{\partial P_3}{\partial \phi_2}(a_1^i, a_2^j, \phi_1, \phi_2; C) \\ & \qquad \times \left[ \omega_2 - f_2(a_1^i, a_2^j, \phi_1, \phi_2; C) \frac{\cos \phi_2}{\omega_2 a_2} \right], \end{aligned} \quad (16) $$

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_1^i, a_2^j)$ is evaluated numerically at the following set of phase angles:

$$ \begin{array}{l} Q_3(a_1^i, a_2^j, \bar{\phi}_1^I, \bar{\phi}_2^J) \quad \text{where} \quad \left\{ \begin{array}{l} \bar{\phi}_1^I = I\pi/N_{\phi_1}, \\ \bar{\phi}_2^J = J\pi/N_{\phi_2}, \end{array} \right. \\ \text{for } I = 1 \dots 2N_{\phi_1}, \quad J = 1 \dots 2N_{\phi_2}. \end{array} $$
Then, a two-dimensional fast Fourier transform (FFT) can be applied to these $2N_{\phi_1} \times 2N_{\phi_2}$ discrete grid point values in order to obtain the two-dimensional Fourier coefficients corresponding to function $Q_3$

$$Q_3(a_1^i, a_2^j, \phi_1, \phi_2) = \sum_{m=-\frac{N_{\phi_1}}{2}}^{\frac{N_{\phi_1}}{2}} \sum_{n=-\frac{N_{\phi_2}}{2}}^{\frac{N_{\phi_2}}{2}} \times \tilde{D}_{mn}^{(i,j)} e^{\sqrt{-1}m\phi_1} e^{\sqrt{-1}n\phi_2}, \quad (17)$$

where the $\tilde{D}$s are the complex version of the D coefficients, as defined in Eq. (13), and $N_{\phi_1}, N_{\phi_2}$ are set to be even. Note that $2N_{\phi_1} \times 2N_{\phi_2}$ grid points in the phase domain are used to evaluate the $N_{\phi_1} \times N_{\phi_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, $\partial Q_3/\partial \phi_1$ and $\partial Q_3/\partial \phi_2$, in Eq. (11) can be efficiently calculated by the two-dimensional 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(\phi_1)F_m(\phi_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 $\text{int}[l/2] + \text{int}[m/2]$ is odd, the corresponding basis function is retained in the expansion, otherwise it is removed (where the operator $\text{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:

$$P_3(a_1^i, a_2^j, \phi_1, \phi_2) \approx \sum_{\alpha=1}^{N_{\alpha}} C_{\alpha}^{(i,j)} T_{\alpha}(\phi_1, \phi_2), \quad (18)$$

where $(a_1^i, a_2^j)$ is a grid point in the amplitude domain of interest, $N_{\alpha}$ is the total number of expansion functions, $T_{\alpha}(\phi_1, \phi_2)$ is a simplified notation for the individual basis functions defined in Eq. (12), with odd harmonic combination order, and $C_{\alpha}^{(i,j)}$ is the corresponding unknown coefficient. Given an initial guess for the $C_{\alpha}$s at all grid points in the amplitude domain, the complex Fourier coefficients, $\tilde{D}$s, 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_1^i, a_2^j)$ as follows:

$$R_3(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) = \frac{\partial Q_3}{\partial a_1}(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \times \left[ -f_1(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \frac{\sin \phi_1}{\omega_1} \right] + \frac{\partial Q_3}{\partial \phi_1}(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \times \left[ \omega_1 - f_1(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \frac{\cos \phi_1}{\omega_1 a_1} \right] + \frac{\partial Q_3}{\partial a_2}(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \times \left[ -f_2(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \frac{\sin \phi_2}{\omega_2} \right] + \frac{\partial Q_3}{\partial \phi_2}(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \times \left[ \omega_2 - f_2(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) \frac{\cos \phi_2}{\omega_2 a_2} \right] + \omega_3^2 P_3(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha) - f_3(a_1^i, a_2^j, \phi_1, \phi_2; C_\alpha). \quad (19)$$

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
<figure><img src="image_3.png" /><figcaption>Fig. 2. A section of one component of the invariant manifold for the system shown in Fig. 1. The modal displacement $\eta_3$ is shown at phase angles $(\phi_1, \phi_2) = (0, 0)$.</figcaption></figure>

basis function is required to vanish

$$ \int_0^{2\pi} \int_0^{2\pi} [T_{\alpha}(\phi_1, \phi_2) \times R_3(a_1^i, a_2^j, \phi_1, \phi_2; C_{\alpha})] d\phi_1 d\phi_2 = 0, \\ \text{for } \alpha = 1 \dots N_{\alpha}, \quad \text{at } \forall (a_1^i, a_2^j), \qquad (20) $$

where $T_{\alpha}$ is the basis function in Eq. (18). Eq. (20) yields a set of non-linear algebraic equations in the unknown coefficients $C_{\alpha}$. These are solved using a non-linear 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 $(\phi_1, \phi_2)=(0, 0)$. The two-dimensional amplitude domain $(a_1, a_2)$ over which the invariant manifold is numerically constructed is taken to be $a_1 \in [0.01, 0.35]$ and $a_2 \in [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

<figure><img src="image_4.png" /><figcaption>Fig. 3. Diagram of another section of a component of the invariant manifold for the system shown in Fig. 1. The modal displacement $\eta_3$ is shown at amplitude $(a_1, a_2) = (0.35, 0.35)$. (a) An illustration for the invariant manifold defined in the torus domain $\{(\phi_1, \phi_2) | \phi_1 \in [0, 2\pi], \phi_2 \in [0, 2\pi]\}$; (b) A cross-sectional view of the invariant manifold at phase angle $\phi_1 = 2/3\pi$.</figcaption></figure>

chosen discretization is sufficient to capture the manifold geometry.

The domain defined by $(\phi_1, \phi_2)$, where $\phi_1 \in [0, 2\pi]$ and $\phi_2 \in [0, 2\pi]$, is a two-dimensional torus. The invariant manifold at any grid point $(a_1^i, a_2^j)$ can be visualized in this torus domain. In Fig. 3(a), the phase angle $\phi_1$ is defined by the angle from vector $\vec{Or}_0$ to vector $\vec{Or}$. At phase angle $\phi_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 $\vec{p}$, which
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 4. Comparison of the time responses of the master coordinates for the system shown in Fig. 1, for initial conditions $a_1(0) = 0.28$, $\phi_1(0) = 0.0$, $a_2(0) = 0.28$, $\phi_2(0) = 0.0$: (—) time simulation of the original system model; (---) simulation of the reduced-order model. Plot (a) shows the time response of $\eta_1$ versus $t$; (b) $\dot{\eta}_1(t)$; (c) $\eta_2(t)$; (d) $\dot{\eta}_2(t)$.</figcaption></figure>

indicates the magnitude of the manifold at phase angle $\phi_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_{\phi_1} = N_{\phi_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_{\alpha}=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_{\alpha 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_{\alpha}$, 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 reduced-order 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)$, $\phi_1(0)$, $a_2(0)$, and $\phi_2(0)$. In Fig. 4, time simulations for the master coordinates are shown using the reduced-order model with two DOFs and the original three-DOF model restricted to the invariant manifold. Based
<figure><img src="image_9.png" /><figcaption>Fig. 5. Comparison of the time responses of the slave coordinates and mass $M_3$ for the system shown in Fig. 1, for initial conditions deduced from the initial values of the master coordinates: $a_1(0) = 0.28$, $\phi_1(0) = 0.0$, $a_2(0) = 0.28$, $\phi_2(0) = 0.0$: (—) time simulation of the original system model; (---) simulation based on the reduced-order model. Plot (a) shows the time response of $\eta_3$ versus $t$; (b) $\dot{\eta}_3(t)$; (c) $X_3(t)$; (d) $\dot{X}_3(t)$.</figcaption></figure>

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 $\eta_3(0)$ and $\dot{\eta}_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 simulation

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 ($\eta_3, \dot{\eta}_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 $\dot{X}_3(t)$, are shown in Figs. 5(c) and (d). Again, excellent agreement is found between the results for the reduced-order and original system models.

This example was used to demonstrate the details of the process. A more substantial example is now considered.

## 4. 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 Pescheck 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:

$$T = \frac{1}{2} \int_{0}^{L} m(\dot{u}^2 + \dot{w}^2) + m\Omega^2(h + x + u)^2 dx, \quad (21)$$

$$U = \frac{1}{2} \int_{0}^{L} EI(w_{,xx})^2 + EA \left( u_{,x} + \frac{1}{2}(w_{,x})^2 \right)^2 dx, \quad (22)$$

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 $\dot{)}$ 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.

$$\int_{t_1}^{t_2} \int_{0}^{L} \left\{ \begin{aligned} & [-m\ddot{w} - EIw_{,xxxx}]\delta w \\ & - \left[ EA \left( u_{,x} + \frac{1}{2}(w_{,x})^2 \right) w_{,x} \right] \delta w_{,x} \\ & + [-m\ddot{u} + m\Omega^2(x + h + u) + EAu_{,xx}]\delta u \\ & - \left[ EA \frac{1}{2}(w_{,x})^2 \right] \delta u_{,x} \end{aligned} \right\} dx dt = 0, \quad (23)$$

where $\delta()$ 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

$$m\ddot{u}_d - m\Omega^2 u_d - EAu_{d,xx} = 0, \quad (24)$$

$$m\ddot{w} + EIw_{,xxxx} - EA(u_{s,x}w_{,xx} + u_{s,xx}w_{,x}) = 0, \quad (25)$$

where $u_d$ is the dynamic component of the axial elongation, $u(x, t)$, defined as

$$u_d(x, t) = u(x, t) - u_s(x). \quad (26)$$

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
<figure><img src="image_10.png" /><figcaption>Fig. 6. Schematic diagram of the rotating beam system. Rotating speed $\Omega$ is constant. The uniform beam has the following material and geometrical parameters: Young's modulus $E$, cross-sectional area $A$, second moment of cross-sectional area $I$, length $L$, and mass per unit length $m$. The hub radius is $h$. (Abridged from Ref. [18].)</figcaption></figure>

to the non-linear system (23) are sought in the form of an expansion as

$$ 
\begin{aligned}
u_d(x, t) &= \sum_{i=1}^{N_{\eta}} \eta_i(t) \bar{U}_i(x), \\
w(x, t) &= \sum_{i=1}^{N_{\xi}} \xi_i(t) \bar{W}_i(x),
\end{aligned} 
\quad (27) $$

where $\bar{U}_i$ and $\bar{W}_i$ are the linear modes corresponding to Eqs. (24) and (25), and the integers $N_{\eta}$ and $N_{\xi}$ 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:

$$ 
\begin{aligned}
\ddot{\eta}_i + \omega_{\eta,i}^2 \eta_i &= f_{\eta,i}(\xi_j) \\
\text{for } i = 1-N_{\eta}, \quad j = 1-N_{\xi},
\end{aligned} 
\quad (28) $$

$$ 
\begin{aligned}
\ddot{\xi}_i + \omega_{\xi,i}^2 \xi_i &= f_{\xi,i}(\xi_j, \eta_k) \\
\text{for } i, j = 1-N_{\xi}, \quad k = 1-N_{\eta},
\end{aligned} 
\quad (29) $$

where $\omega_{\eta,i}$ and $\omega_{\xi,i}$ are the linear modal frequencies associated with the $i$th 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_{\eta,i}$, depend on the transverse motion $\xi_j$, while the quadratic and cubic non-linear forces in the transverse direction $f_{\xi,i}$, are dependent on both axial and transverse motions, $\xi_j$ and $\eta_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_{\eta} = N_{\xi} = 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 \times 10^5$ N m$^2$, $EA=2.23 \times 10^8$ N, $\Omega = 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, $\omega_{\xi,2} \approx 3\omega_{\xi,1}$. The master coordinates are chosen as the state variable pairs, $(\xi_1, \dot{\xi}_1)$ and $(\xi_2, \dot{\xi}_2)$. Polar coordinate transformations, defined in Eq. (3), are applied to these two pairs of state variables, resulting in two amplitude-phase pairs as the transformed master coordinates, $(a_1, \phi_1)$ and $(a_2, \phi_2)$. All the remaining DOFs, including 7 transverse deflection modes and 9 axial modes, form the slave coordinates, which are constrained as follows:

$$ 
\begin{aligned}
\dot{\xi}_i &= P_i(a_1, a_2, \phi_1, \phi_2), \\
\dot{\eta}_i &= Q_i(a_1, a_2, \phi_1, \phi_2), \quad i = 3 \dots 9, \\
\dot{\eta}_i &= P_{i+9}(a_1, a_2, \phi_1, \phi_2), \\
\dot{\eta}_i &= Q_{i+9}(a_1, a_2, \phi_1, \phi_2), \quad i = 1 \dots 9.
\end{aligned} 
\quad (30) $$

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:

$$ \{(a_1, a_2, \phi_1, \phi_2) | a_1 \in [0.01, 0.75], a_2 \in [0.01, 0.4], \\ \phi_1 \in [0, 2\pi], \phi_2 \in [0, 2\pi]\}, $$
<figure><img src="image_11.png" /><figcaption>Fig. 7. A section of one component of the invariant manifold for the system shown in Fig. 6. The displacement constraint relationship for the third transverse deflection mode $\xi_3$, is shown at phase angle $(\phi_1, \phi_2) = (0, 0)$.</figcaption></figure>

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, $(\phi_1, \phi_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 two-dimensional 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_{\phi_1} = N_{\phi_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 \in [0.01, 0.75], a_2 \in [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

$$ \{(a_1, a_2, \phi_1, \phi_2) | a_1 \in [a_1^{\text{low}}, a_1^{\text{up}}], a_2 \in [a_2^{\text{low}}, a_2^{\text{up}}], \\ \phi_1 \in [0, 2\pi], \phi_2 \in [0, 2\pi]\} \qquad (31) $$
as

$$
\begin{equation}
\begin{split}
P_i(a_1, a_2, \phi_1, \phi_2) & \approx \sum_{j=1}^{2} \sum_{k=1}^{2} \sum_{l=1}^{N_{\phi_1}} \sum_{m=1}^{N_{\phi_2}} C_{i,jklm} T_j(a_1) T_k(a_2) \\
& \quad \times F_l(\phi_1) F_m(\phi_2), \quad \text{for } i=3 \dots 18,
\end{split}
\tag{32}
\end{equation}
$$

where $T_j(a_1)$ are piecewise linear functions defined in
the amplitude segments, $a_1 \in [a_1^{\text{low}}, a_1^{\text{up}}]$, as follows:

$$
T_1(a_1) = \frac{a_1 - a_1^{\text{low}}}{a_1^{\text{up}} - a_1^{\text{low}}}, \quad T_2(a_1) = \frac{a_1^{\text{up}} - a_1}{a_1^{\text{up}} - a_1^{\text{low}}}. \quad (33)
$$

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 \in [a_2^{\text{low}}, a_2^{\text{up}}]$. The Fourier terms, $F_l(\phi_1)$ and $F_m(\phi_2)$, are defined in Eq. (14).

In each element, given by Eq. (31), the deduction
of the velocity constraint ($Q_i$) from the correspond-
ing 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 veloc-
ity constraints ($Q_i$) and the residue functions ($R_i$) are
now evaluated at the Gaussian quadrature points for
polynomials in the discretized two-dimensional ampli-
tude domain. A three-by-three-point Gaussian quadra-
ture formula is sufficient in the region, $\{(a_1, a_2) | a_1 \in
[a_1^{\text{low}}, a_1^{\text{up}}]$, $a_2 \in [a_2^{\text{low}}, a_2^{\text{up}}]\}$,
using the inner prod-
uct between the residue functions ($R_i$s) and the basis
functions defined in Eq. (32). This yields

$$
\begin{equation}
\begin{split}
& \int_{a_1^{\text{low}}}^{a_1^{\text{up}}} \int_{a_2^{\text{low}}}^{a_2^{\text{up}}} \int_0^{2\pi} \int_0^{2\pi} [T_j(a_1) T_k(a_2) F_l(\phi_1) F_m(\phi_2) \\
& \qquad \times R_i(a_1, a_2, \phi_1, \phi_2; C)] \, d\phi_1 \, d\phi_2 \, da_2 \, da_1 \\
& = 0,
\end{split}
\tag{34}
\end{equation}
\text{for } j, k = 1, 2; \ l = 1 \dots N_{\phi_1}; \ m = 1 \dots N_{\phi_2}; \\
i = 3 \dots 18.
$$

We set the number of the Fourier terms in expansion
(32) to be $N_{\phi_1} = N_{\phi_2} = 8$. As a result, the total num-
ber 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 to-
tal 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 two-
dimensional amplitude domain given by

$$
\{(a_1, a_2) | a_1 \in [0.01, 0.1157], a_2 \in [0.01, 0.06571]\},
$$

zeros are good initial values, due to the fact that the
non-linearities are weak near the origin. Then, for sub-
sequent elements, which have incremental values in
the a₁ or a₂ directions, the expansion coefficients ob-
tained from the preceding element are used as the ini-
tial values. Once the results for all discretized elements
are obtained, the expansion coefficients from contigu-
ous elements are averaged at their interface. The re-
sulting 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₃ in Eq. (30), is depicted
at the phase angles (φ₁, φ₂) = (0, 0). The amplitude
domain, a₁ ∈ [0.01, 0.75] and a₂ ∈ [0.01, 0.4], is
evenly divided into 7/7 patches, and the invariant man-
ifold appears smooth with this mesh. Note that the in-
variant manifold defined in Eq. (30) is the ensemble of
the displacement and velocity constraint relationships
for all 16 slave coordinates, and Fig. 7 represents sim-
ply the cross-section of one slave coordinate among
the 16.

Time responses for the displacements of the mas-
ter and slave coordinates are shown and compared
in Figs. 8 and 9 using three different simulation ap-
proaches: (i) direct time simulations based on the
36-state reference model, with initial conditions that
satisfy the constraint relationships; (ii) time simu-
lations 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 reduced-
order model obtained by the asymptotic expansion
method described in Pescheck et al. [18], wherein the
invariant manifold and the corresponding reduced-
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 8. Comparison of the time histories of the master coordinates $\xi_1(t)$ and $\xi_2(t)$ for the system shown in Fig. 6, with initial conditions $a_1(0) = 0.63$, $\phi_1(0) = 0.0$, $a_2(0) = 0.38$, $\phi_2(0) = 0.0$. (—) time simulation of the full reference model; (---) simulation of the reduced-order model; (----) results from the asymptotic expansion method. Plot (a) shows the first transverse modal displacement, $\xi_1(t)$; (b) second transverse modal displacement $\xi_2(t)$.</figcaption></figure>

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

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 9. Time histories of the slave coordinates $\xi_9(t)$ and $\eta_3(t)$ for the system shown in Fig. 6, for initial conditions deduced from the initial values of the master coordinates: $a_1(0) = 0.63$, $\phi_1(0) = 0.0$, $a_2(0) = 0.38$, $\phi_2(0) = 0.0$. (—) time simulation of the full reference model; (---) simulation of the reduced-order model; (----) results from the asymptotic expansion method. Plot (a) shows the ninth transverse modal displacement, $\xi_9(t)$; (b) first axial modal displacement $\eta_1(t)$.</figcaption></figure>

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.
<figure><img src="image_16.png" /><figcaption>Fig. 10. Time history of the transverse displacement at the tip of the beam for the system shown in Fig. 6, for initial conditions deduced from the initial values of the master coordinates: $a_1(0) = 0.63$, $\phi_1(0) = 0.0$, $a_2(0) = 0.38$, $\phi_2(0) = 0.0$. (—) time simulation of the reference model; (---) simulation of the reduced-order model; (---) results from the asymptotic expansion method.</figcaption></figure>

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.

## 5. 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-

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 11. Time histories of the master coordinates $a_1(t)$ and $a_2(t)$ for system shown in Fig. 6 for 1 s, under initial conditions $a_1(0)=0.63$, $\phi_1(0)=0.0$, $a_2(0)=0.38$, $\phi_2(0)=0.0$. Plot (a) shows the amplitude of the first transverse modal displacement, $a_1(t)$; (b) amplitude of the second transverse modal displacement $a_2(t)$.</figcaption></figure>

ple system and for the rotating beam system, using a combination of finite difference or finite element discretization schemes in the amplitude domain and two-dimensional 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
<figure><img src="image_19.png" /><figcaption>Fig. 12. Time histories of the master coordinates $a_1(t)$ and $a_2(t)$ for the system shown in Fig. 6 for 12 s, under initial conditions $a_1(0)=0.63$, $\phi_1(0)=0.0$, $a_2(0)=0.38$, $\phi_2(0)=0.0$. Plot (a) shows the amplitude of the first transverse modal displacement, $a_1(t)$; (b) amplitude of the second transverse modal displacement $a_2(t)$.</figcaption></figure>

extended to systems with more complicated non-linear forces, and this is relatively straightforward if the non-linear 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].

## References

[1] R.M. Rosenberg, On nonlinear vibrations of systems with many degrees of freedom, Adv. Appl. Mech. 9 (1966) 155-242.

[2] R.H. Rand, A direct method for non-linear normal modes, Int. J. Non-Linear Mech. 9 (1974) 363-368.

[3] R.H. Rand, C.H. Pak, A.F. Vakakis, Bifurcation of nonlinear normal modes in a class of two degree of freedom systems, Acta Mech. 3 (1992) 129-146.

[4] A.H. Nayfeh, S.A. Nayfeh, On nonlinear modes of continuous systems, J. Vib. Acoust. 116 (1994) 129-136.

[5] A.H. Nayfeh, S.A. Nayfeh, Nonlinear normal modes of a continuous system with quadratic nonlinearities, J. Vib. Acoust. 117 (1995) 199-205.

[6] A.H. Nayfeh, W. Lacarbonara, C. Chin, Nonlinear normal modes of buckled beams: three-to-one and one-to-one internal resonances, Nonlinear Dynam. 18 (1999) 253-273.

[7] M.E. King, A.F. Vakakis, An energy-based formulation for computing nonlinear normal modes in undamped continuous systems, J. Vib. Acoust. 116 (1994) 332-340.

[8] M.E. King, A.F. Vakakis, An energy-based approach to computing resonant nonlinear normal modes, ASME J. Appl. Mech. 63 (1996) 810-819.

[9] A.F. Vakakis, R.H. Rand, Normal modes and global dynamics of a two-degree-of-freedom non-linear system-i. low energies, Int. J. Non-Linear Mech. 27 (1992) 861-874.

[10] S.W. Shaw, C. Pierre, Normal modes for non-linear vibratory systems, J. Sound Vib. 164 (1) (1993) 85-124.

[11] S.W. Shaw, C. Pierre, Normal modes of vibration for non-linear continuous systems, J. Sound Vib. 169 (3) (1993) 319-347.

[12] A.H. Nayfeh, On direct methods for constructing nonlinear normal modes of continuous systems, J. Vib. Control 1 (1995) 389-430.

[13] A.F. Vakakis, L.I. Mikhlin, Pilipchuk, V. Yu, A.A. Zevin, Normal Modes and Localization in Nonlinear Systems, Wiley, New York, 1996.

[14] R. Camillacci, Nonlinear modal analysis for systems with geometric nonlinearities: analytical and experimental methods, Ph.D. Thesis, University of Rome Tre, Rome, Italy, 2003.

[15] E. Pescheck, C. Pierre, S.W. Shaw, A new galerkin-based approach for accurate nonlinear normal modes through invariant manifolds, J. Sound Vib. 249 (5) (2002) 971-993.

[16] D. Jiang, C. Pierre, S.W. Shaw, Large-amplitude non-linear normal modes of piecewise linear systems, J. Sound Vib. 272 (2004) 869-891.

[17] M. Legrand, D. Jiang, C. Pierre, S.W. Shaw, Nonlinear normal modes of a rotating shaft based on the invariant manifold method, Int. J. Rotating Mach. 10 (2004) 1-17.
[18] E. Pescheck, C. Pierre, S.W. Shaw, Model reduction of a nonlinear rotating beam through nonlinear normal modes, ASME J. Vib. Acoust. 124 (2) (2002) 229-236.

[19] A.H. Nayfeh, C. Chin, S.A. Nayfeh, On nonlinear normal modes of systems with internal resonance, ASME J. Vib. Acoust. 118 (1996) 340-345.

[20] M.J.D. Powell, Numerical Methods for Nonlinear Algebraic Equations, Gordon and Breach, London, 1970.

[21] P. Apiwattanalunggarn, S.W. Shaw, C. Pierre, D. Jiang, Finite-element-based nonlinear modal reduction of a rotating beam with large-amplitude motion, J. Vib. Control 9 (2003) 235-263.