Abstract

Nonlinear dynamics analysis is a crucial topic in mechanical and aerospace engineering. The analysis of nonlinear normal modes (NNMs) provides an effective mathematical tool for interpreting complex nonlinear vibration phenomena. Unlike the invariant normal modes of linear systems, NNMs often exhibit frequency–energy dependence and cannot be computed using the traditional eigen-decomposition method. Calculating NNMs relies on numerical methods that involve expensive iterative computations, especially for systems with numerous degrees-of-freedom. To relieve computational costs, this article proposes a mode selection method for the component mode synthesis (CMS) technique to enable compact reduced-order modeling of structures with localized nonlinearities. The reduced-order model is then combined with a numerical continuation scheme to establish a low-cost NNM analysis framework. This framework can efficiently predict the NNMs of high-dimensional finite element models. The proposed framework is demonstrated on an I-shaped cantilever beam with localized nonlinear stiffness. The results show that the proposed approach can identify the key modes in the CMS modeling procedure. The NNMs of the nonlinear I-shaped beam structure can then be analyzed with a significant reduction in computational costs.

1 Introduction

Research on nonlinear dynamics analysis is significant for industries such as mechanical, aerospace, and civil structural engineering. Many systems encounter localized nonlinearity at mechanical joints such as connections [13] and frictional and contact joints [46]. As a result, it is critical to have effective tools for analyzing the dynamic characteristics, such as nonlinear normal modes (NNMs), of these complex systems.

The NNM concept offers a strong theoretical and mathematical framework for understanding complex nonlinear dynamic phenomena. An NNM represents the modal property of nonlinear systems in a way similar to how engineers understand traditional linear normal modes. NNMs constitute a collection of periodic solutions for undamped, unforced mechanical systems [7]. However, in contrast to the invariant normal modes of linear systems, the normal modes of nonlinear structural systems exhibit frequency–energy dependence [8]. Several approaches have been proposed for analyzing NNMs. For example, Rosenberg [9] defined the concept of NNMs as synchronized vibrations within nonlinear systems and formulated a framework for studying the nonlinear modes of multi-degree-of-freedom (DOF) systems. Peeters et al. [10] introduced a numerical continuation framework that combines numerical integration, the shooting method [11], and the pseudo-arclength continuation method [12] to compute the NNMs of nonlinear systems through iterative procedures. Although these early efforts provided effective methods for analyzing NNMs, the need for iterative computations during the analysis process often results in prohibitive computational costs, especially when analyzing high-dimensional finite element models (FEMs). Thus, reducing the computational expense of NNM analysis has become a critical research area in the nonlinear dynamics community.

In 2011, Peeters et al. proposed a reduced-order modeling framework to investigate the NNMs of a full-scale aircraft structure with localized nonlinearities [13]. This framework combines the component mode synthesis (CMS) technique [14,15] with the numerical continuation scheme [10] to enable nonlinear modal analysis of high-dimensional complex structures. In the framework, correct selection of normal modes in the CMS method is critical to ensure the accuracy of NNM analysis. However, identifying appropriate CMS normal modes is a challenging task if a system exhibits complex modal behavior. Therefore, all the normal modes within a frequency range of interest are usually selected to construct CMS models [1620], which can still result in high computational costs, especially when analyzing higher NNMs.

In this article, a mode identification procedure is proposed to identify effective CMS normal modes for constructing compact reduced-order models. The proposed mode identification procedure provides a systematic approach to construct a modal space that can well capture the mode shape of an NNM. Unnecessary CMS normal modes can be eliminated in the identification procedure to compact the reduced-order model and reduce the computational costs required to analyze NNMs of complex structures. The proposed computational framework is demonstrated on an I-shaped cantilever beam that is locally connected to strong nonlinear stiffness. The effects of the CMS mode selection on the beam's NNM analysis are discussed. Finally, the computational savings offered by the proposed computational framework are presented. The remainder of the article is organized as follows: First, the proposed computational framework is introduced. Next, the results of applying the framework to the nonlinear I-shaped beam are discussed. Finally, the conclusions and discussions are presented.

2 Methodology

In this section, a computational framework that combines the CMS technique [14,15], a new mode identification procedure, and a numerical continuation scheme [10] is introduced. This framework is developed for analyzing the NNMs of FEMs with strong localized nonlinearities.

2.1 The Component Mode Synthesis Method.

The first step of the proposed framework involves constructing a reduced-order model for an FEM whose equation of motion can be expressed as follows:
(1)
where M represents the mass matrix, K represents the stiffness matrix, x represents the physical coordinates, and fnl(x) represents the displacement-dependent nonlinear force. If the nonlinear force is applied to a subset of DOFs, then the matrices and vectors given in Eq. (1) can be partitioned into submatrix forms as follows:
(2)
where the subscripts i and b represent the linear portion and the nonlinear portion of the system, respectively. Note that the nonlinear force is only associated with the nonlinear DOF xb. Next, the fixed-interface CMS method [14] is employed to reduce the model by using two groups of modes: constraint modes (θc) and CMS normal modes (θn).
The matrix of constraint modes can be computed using
(3)
Note that the constraint modes in θc represent the static deformations of the system when a unit displacement is successively applied on each of the nonlinear DOF in xb while keeping all the other nonlinear DOFs fixed. The CMS normal modes are then defined as the linear normal modes of the system for fixing all the nonlinear DOFs. The CMS normal modes θn can be obtained by solving the eigenvalue problem associated with the linear portion of the model:
(4)
The CMS modal transformation matrix β can then be established as follows:
(5)
where I represents the identity matrix, θ~n=[θnT,0T]T represents the expanded CMS normal modes that include the fixed nonlinear DOFs, and θ~c=[θcT,IT]T represents the expanded constraint modes that include the unit displacements applied at each of the nonlinear DOFs. Note that only a subset of CMS normal modes is included in θ~n to reduce the dimension of β. The selection of CMS normal modes is discussed in Sec. 2.2. Then, β can be employed to construct a reduced-order model:
(6)
where Mrom=βTMβ,Krom=βTKβ, f¯nl=βTfnl, and q=[ηT;xbT]T. Note that η=[η1,η2,]T in the reduced coordinate q represents the modal coordinates that are associated with the CMS normal modes θ~n.

