Non-linear normal forced vibration modes in systems with internal resonance

Abstract The new method of the forced resonance vibrations construction in mechanical systems with internal resonance is represented. According to this approach, the generalized theory of non-linear normal vibration modes by Shaw and Pierre, the modified Rauscher method and the harmonic balance method are combined with a new iterative computation procedure. The proposed approach is used in analysis of the single-disk rotor system with the isotropic-elastic shaft and the non-linear supports of Duffing type. Gyroscopic effects, asymmetrical disposition of the disk on the shaft and internal resonance are also taken into account. The NNM approach allows reducing the 8-DOF problem of the rotor dynamics to the 2-DOF non-linear system for each non-linear normal mode. Both the model of massless supports and the model of supports with inertial effects are considered. It is shown that in last case all resonance regimes are separated into two different kinds. First kind corresponds to cyclic symmetric trajectories in a system's configuration space; the second kind corresponds to centrally symmetric ones. Regimes of the first kind can be evaluated by the use of the simplified mathematical model proposed in this work. Simplified model consists only of four generalized coordinates instead of the eight initial ones.


Introduction
Non-linear normal vibrations modes (NNMs) are a generalization of the normal vibrations in linear systems. In the normal mode, a finite-dimensional non-linear system behaves like a one having a single degree of freedom. Two principal concepts of nonlinear normal vibrations modes were developed by Kauderer and Rosenberg (NNMs in conservative systems) and by Shaw and Pierre (NNMs in dissipative systems).
The Kauderer-Rosenberg concept of NNMs [1][2][3] is based on determination of trajectories in the non-linear system's configuration space. If some generalized coordinate is chosen as the independent variable, for example, x 1 ≡x, then the following functions must be obtained: These relations define the NNMs by Kauderer-Rosenberg. The relations (1) describe trajectories of some synchronous periodic motions during which all positional coordinates of the system vibrate equiperiodically. For conservative system trajectories in the system configuration space can be obtained directly from the equations of motion using the integral of energy [1][2][3]. In general, the NNM trajectories are curvilinear, in contrast with the rectilinear trajectories of linear systems.
Shaw and Pierre reformulated the concept of non-linear normal modes for a general class of non-linear dissipative systems [4,5]. Their analysis is based on the computation of invariant manifolds of motion on which the NNMs take place. In this case NNMs can be represented by dependencies between phase coordinates of the system.
Numerous publications describe problems of the NNMs existence, methods of the regimes construction, an analysis of the NNMs stability, a generalization of the NNMs concept to systems with non-smooth characteristics, etc. The NNMs theory is applied to analyze dynamics of finite-DOF mechanical systems, such as floating offshore platforms, rotors, piece-wise linear systems. The NNMs are used to describe dynamics of beams, plates and shallow shells. Target energy transfer and localization of structures motions in light of NNMs theory are treated too. Detailed reviews on theoretical and applied aspects of the NNM s approach can be found in [3,11].
In the present paper Shaw and Pierre concept of non-linear normal modes is widely used, so some detailed information about the Shaw-Pierre NNMs is represented here. Namely, one considers the finite-degree-of-freedom mechanical system, where the damping is taken into account. This system can be presented in a phase-space of the form: where s ¼ ½s 1 ; s 2 ; …; s n T . One chooses a couple of phase variables (u, v) of the non-linear dynamical system, so-called "master coordinates" (or active coordinates), where u is some dominant generalized coordinate, and v is the corresponding generalized velocity. By the Shaw-Pierre approach, the non-linear normal mode is such regime when all generalized coordinates and velocities are univalent functions of the selected couple of master variables. The master coordinates can be chosen as new independent ones instead of time. Denoting, for example, these master coordinates as the coordinate and velocity with the index 1 (that is, q 1 ¼ u; s 1 ¼ v), one writes the non-linear normal mode as where Q 1 ðu; vÞ ¼ u ; S 1 ðu; vÞ ¼ v. The selected pair of phase coordinates is used as new independent coordinates instead of time, and the functions Q i ðu; vÞ ; S i ðu; vÞ can be determined from the following system of non-linear partial differential equations: Solutions of Eq. (4) are obtained, in particular, in the form of the power series by u and v.
If we have both external and internal resonances, the method of NNMs is supplemented by the modified Rauscher method. The Rauscher method is first proposed for the single-DOF system in [12]. This method is first used to construct forced NNMs for some special system in paper [6]. General procedures of the modified Rauscher method utilization to construct NNMs in N-DOF nonautonomous systems are described in [7,2,13].
In a case of internal resonance it can observe an interaction of two NNMs. So, four phase coordinates are active, and they must be chosen as new independent variables. A general description of the Shaw-Pierre NNMs use in a case of internal resonance is presented in [14]. In this resonance case all other phase coordinates are determined as univalent functions of the four chosen master coordinates. Namely such situations may occur in the problems of the forced vibrations in rotor dynamics.
Rotor systems are very important elements of machines and mechanisms. It is well known that rotor systems show complicated behavior due to different non-linear effects. Moreover, internal resonances in the rotor systems dynamics must be taken into account. Some important publications on the rotor non-linear dynamics are selected below. Asymptotic method in non-linear dynamics of rotating shaft is first suggested in [15]. Bolotin [16] took into account a non-linear inertia in a model of the one disk rotor. Different models of rotor vibrations and analysis of motion stability are treated in the book by Tondl [17]. Different non-linear effects in the rotor dynamics are discussed in [18]. The mode locking of the whirling motions of rotor is considered in [19]. In [20] the restoring forces are expressed so as to obtain a Duffing-type equation of motion for an axi-symmetrical, two-DOF rotor. The stability of synchronous motions is analytically investigated while a numerical simulation makes possible to determine a more complete non-linear behavior. An analysis of the dynamic behavior of a rigid rotor with non-linear elastic restoring forces is carried out in [21]. It is shown that, in addition to synchronous solutions, relatively small values of the damping forces made possible the onset of the rotor precession motions which are periodic or quasiperiodic. Chaotic rotor motions are obtained for high level of unbalance and relatively high damping in the system. An experimental confirmation of the theoretical data is sought. Comparison of dynamics of linear and non-linear rotor models is presented in [22]. The periodic and chaotic vibrations under condition of internal resonance in the 2-DOF model of the Jeffcott rotor are investigated in [23] by using an asymptotic approach. Different problems of the rotor dynamics are considered in [24]. Non-linear oscillations of a rotor-magnetic bearing system under superharmonic resonance conditions are studied in [25]. Numerical simulation is used to analyze the symmetrical single-disc flexible rotor-bearing system in [26]. The 4-DOF non-linear model of the rotor dynamics is considered in [27] by using the multiple scales method.
Note that in many publications mostly the simplest models, for example, the Jeffcott rotor, are considered due to a complexity of the rotor system dynamics. The analysis of more complex models becomes a non-trivial task because of a large number of degrees of freedom and, therefore, much more complex non-linear equations of motion. It seems that one of the most appropriate approaches to analysis of the rotor steady-state dynamics is the NNMs approach.
The non-linear normal modes (NNMs) are constructed here for the rotor systems with the internal resonance. This situation is always realized in the rotor dynamics with the isotropic-elastic shaft and the isotropic-elastic supports. Two such models are considered: single disk rotor models with and without inertial forces in supports taken into account. Gyroscopic effects, nonlinear behavior of supports and an asymmetrical disposition of the disk in the shaft are taken into account too. The Shaw-Pierre NNM approach and the modified Rauscher method allow reducing the 8-DOF problem to the 2-DOF non-linear system for each nonlinear normal mode of forced vibrations.
The paper is organized as follow. A generalization of the Shaw-Pierre concept of NNMs for a case of internal resonance is presented in Section 2. An application of this concept to forced vibrations analysis is considered. The iteration procedure which combines the generalized NNMs approach and the modified Rauscher method is described. In Section 3 the basic 8-DOF model of the rotor non-linear dynamics is presented. Determination of forced vibration modes for a case of massless supports is considered in Section 4. In Section 5 inertial effects in supports are taken into account in determination of the forced vibration modes. It is shown in Section 6 that in region where resonance regimes of the rotor precession with cyclic symmetric trajectories in the system configuration space are unstable, the other solutions having center symmetric trajectories exist. For a case of the cyclic symmetry of trajectories it is possible to utilize the simplified mathematical model with four generalized coordinates instead of the eight initial ones (Section 7). A procedure of partial reduction to principal coordinates, which permits to simplify analytical and numerical procedures, is described in Section 8.