2.2 Mode Identification Procedure.

To construct an effective reduced-order model, all CMS normal modes within a frequency range of interest are usually selected to construct β [1620]. However, many of the CMS normal modes might contribute to the motion of an NNM, and the inclusion of these modes results in unnecessary computational costs. Thus, to construct a more compact and accurate reduced-order model, an iterative mode identification procedure is proposed herein to identify the CMS normal modes that have significant effects on the NNM analysis. The proposed mode identification procedure aims to construct a compact β and consists of multiple steps:

  • All the constraint modes θ~c are placed into a temporary transformation matrix β^(0) at the beginning of the procedure. Furthermore, all CMS normal modes of an FEM are used to construct a library of candidate mode shapes [θ~n(1),,θ~n(j),,θ~n(N)]. The proposed method then starts an iterative identification procedure with a starting iteration index m=1.

  • At the mth round of the identification procedure, all the candidate CMS normal modes (θ~n(1),,θ~n(p),,θ~n(N+1m)) are successively included into β^(m1) to construct candidate transformation matrices (β^(1)(m),,β^(p)(m),,β^(N+1m)(m)), that is, β^(p)(m)=[β^(m1),θ~n(p)].

  • A target NNM (the kth NNM) is then projected onto the column space of each of the candidate transformation matrices. Although the mode shapes of NNMs of a structure may vary as vibration energy increases, the kth NNM at a low-energy level can be well approximated by the kth linear normal mode of the structure by neglecting the nonlinear force fnl. Note that the kth linear normal mode Φ(k) can be obtained by solving for the kth eigenvector of [M,K]. Therefore, Φ(k) is projected onto each of (β^(1)(m),,β^(p)(m),,β^(N+1m)(m)) in this step to obtain projected mode shapes (Φ^(1)(m),,Φ^(p)(m),,Φ^(N+1m)(m)). The projection can be carried out using
    (7)
  • The similarity between Φ(k) and each of the projected mode shapes Φ^(p)(m) is computed using the modal assurance criterion (MAC) [21]. The MAC value between Φ(k) and Φ^(p)(m) can be calculated by
    (8)

The MAC value is the quantity between [0,1], where a higher value indicates a higher similarity between the two mode shapes.

  • (v) A transformation matrix β^(p)(m) that provides the highest MAC value between Φ(k) and Φ^(p)(m) is selected to be β^(m) in the mth-round identification procedure. The iteration index is then updated to m=m+1, and the CMS normal mode θ~n(p) identified in the mth- round identification, i.e., the CMS normal mode corresponding to the matrix β^(p)(m) is removed from the candidate library for the (m+1)th-round identification. Note that if the MAC value between a pair of Φ(k) and Φ^(p)(m) is close to 1, it implies that Φ(k) can be well approximated by the linear combination of the mode shapes in β^(p)(m) because Φ(k) lies on a vector space that is close to the column space of β^(p)(m).

Steps (ii) to (v) can be repeated until the MAC value converges to 1. In the iteration procedure, the dimension of β^(m) is larger than β^(m1) by 1; thus, the proposed identification procedure continues to expand the column space of the transformation matrix until the MAC value satisfies a specific convergence criterion (i.e., the MAC value is close to 1). Finally, β^(m) obtained in the final round of the identification procedure can be used as β to construct the reduced-order model in Eq. (6). The proposed mode identification procedure is summarized in Fig. 1.

Fig. 1
Mode identification procedure for analyzing the kth NNM
Fig. 1
Mode identification procedure for analyzing the kth NNM
Close modal

2.3 Numerical Continuation Scheme.

The numerical continuation scheme proposed by Peeters et al. [10] is combined with the reduced-order model Eq. (6) to analyze the NNMs of complex structures. The first step of the scheme is to find the initial conditions z0 and period T that satisfy the periodicity condition z(t+T,z0)=z(t,z0). Next, the shooting method [11] is employed to solve the following shooting function:
(9)

The Newton–Raphson method is iteratively applied to determine the values of z0 and T that minimize the residual of Eq. (9) when t=0. Note that z in Eq. (9) represents the state variables of the reduced-order model and can be expressed as z=[qT,q˙T]T. Once a set of solutions (z0,T) is obtained, the pseudo-arclength continuation approach is then utilized to find the corresponding NNM family that reflects the frequency–energy dependence of the dynamic mode.

The pseudo-arclength continuation method consists of two main stages: predictor step and corrector step. The predictor step initiates by determining an anticipated solution (q~0,T~):
(10)
where s represents the step size of the predictor and p=[pqT;pT]T represents the tangent vector along the solution curve in the solution space. Note that q0 in Eq. (10) represents the initial modal displacement. Also note that the initial modal velocity q˙0 is set to zero as a phase condition in the numerical continuation process [10]. The tangent vector p can be obtained by solving the following equation:
(11)
Equation (11) is an overdetermined system. As a result, p can be determined by calculating the pseudo-inverse of J. Following the predictor step, the procedure moves on to the corrector step. During this stage, the Newton–Raphson iteration is employed through a shooting technique until the convergence criterion of Eq. (9) is satisfied. The iterative procedure for the corrector step can be formulated as follows:
(12)
where the superscript l is the iteration index for the corrector step. The corrector (Δq0(l),ΔT(l)) in Eq. (12) can be obtained by solving the following equation:
(13)

Once the convergence criterion of the iteration procedure Eq. (12) is satisfied, the final solution q0 represents the initial modal displacement that lead to periodic motion with the period T, and the elastic potential energy of the system can then be evaluated using the initial state variables q0. After the correction procedure (Eq. (12)) has been completed, the computational procedure moves back to the prediction procedure (Eq. (9)) to find another solution set (q0,T) that represents the periodic solution of another energy level. Finally, an NNM branch can be computed by repeating Eqs. (9)(13) to find the energy frequency dependence of the NNM. A detailed description of the numerical continuation scheme can be found in Refs. [10,22].

3 Results

The proposed computational framework is applied to analyze the NNMs of an I-shaped cantilever beam, as shown in Fig. 2. The beam is modeled using an FEM with 1224 DOFs, and the full FEM is constructed using the commercial software ansys. The material and geometric properties of the beam are listed in Table 1. To test the proposed computational framework, a nonlinear spring with cubic spring constant α=108 is attached to a node at the beam's free end in the z-direction as shown in Fig. 2. Note that the nonlinear spring is intentionally placed far from the geometric center of the beam to complicate modal behavior. The mass and stiffness matrices are then extracted from ansys and converted into the script in the matlab environment. The model reduction procedure and NNM analysis are then carried out using matlab.

Fig. 2
FEM of the I-shaped cantilever beam
Fig. 2
FEM of the I-shaped cantilever beam
Close modal
Table 1

Geometric and material properties of the I-shaped cantilever beam

L4.5m
W0.5m
H0.8m
t0.05m
Young's modulus7×1010Pa
Density2700kg/m3
L4.5m
W0.5m
H0.8m
t0.05m
Young's modulus7×1010Pa
Density2700kg/m3

In the model reduction procedure, the constraint mode θc is first computed using Eq. (3). This system has only one constraint mode because it contains only a single nonlinear DOF (i.e., the DOF where the nonlinear spring is connected). This constraint mode is plotted in Fig. 3(a). Note that this constraint mode reflects the static deformation of the structure when unit displacement is applied at the nonlinear DOF. Next, the first nine CMS normal modes of the structure of fixing the nonlinear DOF are obtained by solving for Eq. (4), and these mode shapes are illustrated in Figs. 3(b)3(j). Note that these motions represent the linear dynamic modes when the nonlinear DOF (i.e., the z-component of the dot in Fig. 3) is fixed. Many of the mode shapes shown in Fig. 3 combine bending and torsional motions; thus, it is difficult to intuitively identify the key CMS normal modes that are essential in specific NNM analysis. To evaluate the effectiveness of the mode identification procedure proposed in this work, the first four NNMs of the structure are then analyzed using the proposed computational framework.

Fig. 3
Mode shapes used to construct the reduced-order model: (a) constraint mode and (b)–(j) first to ninth CMS normal modes
Fig. 3
Mode shapes used to construct the reduced-order model: (a) constraint mode and (b)–(j) first to ninth CMS normal modes
Close modal

3.1 Analysis of the First Nonlinear Normal Mode.

To investigate the convergence of NNM analysis by using a traditional CMS modeling strategy, CMS normal modes of the first few low frequencies are cumulatively included in β to construct five sizes of reduced-order models: 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF. Note that the inclusion of CMS normal modes in these models is based on the order of their frequencies (the traditional frequency-based strategy) rather than the proposed mode identification procedure. The 2DOF model is constructed using the constraint mode and the first CMS normal mode shown in Figs. 3(a) and 3(b). The 3-DOF model is constructed using the constraint mode and the first two CMS normal modes shown in Figs. 3(a)3(c). The 4DOF model is constructed using the constraint mode and the first three CMS normal modes shown in Figs. 3(a)3(d). Similarly, the 15DOF model is constructed using the constraint mode and the first 14 CMS normal modes, and the 20DOF model is constructed using the constraint mode and the first 19 CMS normal modes. The frequency–energy curves of the first NNM of each of the reduced-order models are then computed using the numerical continuation scheme and compared in Fig. 4(a). In Fig. 4(a), the frequency–energy curves show how the first natural frequency of each of these reduced-order models varies as the vibration energy changes. Figure 4(a) shows that the NNM curves exhibit a sudden change as the DOF of the reduced-order models increases from 3 to 4. However, the curves remain almost unchanged as the DOF increases from 2 to 3 and from 4 to 20. This implies that only a subset of CMS normal modes has a significant influence on the nonlinear modal behavior, emphasizing the importance of developing a mode identification procedure for NNM analysis. Furthermore, the sudden change in modal curves results in unpredictability for convergence study when the traditional frequency-based CMS modeling strategy is used for NNM analysis.

Fig. 4
(a) Frequency-energy curves of the first NNM. The modal curves of the 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF models are represented using blue, orange, green, red, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency-energy curves of the first NNM computed using the 2DOF model (blue), 3DOF model (orange), and 15DOF model (black). The 2DOF and 3DOF models are constructed using the proposed mode identification process (Color version online.)
Fig. 4
(a) Frequency-energy curves of the first NNM. The modal curves of the 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF models are represented using blue, orange, green, red, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency-energy curves of the first NNM computed using the 2DOF model (blue), 3DOF model (orange), and 15DOF model (black). The 2DOF and 3DOF models are constructed using the proposed mode identification process (Color version online.)
Close modal