Iteration procedure to construct forced non-linear normal vibration modes in a case of internal resonance
An analysis of the resonance solutions in multi-degree-offreedom systems is made in numerous publications due to exclusive importance of the problem for theory and applications. In particular, the multiple-scale method allows us to analytically solve the equations of motion and recognize resonances [28][29][30]. It is known that all analytical transformations are essentially complicated with an increase of the number of the degrees of freedom of the system under consideration. Authors of the paper suppose that the method of NNMs in combination with some other methods has some benefits in comparison to the existing perturbation methods, in particular, in multi-DOF systems having external and internal resonances. Really, the NNMs approach allows to reduce the non-linear dynamical system to a single-DOF system when the internal resonance is absent, or to two-DOF one if the internal resonance is present. Besides, the NNMs approach is very convenient for realization of automated analytical transformations.
The procedure described in this section is intended to obtain solutions describing forced near-resonance steady-state vibrations of mechanical systems. It is iterative and combines the NNMs, Rauscher, and harmonic balance methods.
One considers the non-linear dynamical system under an external periodical excitation, reduced to principal coordinates and written in the following canonical form: Here q ¼ fq 1 ; q 2 ; …; q N g T and s ¼ fs 1 ; s 2 ; …; s N g T are general coordinates and corresponding velocities. It is assumed that the functions f i ðq; sÞ are analytical with respect to their arguments. It is assumed too that two eigenfrequencies ν 1 and ν 2 are close to the excitation frequency, Ω, that is Ω≈ν 1 ≈ν 2 . In this case two coordinates, q 1;2 , and two corresponding velocities, s 1;2 , may be treated as master (active) ones and may be taken as independent ones to construct forced NNMs. It is assumed that there is some initial approximate representation of the master coordinates in the form of the Fourier series when the other coordinates are essentially smaller than the master ones. Note that this approximate representation can be obtained, for example, by the harmonic balance method from the 2-DOF autonomous system corresponding to four phase master coordinates, or by using other possibilities. One has Presenting terms on right-hand sides of the equalities (6) as powers of the cos ðΩtÞ and sin ðΩtÞ, then using some trigonometric and algebraic manipulations, one has the following: Note that a similar relation can be obtained for the function sin ðΩtÞ. By using the relation (7), the next N-DOF "pseudo-autonomous" system is obtained instead of the system (5): Transformations above correspond to the principal idea of the Rauscher method. In the obtained autonomous system the NNMs can be constructed as functions of four new independent variables, q 1 ; s 1 ; q 2 ; s 2 , of the form: The relations (9) determine the Shaw-Pierre NNMs of the autonomous system (8) in a case of the internal resonance. One writes the system (8) in the matrix form: Passing on to the new independent variables, one presents the differentiation in time t of the form of the linear differential operator in partial derivatives L: Then the system (8) can be rewritten as One has, in the scalar form, Eliminating the time t by using the first four equations of the system (13), the partial differential equations are obtained from the other N-4 equations. Then the relations (9) can be determined as solutions of this system of PDEs, for example, in power series: Substituting the series (14) to the PDEs system and equating terms of the same powers on q 1 ; s 1 ; q 2 ; s 2 , one obtains algebraic equations with respect to unknown coefficients of the series (14), that is, a n ð Þ j ; b n ð Þ j . Note that the coefficients a n ð Þ 1 ; a n ð Þ 2 ; b n ð Þ 1 ; b n ð Þ 2 satisfy a system of non-linear algebraic equations. In general, an analysis of such systems is not a simple problem due to existence of many solutions. Some arguments permit to simplify this analysis. In fact, the active variables q 1 ; s 1 ; q 2 ; s 2 in corresponding regime of NNM have the largest amplitudes, hence the relations (9) determine sloping surfaces in the system phase space. So, values of the coefficients a n ð Þ 1 ; a n ð Þ 2 ; a n ð Þ 3 ; a n ð Þ are small, and the coefficients can be easily determined near zero by effective numerical procedure. All other coefficients satisfy a system of recurrent linear algebraic equations: using obtained values a n ð Þ 1 ; a n ð Þ 2 ; a n ð Þ 3 ; a n ð Þ After determination of the relations (9) the N-DOF system (8) is reduced to the two-DOF one for each resonance non-linear normal mode. Four master phase coordinates can be obtained from this reduced system in a form of the more precise Fourier series instead of the representation (6). As a result, the iterative process can be constructed, and the pointed out series of operations can be repeated some times to reach a necessary exactness. Numerical simulation, realized in different N-DOF non-linear systems, confirms a good exactness of the proposed approach. Note that an approach presented in this section is described in [17] for a more simple case when the internal resonance is absent.
As an example, forced vibrations of a system of four oscillators, connected by elastic springs, one of them is non-linear, is considered (Fig. 1).
Equations of motion are the following: Calculations are realized for the next values of the system parameters (dimension of the coefficients is not presented here for simplification, because this example is illustrative): For these values of parameters the system (15) can be written in principal coordinates as So, it is realized only the simple external resonance if the external excitation frequency is close to the first, or second natural frequencies (resonances 1,2); but if the external excitation frequency is close to the third (resonance 3), or fourth natural frequencies (resonance 4) one has a realization of two resonance conditions simultaneously; namely, both the external resonance and the internal one occur.
The proposed approach of the NNM determination is used to obtain the system frequency responses. The frequency responses are presented for the first harmonics of original coordinates x 1 ; x 2 ( Fig. 2) and x 3 ; x 4 ( Fig. 3).
In the case of the simple resonance (resonances 1,2 -

Principal model of the rotor dynamics
In the this section a model of the rotor dynamics with an asymmetrical disposition of the disk on the shaft is considered (Fig. 4). Gyroscopic effects, non-linearity and inertial forces in supports are taken into account. The fixed and moving coordinate systems and positional angles of the disk are shown in Fig. 5.
Equations of the rotor motion are the following: y are similar coefficients for the right support; β is a coefficient of damping in supports; ρ 1 ; ρ 2 are coefficients of damping during the disk motion; m is the disk mass; ε is an eccentricity of the disk mass center. Note that the cubic non-linearity of Duffing type can be considered as an acceptable approximation for different types of the support restoring forces [18,22,23,28,30].
One has the 8-DOF non-linear system (17), describing the displacements and rotations of the disc, and displacements of the non-linear supports. In subsequent sections two cases are considered: rotor with massless supports (Section 4) and rotor with massive supports (Sections 5 and 7). Forced resonance oscillations are obtained by use of the iterative procedure described above and compared with results obtained by use of other methods.
Consider Eq. (17) under condition m 1 ¼ m 2 ¼ 0 in the following form: In the matrix form Eq. (18) can be rewritten as One has from the second vector equation of the system (19) that Introducing the relation (20) to the first equation of the system (19), one obtains the following: Note that the linear part of Eq. (21) does not contain components of the vector x s . Eq. (21) and the second equation of the system (19) form the following modified system: Here ½K m ¼ ½c mm À½c ms ½c ss À1 ½c sm ,F m ¼ F m ðx m ; x s ; _ x m ; _ x s ; tÞÀ ½c ms ½c ss À1 F s ðx m ; x s ; _ x m ; _ x s ; tÞ. If ½U is a matrix formed by the eigenvectors of the matrix ½K m , then, making the next change of variables, x m ¼ ½Uq m , and multiplying from the left the first Eq. (22) to ½U À1 , one obtains the next system written in principal coordinates: Taking into account that ½U À1 ½U ¼ ½E and ½U À1 ½K m ½U ¼ Diagðν 1 2 ; ν 2 2 ; …; ν N 2 Þ ¼ ½ν, and making the formal change x j ¼ q j ðj ¼ k þ 1; NÞ, one obtains the following final variant of the system in principal coordinates: 1 ; q 2 ; …; _ q 1 ; _ q 2 ; …; tÞ ¼ 0 r kþ1;1 q 1 þ r kþ1;2 q 2 þ r kþ1;3 q 3 þ ⋯ þ φ kþ1 ðq 1 ; q 2 ; …; _ q 1 ; _ q 2 ; …; tÞ ¼ 0 ⋮ r N;1 q 1 þ r N;2 q 2 þ r N;3 q 3 þ ⋯ þ φ N ðq 1 ; q 2 ; …; _ q 1 ; _ q 2 ; …; tÞ ¼ 0 Then the procedure described earlier for the system (8) is used under the resonance condition Ω≈ν 1 ≈ν 2 . Two positional coordinates, q 1;2 , and two corresponding velocities, s 1;2 , are taken as master (active) coordinates to construct forced NNMs. The PDEs similar to Eq. (13) can be derived, and their solutions of the form of the powers series (14) can be obtained. When coefficients of the series (14) are known, the series must be introduced to the system (24). It permits to consider only two first equations of the system. From this reduced system new, more exact, presentations of the form (6) for four active phase coordinates can be found. So, the iteration procedure is constructed to obtain resonance forced vibrations with required exactness. Results presented in Figs. 6-8 are obtained for the next parameters of the rotor system: m ¼ 10 kg; I p ¼ 0:0405 kg m 2 ; I e ¼ 0:02858 kg m 2 ; 11 Pa, radius of the cross-section of the shaft R ¼0.018 m. A ratio of the external excitation frequency to the first fundamental frequency is equal to ω ¼ Ω=v 1 ¼ 1:0021.
The relation of initial and principal coordinates is the following: where the following designation is used: The four last lines correspond to displacements of supports. In Fig. 6 projections of the forced NNM trajectory on the planes of generalized coordinates q 1 ; q 2 (Fig. 6a); q 1 ; q 3 (Fig. 6b); q 1 ; q 5 ( Fig. 6c) are shown, where q 1 ; q 2 are active (master) coordinates. In Fig. 7 the generalized coordinates q 1 ; q 3 versus the dimensionless time τ (τ ¼ ν 1 t) are shown. Results obtained by the NNMs method, the harmonic balance method, and the verifying numerical simulation, are compared. In Fig. 8 the frequency response of the rotor system for the generalized coordinate q 1 is presented. Results obtained by the NNMs and harmonic balance method, are compared. The dimensionless frequency ω ¼Ω/ν 1 is postponed on the horisontal axe and amplitudes of the first and third harmonics of the normalized solution for the coordinate q 1 are postponed on the vertical axe.

Construction of the rotor forced vibration modes taking into account inertial forces in supports
One considers the complete one-disk rotor dynamics presented in Fig. 4 taking into account inertial forces in supports. The full 8-DOF model (17) is investigated in this case. As in the precedent sections, the generalized coordinates q 1 ; s 1 ; q 2 ; s 2 are chosen as Fig. 6. Projections of the forced NNMs trajectory on the planes q 1 ; q 2 (a); q 1 ; q 3 (b); q 1 ; q 5 (c) for the case of the inertialess supports. Lines correspond to checking numerical simulation, points correspond to results, obtained by the harmonic balance method, and cicles correspond to results obtained by the NNMs method. Fig. 7. The coordinates q 1 (a) and q 3 (b) in regime of the resonance NNMs versus dimensionless time τ (in a system with inertialess supports). Lines correspond to cheching numerical calculations, points correspond to resuls obtained by the harmonic balance method, and circles correspond to results obtained by the NNMs method. master coordinates in regime of the resonance NNM. The proposed previously procedure which joints the NNMs approach and the modified Rauscher method is used.
Results of numerical calculation presented below are obtained for the next values of the system parameters: m ¼ 18 kg; m 1 ¼ m 2 ¼ 1:8 kg; I p ¼ 0:36 kg m 2 ; I e ¼ 0:195 kg m 2 ; l ¼ 1 m; l 1 ¼ 0:3 m; 11 Pa, and a radius of the cross-section of the shaft is equal to 0.015 m. The first fundamental dimensionless frequency is equal here to 144.27.
Frequency responses of principal coordinates near the first resonance are presented in Fig. 9: Figures (a)  The relation between initial and principal coordinates is the following: where fx; y; θ 1 ; θ 2 ; x 1 ; y 1 ; x 2 ; y 2 g T ¼ fx 1 ; x 2 ; x 3 ; x 4 ; x 5 ; x 6 ; x 7 ; x 8 g T ¼ x.
Trajectories of the resonance vibrations in the system configuration space are presented in Fig. 10. Comparison of results obtained analytically and by numerical simulation shows a good efficiency of the proposed NNM approach. The corresponding time response is shown in Fig. 11, where results obtained by the NNM approach are compared with results obtained by numerical simulation. A ratio of the external excitation frequency to the first fundamental frequency is equal to ω ¼ Ω=v 1 ¼ 1:0298.

Stability and bifurcations of the forced vibration modes in 8-DOF rotor system
The frequency responses shown in Fig. 12 are obtained by the NNMs approach near the first resonance for the following parameters of the system: I e ¼ 0:1225 kg m 2 ; l ¼ 0:8 m; l 1 ¼ 0:24 m; Amplitudes of the first harmonics of the disk center displacements are shown.
Analysis of the forced NNMs stability, which is made by calculation of multiplicators, gives the unexpected results which are shown in Fig. 12. It follows from Fig. 12 that in some frequency range obtained regimes are unstable.
Additional analysis shows that in this frequency range a pair of new solutions bifurcates. It can be found by the NNMs approach, as well by the harmonic balance method. More precise frequency response for the first harmonic of the disk displacement x is shown in Fig. 13a, where solutions obtained by the harmonic balance method are presented by lines and ones obtained by the NNMs method are presented by points. Analysis of the new solutions stability gives results presented in Fig. 13b.  One considers the obtained new regimes in details. One chooses some value of the frequencies ratio, namely, Ω=ν 1 ¼ 1:02. Regimes, denoted as A, B and С (Fig. 13b), correspond to the chosen frequency value. Space representation of the rotor precession, corresponding to the regime A, is shown in Fig. 14. Besides, trajectories, which describe a motion of the disk center (point 1) and motion of the left and right supports (points 2 and 3 respectively), are too shown in Fig. 14.
A form of trajectories presented in Fig. 14 is typical for all regimes which are shown in Fig. 12, namely, trajectories of motion as the disk center as well supports are cyclic symmetric lines. But for the rotation frequency chosen earlier the regime A is unstable, and it is confirmed by numerical simulation.
Trajectories corresponding to regimes of the rotor precession B and C are not cyclic symmetric; they are center symmetric. Space representation of these trajectories is presented in Figs. 15 and 16, respectively, where it is shown trajectories describing a motion of the disk center, left and right supports (points 1, 2, 3 respectively). In these regimes the disk center circumscribes a trajectory close to the ellipse. In the regime В the vibration amplitudes in direction of the axis ОХ are essentially more than amplitudes of motion on direction OY; for the regime С everything is inverse. Each of these regimes is stable; results of checking numerical calculations are shown in Figs. 15 and 16.

Use of the model having half dimension
Preceding results show that there are two kinds of the system forced resonance solutions, namely: first one corresponds to cyclic symmetrical motions of the disk and supports (regime А), the second one corresponds to motions when trajectories of the disk and supports are central symmetric (regimes B,C) (note that the supports of the rotor considered above have identical stiffness in X and Y directions ). The cyclic symmetry of the trajectories in regime A in the system configuration space permits to find such regimes by using the simplified mathematical model.
Introducing the following definition, f ¼ fx; θ 1 ; x 1 ; x 2 g T ; f n ¼ fy; θ 2 ; y 1 ; y 2 g T , it is possible to verify that the rotor system (17) is invariant with respect to the transformations fðtÞ-f n ðt þ T=4Þ; f n ðtÞ-Àfðt þ T=4Þ under the simultaneous change t-t þ T=4 (here T ¼ 2π=Ω is the period of rotation). Such symmetry is also observed in analytical solution and corresponding trajectories in Fig. 14. So, the symmetry conditions permit to use the simplified model for regimes of the type A; a number of DOF is reduced by half. The corresponding reduced equations of motion can be written in the following form: The harmonic balance method can be directly used to solve Eq. (25). The utilization of the NNMs approach is possible under some conditions. At first, the proposed variant of the NNMs approach needs in transformation to principal coordinates. But in the regime A principal coordinates may also be separated to two groups similar to introduced earlier vectors f and f n . This allows to obtain four equations in principal coordinates corresponding to Eq. (25). Second, among these four coordinates only one principal coordinate will be active (for example, q 1 ). So, the Shaw-Pierre NNMs must be determined of the form q i ¼ q i ðq 1 ; s 1 Þ; s i ¼ s i ðq 1 ; s 1 Þ, where s 1 ¼ _ q 1 ; it essentially facilitates all calculations. At third, it is possible to assume that in zero approximation a trajectory of the motion in the plane ðq 1 ; s 1 Þ is elliptic, that is, one has q 1 ¼ A cos ðΩt þ ψÞ; s 1 ¼ ÀΩA sin ðΩt þ ψ Þ. Such representation results in the following relations: q 1 ðt þ ðπ=2ΩÞÞ ¼ ÀA sin ðΩt þ ψÞ ¼ s 1 ðtÞ=Ω; s 1 ðt þ ðπ=2ΩÞÞ ¼ ÀAΩ cos ðΩt þ ψÞ ¼ ÀΩq 1 ðtÞ. The last relations allow make a change of active variables in the reduced system, then the system of PDE of the form (9) can be constructed.
The relations q i ðt þ ðπ=2ΩÞÞ ¼ q i ðs 1 =Ω; ÀΩq 1 Þ; s i ðt þ ðπ=2ΩÞÞ ¼ s i ðs 1 =Ω; ÀΩq 1 Þ must be substituted into the series corresponding to NNM being calculated. Other aspects of the construction of the NNMs are the same as earlier. Checking numerical simulation shows that frequency responses near the resonance Ω≈Ω 1 cr: (Ω 1 cr: is the first critical speed of the rotor), constructed by the harmonic balance method applied both to complete initial system (17) and to the reduced system (25), are not different.

Simplification of the NNMs determination
The presented above procedure of the NNMs determination leads to complicated analytical transformations to construct the system (13) and to obtain coefficients a n ð Þ j ; b n ð Þ j of the series (14) in a case of internal resonance. As a result, a computer realization of the procedure needs in a long-time calculation and consumes a large amount of the computer RAM resources. It is possible to propose a modified version of the procedure, stated above. By implementing this approach one can get analytical expressions of the NNM method that are essentially simpler than original ones. Note that the system is the standard form (5), as a rule, is not the initial mechanical one. It can be obtained from the system,  which is presented in the following matrix form: where ½M and ½C are inertial and stiffness matrixes respectively, the vector Φ contains non-linear, dissipative and gyroscopic terms. Eq. (27) is obtained from (26) by multiplication with (from the left) M ½ À1 . The initial system (27) is associated with the system in principal coordinates (5) by linear change of variables, x ¼ ½Uq; _ x ¼ ½Us, where U ½ is a matrix of the eigenvectors of the linearized system. To transfer to the system in principal coordinates it is necessary to multiply (from the left) the system (27) to U ½ À1 and to make corresponding change of variables. As a result, one has the following system, which is the matrix form of the system (5): Then a modification of the procedure of NNMs construction will be presented.
The initial equation (27) is essentially simpler than Eq. (28), which are used to construct NNMs, but among variables fx 1 ; x 2 ; …; x N g it is not possible to select dominant (active) ones. On the other hand all components of the vector q can be expressed by four active phase variables (for example, q 1 ; s 1 ; q 2 ; s 2 ) as power series. But the vector x linearly depends on the vector q; so, components of the vector x can be presented as power series by the same active phase variables of the system (28). Thus NNMs can be described both by dependences (9) and by dependences x i ¼ x i ðq 1 ; s 1 ; q 2 ; s 2 Þ; y i ¼ y i ðq 1 ; s 1 ; q 2 ; s 2 Þ; ð29Þ where y i are generalized velocities of the initial system.  The main idea of the modification of the NNMs determination is in utilization of the standard procedure to use expansions x i ðq 1 ; s 1 ; q 2 ; s 2 Þ; y i ðq 1 ; s 1 ; q 2 ; s 2 Þ for phase variables of the more simple initial system (27).
Thus, together with variables of this initial system, x ¼ fx 1 ; x 2 ; …; x N g T , and generalized variables q ¼ fq 1 ; q 2 ; …; q N g T , describing the system in principal coordinates, one considers the set of variables x n ¼ fq 1 ; q 2 ; x 3 ; x 4 ; …; x N g T corresponding to some intermediate system, which is partially reduced to principal coordinates: or in matrix form: € x n þ ½K n x n þ F n ðx n ; _ x n Þ ¼ 0 Such system is fully equivalent to the systems (26) and (27), and is connected with these systems by corresponding linear transformations which are presented in details in Appendix A.
Supposing that all variables x i ði ¼ 3; NÞ in regime of the NNM can be univalently presented of the form (29), than this dependence can be determined by procedure described in Section 2, when the corresponding PDE system is solved. This system is the following: fLx n ¼ y n ; Ly n þ ½K n x n þ F n ðx n ; y n Þ ¼ 0 ð32Þ A solution is presented as power series x n n ¼ a n n ð Þ 1 q 1 þ a n n ð Þ 2 s 1 þ a n n ð Þ 3 q 2 þ a n n ð Þ 4 s 2 þ a n n ð Þ 5 q 2 1 þ ⋯; y n n ¼ b A determination of the series coefficients is described in Section 2.
It is shown in Appendix A that the pair of vectors x n and q, and the pair y n and s, are interconnected by linear relations, x n ¼ ½Pq; y n ¼ ½Ps. So, the system (32) is fully equivalent to the system described in principal coordinates.
Remark: Coefficients under linear terms of the expansions (33) must be found by numerical calculation of the system of nonlinear algebraic solutions (Section 2). In fact, a determination of these coefficients is not simple because it is not possible to select dominant variables among the phase coordinates of the system (5). So, the coefficients a n n