Next, the proposed mode identification procedure described in Sec. 2.2 is applied to identify the CMS normal modes that have significant effects on the first NNM. Note that all the CMS normal modes (there are 1223 CMS normal modes for this model) are considered candidate mode shapes and included in the library. In this case study, the iterative mode identification procedure is conducted for two rounds. The CMS normal modes that provide the highest MAC values in each identification round and the corresponding MAC values are listed in Table 2. Note that the CMS normal modes that provide the second- and third-largest MAC values in each round of the identification procedure are also listed in Table 2 for reference. In the first-round identification procedure (m=1), the first CMS normal mode θ~n(1) is identified as the most important mode, and the transformation matrix β(1)=[θ~c,θ~n(1)] is formed to construct a 2-DOF reduced-order model. Note that this transformation matrix results in a MAC value of 0.9995 between the linear mode Φ(1) and its projection on β(1). In the second-round identification procedure (m=2), the third CMS normal mode θ~n(3) is identified and included in the transformation matrix β(2)=[θ~c,θ~n(1),θ~n(3)] to construct a 3DOF reduced-order model. This transformation matrix results in an MAC value of 1.000 between the linear mode Φ(1) and its projection on β(2). Then, the frequency–energy curves of the first NNM of the 2DOF and 3DOF reduced-order models are computed and compared with that of the 15DOF model, constructed using the frequency-based strategy, in Fig. 4(b). Figure 4(b) shows that the solution curve of the 3DOF model has excellent agreement with the 15DOF model, which is considered the exact solution in this case study. This observation implies that only θ~n(1) and θ~n(3) have significant effects on the first NNM, and these CMS normal modes can be identified in two rounds of identification procedure. The proposed method successfully excludes unnecessary CMS mode θ~n(2), which will be included in β if the traditional frequency-based strategy is used. The result shows that the proposed identification procedure speeds up the convergence rate of the CMS reduced-order modeling process.

Table 2

Mode identification procedure for the first NNM analysis

2DOF model constructed in the first-round identification procedure3DOF model constructed in the second-round identification procedure
CMS modesMACCMS modesMAC
θ~c,θ~n(1)_0.9995θ~c,θ~n(1),θ~n(3)_1.000
θ~c,θ~n(3)0.9268θ~c,θ~n(1),θ~n(6)0.9996
θ~c,θ~n(6)0.9244θ~c,θ~n(1),θ~n(12)0.9996
2DOF model constructed in the first-round identification procedure3DOF model constructed in the second-round identification procedure
CMS modesMACCMS modesMAC
θ~c,θ~n(1)_0.9995θ~c,θ~n(1),θ~n(3)_1.000
θ~c,θ~n(3)0.9268θ~c,θ~n(1),θ~n(6)0.9996
θ~c,θ~n(6)0.9244θ~c,θ~n(1),θ~n(12)0.9996

Note: The most important CMS normal modes identified in each round of the procedure are presented in bold and underline. Only the modes that provide the top three MAC values are listed in the table.

To understand the modal behavior of the first NNM, the mode shapes at its two energy levels, indicated by () and () in Fig. 4(b), are plotted in Figs. 5(a) and 5(b), respectively. These plots show the maximum deflections of the first NNM at its low- and high-energy levels. It is observed that the first NNM behaves like the first linear bending mode at the low-energy level, while it exhibits more torsional deflection as the energy increases.

Fig. 5
(a) Mode shape of the first NNM at its low-energy level () and (b) mode shape of the first NNM at its high-energy level ()
Fig. 5
(a) Mode shape of the first NNM at its low-energy level () and (b) mode shape of the first NNM at its high-energy level ()
Close modal

3.2 Analysis of the Second Nonlinear Normal Mode.

In this section, the second NNM of the structure is analyzed using the same procedure. The traditional frequency-based strategy is used to construct five sizes of reduced-order models: 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF. The frequency–energy curves of the second NNM of each of the reduced-order models are then computed and shown in Fig. 6(a). Similar to the convergence study of the first NNM, Fig. 6(a) shows that the NNM curves exhibit a sudden change as the DOF of the reduced-order models increases from 3 to 4. However, the curves remain almost unchanged as the DOF increases from 2 to 3 and from 4 to 20. Again, this implies that only a subset of CMS normal modes has a significant influence on the second NNM.

Fig. 6
(a) Frequency–energy curves of the second NNM. The solution curves of the 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF models are represented using blue- orange, green, red, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency–energy curves of the second NNM computed using the 2DOF model (blue), 3DOF model (orange), and 15DOF model (black). The 2DOF and 3DOF models are constructed using the proposed mode identification process (Color version online.)
Fig. 6
(a) Frequency–energy curves of the second NNM. The solution curves of the 2DOF, 3DOF, 4DOF, 15DOF, and 20DOF models are represented using blue- orange, green, red, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency–energy curves of the second NNM computed using the 2DOF model (blue), 3DOF model (orange), and 15DOF model (black). The 2DOF and 3DOF models are constructed using the proposed mode identification process (Color version online.)
Close modal

Next, the proposed mode identification procedure is applied to identify the CMS normal modes that have effects on the second NNM. Again, all the CMS normal modes are considered candidate mode shapes in the mode identification procedure. In this case study, the iterative mode identification procedure is conducted for two rounds. The CMS normal modes that provide the highest MAC values in each round of identifications and the corresponding MAC values are listed in Table 3. Note that the CMS normal modes that provide the second- and third-largest MAC values in each round of the identification procedure are also listed in Table 3 for reference. In the first-round identification procedure (m=1), the first CMS normal mode θ~n(1) is identified as the most important mode, and the transformation matrix β(1)=[θ~c,θ~n(1)] is formed to construct a 2DOF reduced-order model. Note that this transformation matrix results in a MAC value of 0.9986 between the linear mode Φ(2) and its projection on β(1). In the second-round identification procedure (m=2), the third CMS normal mode θ~n(3) is identified and included in the transformation matrix β(2)=[θ~c,θ~n(1),θ~n(3)] to construct a 3DOF reduced-order model. This transformation matrix results in a MAC value of 1.000 between the linear mode Φ(2) and its projection on β(2). Then, the frequency–energy curves of the second NNM of these reduced-order models are computed and compared with that of the 15DOF model, constructed using the frequency-based strategy, in Fig. 6(b). Figure 6(b) shows that the 3DOF model constructed using the proposed mode identification method has excellent agreement with the exact 15DOF model. This observation implies that the proposed method is able to identify θ~n(1) and θ~n(3) that have significant effects on the second NNM and exclude the unnecessary mode θ~n(2).

Table 3

Mode identification procedure for the second NNM analysis

2DOF model constructed in the first-round identification procedure3DOF model constructed in the second-round identification procedure
CMS modesMACCMS modesMAC
θ~c,θ~n(1)_0.9986θ~c,θ~n(1),θ~n(3)_1.000
θ~c,θ~n(3)0.2249θ~c,θ~n(1),θ~n(6)0.9989
θ~c,θ~n(6)0.1121θ~c,θ~n(1),θ~n(12)0.9988
2DOF model constructed in the first-round identification procedure3DOF model constructed in the second-round identification procedure
CMS modesMACCMS modesMAC
θ~c,θ~n(1)_0.9986θ~c,θ~n(1),θ~n(3)_1.000
θ~c,θ~n(3)0.2249θ~c,θ~n(1),θ~n(6)0.9989
θ~c,θ~n(6)0.1121θ~c,θ~n(1),θ~n(12)0.9988

Note: The most important CMS normal modes identified in each round of the procedure are presented in bold and underline. Only the modes that provide the top three MAC values are listed in the table.

Finally, the mode shapes of the second NNM at two energy levels, indicated by () and () in Fig. 6(b), are plotted in Figs. 7(a) and 7(b), respectively. Figure 7 shows that the geometric center of the I-beam at the free end remains stationary and behaves like the second linear torsion mode at the low-energy level, while this point exhibits translational motion in the horizontal direction as the energy increases.

Fig. 7
(a) Mode shape of the second NNM at its low-energy level () and (b) mode shape of the second NNM at its high-energy level ()
Fig. 7
(a) Mode shape of the second NNM at its low-energy level () and (b) mode shape of the second NNM at its high-energy level ()
Close modal

3.3 Analysis of the Third Nonlinear Normal Mode.

In this section, the third NNM of the structure is discussed. Again, the frequency-based strategy is employed to construct five sizes of reduced order: 3DOF, 4DOF, 5DOF, 15DOF, and 20DOF. The frequency–energy curves of the third NNM of each of these reduced-order models are then computed and plotted in Fig. 8(a). Figure 8(a) shows that the frequencies of the five reduced-order models do not change with increasing energy, and all the modal curves are almost identical. This result indicates that the frequency of the third NNM is not affected by the nonlinear spring and the NNMs exhibit good convergence starting from the 3DOF model.

Fig. 8
(a) Frequency–energy curves of the third NNM. The solution curves of the 3DOF, 4DOF, 5DOF, 15DOF, and 20DOF models are represented using red, orange, green, blue, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency–energy curves of the third NNM computed using the 2DOF model (red) and 20DOF model (black). The 2DOF models are constructed using the proposed mode identification process (Color version online.)
Fig. 8
(a) Frequency–energy curves of the third NNM. The solution curves of the 3DOF, 4DOF, 5DOF, 15DOF, and 20DOF models are represented using red, orange, green, blue, and black, respectively. The inclusion of CMS normal modes in these models is based on the order of their frequencies. (b) Frequency–energy curves of the third NNM computed using the 2DOF model (red) and 20DOF model (black). The 2DOF models are constructed using the proposed mode identification process (Color version online.)
Close modal

Next, the proposed mode identification procedure is conducted to identify the key CMS normal modes. The CMS normal mode that provides the highest MAC values in single-round identification and the corresponding MAC value is listed in Table 4. The CMS normal modes that provide the second- and third-largest MAC values in the identification procedure are also listed in Table 4 for reference. Table 4 shows that the second CMS normal mode θ~n(2), which is excluded in the first and second NNM analyses, is identified as the key CMS normal mode in the third NNM analysis. The transformation matrix β=[θ~c,θ~n(2)] constructed using this CMS normal mode results in a MAC value of 1.0000 between Φ(3) and its projection on β. Then, the frequency–energy curve of the third NNM of the 2DOF reduced-order model constructed using β is computed and compared with that of the 20DOF model, constructed using the frequency-based strategy, in Fig. 8(b). Figure 8(b) shows that the 2DOF model constructed using the proposed method agrees well with the exact 20DOF model. This observation implies that only θ~n(2) has significant effects on the third NNM, and this mode can be identified using the proposed method.

Table 4

Mode identification procedure for the third NNM analysis

2DOF model constructed in the first-round identification procedure
CMS modesMAC
θ~c,θ~n(2)_1.000
θ~c,θ~n(41)0.0324
θ~c,θ~n(5)0.0165
2DOF model constructed in the first-round identification procedure
CMS modesMAC
θ~c,θ~n(2)_1.000
θ~c,θ~n(41)0.0324
θ~c,θ~n(5)0.0165

Note: The most important CMS normal modes identified is presented in bold and underline.

Finally, the mode shape of the third NNM is plotted in Fig. 9 to investigate its modal behavior. It can be observed that the third NNM is a bending mode in the y-direction, which is unaffected by the nonlinear spring applied in the z-direction, making the third NNM a linear mode. Additionally, θ~n(2), identified by the proposed modal selection procedure, is also bending in the y-direction, as shown in Fig. 3(c), consistent with the modal behavior of the third NNM. This also explains why θ~n(2) alone can fully capture the motion of the third NNM.

Fig. 9
Mode shape of the third NNM
Fig. 9
Mode shape of the third NNM
Close modal

3.4 Analysis of the Fourth Nonlinear Normal Mode.

In this section, the fourth NNM of the structure is analyzed. Since the fourth NNM requires many CMS normal modes to capture its modal behavior, the convergence analysis on its modal curves is conducted with the proposed mode identification procedure rather than the traditional frequency-based method. The proposed mode identification procedure is conducted for several rounds to identify the key CMS normal modes. Five sizes of reduced-order models are constructed: 9DOF, 15DOF, 21DOF, 29DOF, and 37DOF. The CMS normal modes used in these models and the MAC values between the linear mode Φ(4) and its projection on each of β(m) are detailed in Table 5. The frequency–energy curves of the fourth NNM of each of the reduced-order models are then computed and shown in Fig. 10.

Fig. 10
(a) Frequency–energy curves of the fourth NNM. (b) Enlarged plot of the rectangular area. The solution curves of the 9DOF, 15DOF, 21DOF, 29DOF, and 37DOF models are represented using blue, magenta, green, yellow, and red, respectively. These models are constructed using the proposed mode identification process (Color version online.)
Fig. 10
(a) Frequency–energy curves of the fourth NNM. (b) Enlarged plot of the rectangular area. The solution curves of the 9DOF, 15DOF, 21DOF, 29DOF, and 37DOF models are represented using blue, magenta, green, yellow, and red, respectively. These models are constructed using the proposed mode identification process (Color version online.)
Close modal
Table 5

CMS normal modes identified for the fourth NNM analysis

ModelModesMAC
9DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12)_,θ~n(10),θ~n(8)0.9999
15DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19)_,θ~n(15),θ~n(14)1.0000
21DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39)_,θ~n(36)1.0000
29DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39),θ~n(36),θ~n(5),θ~n(18),θ~n(42),θ~n(30),θ~n(37),θ~n(43),θ~n(31),θ~n(50)_1.0000
37DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39),θ~n(36),θ~n(5),θ~n(18),θ~n(42),θ~n(30),θ~n(37),θ~n(43),θ~n(31),θ~n(50),θ~n(52),θ~n(47),θ~n(41),θ~n(32),θ~n(54),θ~n(53),θ~n(44),θ~n(59)_1.0000
ModelModesMAC
9DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12)_,θ~n(10),θ~n(8)0.9999
15DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19)_,θ~n(15),θ~n(14)1.0000
21DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39)_,θ~n(36)1.0000
29DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39),θ~n(36),θ~n(5),θ~n(18),θ~n(42),θ~n(30),θ~n(37),θ~n(43),θ~n(31),θ~n(50)_1.0000
37DOFθ~c,θ~n(3),θ~n(4),θ~n(1),θ~n(6),θ~n(7),θ~n(12),θ~n(10),θ~n(8),θ~n(16),θ~n(11),θ~n(17),θ~n(19),θ~n(15),θ~n(14),θ~n(22),θ~n(24),θ~n(20),θ~n(25),θ~n(39),θ~n(36),θ~n(5),θ~n(18),θ~n(42),θ~n(30),θ~n(37),θ~n(43),θ~n(31),θ~n(50),θ~n(52),θ~n(47),θ~n(41),θ~n(32),θ~n(54),θ~n(53),θ~n(44),θ~n(59)_1.0000

Note: The order of CMS normal modes listed for each reduced-order model is based on the order of identification in the iterative mode identification procedure rather than the order of mode index. The highest CMS normal modes used in each of the reduced-order models are presented in bold and underline.

Figure 10 shows that the primary frequency–energy dependency of the NNM exhibits good convergence across these five models; the MAC values of these models are all greater than 0.9999. Additionally, all five models show the “tongue” phenomenon, which represents the internal resonance behavior of the nonlinear system. However, this “tongue” phenomenon has not yet converged at 9DOF. When the model is increased to 15DOF, the original tongue shifts to a lower frequency, and an additional tongue appears at a higher frequency. When the model is increased to 21DOF, both tongues shift to lower frequencies. As the model is further increased to 29DOF and 37DOF, only one tongue remains, and it finally achieves stable convergence. The result indicates two points:

  1. An MAC value greater than 0.9999 is a good indicator for the convergence of primary frequency–energy dependency.

  2. To fully capture internal resonance, more CMS normal modes are still required.

To further investigate the internal resonance, the time history of the periodic motion at the tongue area is plotted in Fig. 11. Figure 11 shows that the internal resonance possesses a 1:3 frequency ratio.

Fig. 11
Time histories of the displacements at internal resonance of the 29DOF model: (a) nonlinear DOF xb, (b) first modal DOF η1, (c) second modal DOF η2, and (d) third modal DOF η3
Fig. 11
Time histories of the displacements at internal resonance of the 29DOF model: (a) nonlinear DOF xb, (b) first modal DOF η1, (c) second modal DOF η2, and (d) third modal DOF η3
Close modal

Next, the computational savings provided by the proposed mode identification procedure is investigated. In this case study, the 2-DOF model has reached a good convergence, as shown in Fig. 10, and the highest CMS normal mode used in this model is θ~n(50). Therefore, the modal curve of this proposed 29DOF model is compared with that of the frequency-based 51DOF model (whose highest CMS normal mode is also θ~n(50)) and the frequency-based 29DOF model shown in Fig. 12. Figure 12 shows that the modal curve of the 29DOF model constructed using the proposed mode identification method agrees well with the 51DOF model constructed using the frequency-based strategy. Additionally, the 29DOF model constructed using the proposed mode identification method shows much better convergence in the internal resonance region compared to the 29DOF model constructed using the traditional frequency-based strategy. This result shows that using the proposed method can save 22DOFs, significantly improving the efficiency of reduced-order modeling. The computational costs required for obtaining the modal curves for the proposed 29DOF model and the frequency-based 51DOF model are compared in Table 6. Table 6 shows that 61% of computational time can be saved as the model size is reduced from 51DOFs to 29DOFs since the proposed mode identification procedure is able to eliminate unnecessary CMS normal modes in the reduced-order modeling process.

Fig. 12
(a) Frequency–energy curves of the fourth NNM. (b) Enlarged plot of the rectangular area. The modal curves of the proposed 29DOF model, frequency-based 51DOF model, and frequency-based 29DOF model are represented using yellow, black, and blue, respectively (Color version online.)
Fig. 12
(a) Frequency–energy curves of the fourth NNM. (b) Enlarged plot of the rectangular area. The modal curves of the proposed 29DOF model, frequency-based 51DOF model, and frequency-based 29DOF model are represented using yellow, black, and blue, respectively (Color version online.)
Close modal
Table 6

Computational costs of the proposed 29DOF model and the frequency-based 51DOF model

Elapsed time (s) (350 solution points)
29-DOF3409
51-DOF8699
Elapsed time (s) (350 solution points)
29-DOF3409
51-DOF8699

Finally, the mode shapes of the fourth NNM at three energy levels, indicated by (), (), and () in Fig. 12, are plotted in Figs. 13(a)13(c), respectively. These plots show the maximum deflections of the fourth NNM in the low-energy, high-energy, and internal resonance regions. The fourth NNM behaves like the fourth linear bending mode at the low-energy level, while it exhibits more torsional deflection as the energy increases. In the internal resonance region, the motion becomes more involved since the dynamics combines the bending and torsional motions with different frequency components, as illustrated in Fig. 11.

Fig. 13
(a) Mode shape of the fourth NNM at its low-energy level (), (b) mode shape of the fourth NNM at its high-energy level (), and (c) mode shape of the fourth NNM at its internal resonance region ()
Fig. 13
(a) Mode shape of the fourth NNM at its low-energy level (), (b) mode shape of the fourth NNM at its high-energy level (), and (c) mode shape of the fourth NNM at its internal resonance region ()
Close modal

4 Conclusions and Discussion

In this article, an efficient computational framework is presented for analyzing the NNMs of complex structures with localized nonlinearities. The framework incorporates a mode identification procedure that combines the modal assurance criterion with a novel iterative selection process to identify the key CMS normal modes. These CMS modes are then used to construct an effective CMS model, and the NNMs of the CMS model are computed using a numerical continuation scheme. This framework provides an efficient tool for analyzing NNMs of high-dimensional localized nonlinear systems.

The effectiveness of the proposed framework is demonstrated on a finite element model of an I-shaped cantilever beam with strong localized nonlinear stiffness. Application of the framework to this nonlinear beam structure successfully captures the system's nonlinear behavior, including internal resonance phenomena, while achieving significant computational savings. The key contributions of this article include:

  1. A standardized procedure is proposed to eliminate unnecessary CMS normal modes during the reduced-order modeling process for NNM analyses, which enables a more compact model compared to the traditional frequency-based mode selection method.

  2. The proposed mode identification method iteratively applies the modal assurance criterion to identify effective CMS normal modes by comparing the similarity between an NNM's motion, at its low-energy level, and its projection onto the reduced space. The proposed iterative procedure is computationally efficient, as it utilizes linear techniques.

  3. Unpredictability in NNM convergence studies (i.e., the sudden change in modal curves in the traditional frequency-based method) can be avoided using the proposed mode identification procedure, as unnecessary CMS modes are eliminated in the proposed iterative identification procedure.

Future work will focus on:

  1. Enhancing the mode identification method to further reduce the number of CMS normal modes required to capture internal resonance motions.

  2. Expanding the mode identification method to multicomponent structures with complex nonlinear interfaces.

Footnote

1

An earlier version of this work was presented at the ASME 2024 International Design Engineering Technical Conferences.

Funding Data

  • The National Science and Technology Council (Taiwan, R.O.C) (Grant No. 112-2222-E-007-004-MY2).

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

l =

iteration index of the corrector step

s =

step size of the predictor step

p =

tangent vector of the solution curve

q =

reduced-order coordinate

x =

physical coordinate

H =

shooting function

I =

identity matrix

K =

stiffness matrix

M =

mass matrix

fnl =

nonlinear force

xi =

physical displacement of linear DOFs

xb =

physical displacement of nonlinear DOFs

Krom =

reduced-order stiffness matrix

Mrom =

reduced-order mass matrix

(q0,T) =

initial value in the predictor step

(Δq0,ΔT) =

corrector

(q~0,T~) =

predicted solutions in the predictor step

β =

CMS transformation matrix

β^ =

temporary transformation matrix

θ~c =

matrix of constraint modes

θ~n =

matrix of CMS normal modes

η =

modal coordinate

Φ =

linear normal mode

Φ^ =

projected linear normal mode

References

1.
Jaumouillé
,
V.
,
Sinou
,
J.-J.
, and
Petitjean
,
B.
,
2010
, “
An Adaptive Harmonic Balance Method for Predicting the Nonlinear Dynamic Responses of Mechanical Systems—Application to Bolted Structures
,”
J. Sound Vib.
,
329
(
19
), pp.
4048
4067
.
2.
Lacayo
,
R.
,
Pesaresi
,
L.
,
Groß
,
J.
,
Fochler
,
D.
,
Armand
,
J.
,
Salles
,
L.
,
Schwingshackl
,
C.
,
Allen
,
M.
, and
Brake
,
M.
,
2019
, “
Nonlinear Modeling of Structures With Bolted Joints: A Comparison of Two Approaches Based on a Time-Domain and Frequency-Domain Solver
,”
Mech. Syst. Signal Process.
,
114
(
1
), pp.
413
438
.
3.
Zucca
,
S.
,
Firrone
,
C. M.
, and
Gola
,
M. M.
,
2012
, “
Numerical Assessment of Friction Damping at Turbine Blade Root Joints by Simultaneous Calculation of the Static and Dynamic Contact Loads
,”
Nonlinear Dyn.
,
67
(
3
), pp.
1943
1955
.
4.
Allara
,
M.
,
2009
, “
A Model for the Characterization of Friction Contacts in Turbine Blades
,”
J. Sound Vib.
,
320
(
3
), pp.
527
544
.
5.
Krizak
,
T.
,
Kurstak
,
E.
, and
D’Souza
,
K.
,
2023
, “
Experimental and Computational Study of a Rotating Bladed Disk With Under-Platform Dampers
,”
AIAA J.
,
61
(
10
), pp.
4717
4727
.
6.
Lawal
,
I. G.
, and
Brake
,
M. R.
,
2023
, “
Energy Dissipation on an Elastic Interface as a Metric for Evaluating Three Friction Models
,”
ASME J. Appl. Mech.
,
90
(
8
), p.
081002
.
7.
Renson
,
L.
,
Kerschen
,
G.
, and
Cochelin
,
B.
,
2016
, “
Numerical Computation of Nonlinear Normal Modes in Mechanical Engineering
,”
J. Sound Vib.
,
364
(
3
), pp.
177
206
.
8.
Vakakis
,
A. F.
,
2001
,
Normal Modes and Localization in Nonlinear Systems
,
Springer
,
New York
.
9.
Rosenberg
,
R. M.
,
1962
, “
The Normal Modes of Nonlinear n-Degree-of-Freedom Systems
,”
ASME J. Appl. Mech.
,
29
(
1
), pp.
7
14
.
10.
Peeters
,
M.
,
Viguié
,
R.
,
Sérandour
,
G.
,
Kerschen
,
G.
, and
Golinval
,
J.-C.
,
2009
, “
Nonlinear Normal Modes, Part II: Toward a Practical Computation Using Numerical Continuation Techniques.
,”
Mech. Syst. Signal Process.
,
23
(
1
), pp.
195
216
.
11.
Nayfeh
,
A. H.
, and
Balachandran
,
B.
,
2008
,
Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods
,
John Wiley & Sons
,
Hoboken, NJ
.
12.
Chan
,
T. F.
,
1984
, “
Newton-Like Pseudo-Arclength Methods for Computing Simple Turning Points
,”
SIAM J. Sci. Stat. Comput.
,
5
(
1
), pp.
135
148
.
13.
Peeters
,
M.
,
Kerschen
,
G.
,
Golinval
,
J. C.
,
Stéphan
,
C.
, and
Lubrina
,
P.
,
2011
, “
Nonlinear Normal Modes of a Full-Scale Aircraft
,”
Modal Analysis Topics, Volume 3: Proceedings of the 29th IMAC, a Conference on Structural Dynamics
,
Jacksonville, FL
,
Jan. 31–Feb. 3
,
Springer
, pp.
223
242
.
14.
Craig
,
R. R.
, and
Bampton
,
M. C.
,
1968
, “
Coupling of Substructures for Dynamic Analyses
,”
AIAA J.
,
6
(
7
), pp.
1313
1319
.
15.
Hurty
,
W. C.
,
1965
, “
Dynamic Analysis of Structural Systems Using Component Modes
,”
AIAA J.
,
3
(
4
), pp.
678
685
.
16.
Kuether
,
R. J.
, and
Allen
,
M. S.
,
2013
, “
Nonlinear Modal Substructuring of Systems With Geometric Nonlinearities
,”
54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
,
Boston, MA
,
Apr. 8–11
.
17.
Kuether
,
R. J.
, and
Allen
,
M. S.
,
2014
, “
Substructuring with Nonlinear Reduced Order Models and Interface Reduction with Characteristic Constraint Modes
,”
55th AIAA/ASMe/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference
,
National Harbor, MD
,
Jan. 13–17
.
18.
Kuether
,
R. J.
, and
Allen
,
M. S.
,
2014
, “
Craig-Bampton Substructuring for Geometrically Nonlinear Subcomponents
,”
Dynamics of Coupled Structures, Volume 1: Proceedings of the 32nd IMAC, A Conference and Exposition on Structural Dynamics, 2014
,
Orlando, FL
,
Feb. 3–6
,
Springer
, pp.
167
178
.
19.
Tien
,
M.-H.
,
Hu
,
T.
, and
D'Souza
,
K.
,
2018
, “
Generalized Bilinear Amplitude Approximation and X-Xr for Modeling Cyclically Symmetric Structures With Cracks
,”
ASME J. Vib. Acoust.
,
140
(
4
), p.
041012
.
20.
Tien
,
M.-H.
,
Hu
,
T.
, and
D’Souza
,
K.
,
2019
, “
Statistical Analysis of the Nonlinear Response of Bladed Disks With Mistuning and Cracks
,”
AIAA J.
,
57
(
11
), pp.
4966
4977
.
21.
Pastor
,
M.
,
Binda
,
M.
, and
Harčarik
,
T.
,
2012
, “
Modal Assurance Criterion
,”
Procedia Eng.
,
48
, pp.
543
548
.
22.
Huang
,
S.-C.
, and
Tien
,
M.-H.
,
2024
, “
A Hybrid Continuation Framework for Analyzing Nonlinear Normal Modes of Systems With Contact Nonlinearity
,”
ASME J. Comput. Nonlinear Dyn.
,
19
(
7
), p.
071008
.