Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions

Nonlinear Dynamics, Aug 2016

The analysis of whole engine rotordynamic models is an important element in the design of aerojet engines. The models include gyroscopic effects and allow for rubbing contact between rotor and stator components such as bladed discs and casing. Due to the nonlinearities inherent to the system, bifurcations in the frequency response may arise. Reliable and efficient methods to determine the bifurcation points and solution branches are required. For this purpose, a multi-harmonic balance approach is presented that allows a numerically efficient detection of bifurcation points and the calculation of both continuous and isolated branches of the frequency response functions. The method is applied to a test case derived from a commercial aeroengine. A bifurcation structure with continuous and isolated solution branches is observed and studied in this paper. The comparison with time marching based on simulations shows both accuracy and numerical efficiency of the newly developed approach.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions

Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions Loïc Salles 0 Bernard Staples 0 Norbert Hoffmann 0 Christoph Schwingshackl 0 0 B. Staples Rolls-Royce plc , PO Box 31, Derby DE24 8BJ , UK The analysis of whole engine rotordynamic models is an important element in the design of aerojet engines. The models include gyroscopic effects and allow for rubbing contact between rotor and stator components such as bladed discs and casing. Due to the nonlinearities inherent to the system, bifurcations in the frequency response may arise. Reliable and efficient methods to determine the bifurcation points and solution branches are required. For this purpose, a multi-harmonic balance approach is presented that allows a numerically efficient detection of bifurcation points and the calculation of both continuous and isolated branches of the frequency response functions. The method is applied to a test case derived from a commercial aeroengine. A bifurcation structure with continuous and isolated solution branches is observed and studied in this paper. The comparison with time marching based on simulations shows both accuracy and numerical efficiency of the newly developed approach. Rubbing; Bifurcation; Harmonic balance method; Finite element; Turbomachinery; Continuation method 1 Introduction The necessity of improving fuel efficiency of aeroengines generally requires a further decrease in clearances between rotor and stator components. This unfortunately increases the risk of contact and rubbing when operational conditions yield a disbalance of the rotor. While the design of the engine ensures well-behaved dynamic response for normal operational conditions, contact and rubbing need additional strategies to deal with the dynamic response: the interaction between rotor and stator components may generate a very complex dynamical behaviour, and as already pointed out by Sinha [ 1 ], aeroengine manufacturers have to determine the loads not only during the initial contact event but also during any continued rotation of the rotor after engine shutdown. An accurate and efficient numerical tool is necessary for this due to the extremely high cost and complexity of laboratory testing. In this paper, we present an approach to determine the vibration response for whole aeroengine rotordynamic models (WEM) in which contact and friction interactions may take place. The goal of the study is to devise a method that is numerically efficient, but at the same time allows to unfold complex bifurcation structures to a reasonable level of accuracy, including imperfect bifurcations and isolated solution branches. Many approaches are available to model nonlinear rotordynamics similar to the one under investigation here. The most common approaches are harmonic balance and transient simulation. The first one allows to quickly obtain periodic solutions, like for frequency response of the system, while the second also allows for analysis of more complex, e.g., irregular behaviour. Several studies have been devoted to the rotordynamics of systems for which contact between blades and casing occurs. Petrov [ 2 ] proposed a model with rigid blades. Aidanpaa [ 3 ] proposed a model where the blade deflection is calculated by an equivalent beam deflection model and the casing is rigid. Sinha [ 4 ] developed a model for the analysis of a system composed of blades attached to a rigid disc mounted on a shaft. His approach was also used by Lesaffre [ 5 ] and extended to a dual shaft by Gruin [ 6 ]. Parent [ 7 ] used the same model to analyse the divergence and flutter of a whole engine model with interaction in a bladed rotor-tostator contact. Choy and Padovan [ 8,9 ] studied the transient response of the bladed rotor interaction with a flexible casing. An important aspect of modelling contact with friction, which we will sometimes call rubbing in the following, is the manner of dealing with the normal contact itself. On the one hand, contact can be considered numerically as a pure unilateral contact and solved rigorously using Lagrangian multipliers [ 10 ]. On the other hand, experimental results suggest that the interaction between the rotor and an abradable lining material along the casing is more complicated and contact compliance might play a role. Ahrens et al. [ 11 ] have experimentally measured a contact normal stiffness. Kascak and Tomko [ 12 ] compared different rub models, smearing rub and abradable rub. For small penetration, both models behave in the same way. They concluded that both models had a threshold after which they caused the rotor to proceed into backward whirl. When the abradable model goes into backward whirl, it can result in the model encountering an instability, whereas the smearing model shows more benign behaviour. Following these results, the model implemented in this study is based on a compliant contact, and a contact normal stiffness is included. For simplicity, the tangential forces occurring during contact are modelled by a Coulomb friction model considering pure sliding. In the present study, we restrict ourselves to determine the steady-state periodic response of the system for periodic forcing. A harmonic balance method (HBM) is employed. Due to its numerical efficiency, this and derived methods can now be found more and more often in various application domains. Several studies have been done on harmonic balance in rotordynamics and rubbing modelling. Kim and Noah [ 13 ] developed harmonic balance with an alternating frequency/time (AFT) technique to obtain synchronous and sub-harmonic whirling motions of a horizontal Jeffcott rotor with bearing clearances. They also studied the stability of the steady-state motions using perturbation equations. The Floquet multipliers of the associated monodromy matrix were determined using a discrete HBM/AFT method. Quasi-periodic and chaotic behaviour was also observed for their Jeffcott rotor model. Peletan et al. [ 14 ] have studied the use of the harmonic balance method for rub-impact problems. They have also discussed advantages and limitations of the method compared to time marching. They concluded that HBM can predict the initiation of the quasiperiodic partial rub regime, but the partial rub and backward whirl and whip motions cannot be described by HBM. Kim and Choi [ 15 ] have developed a multiple harmonic balance method to obtain quasi-periodic response of a horizontal Jeffcott rotor with a bearing clearance. In this paper, we focus on bifurcations and multiplicity of periodic solutions. Asynchronous vibration is not considered since the analysis of the targeted whole engine model has not shown any corresponding bifurcations, such as the Neimark–Sacker type. A similar approach has, e.g., been followed by [ 16 ]. HBM has also been applied to study cyclic symmetric systems with geometric nonlinearities, and bifurcation structures have been observed and calculated [ 17,18 ]. Here we follow and expand on these ideas and apply them to a whole engine model as used in industry. It will be shown that, due to symmetry breaking, many of the bifurcations are imperfect in a weaker or stronger sense. This spoils many traditional approaches to determine bifurcation points. A strategy is thus proposed to deal with these imperfect bifurcations. The task is twofold: first, the imperfect bifurcation has to be located in parameter space, and second, detached branches, not continuously linked to the originally followed branches, have to be found. One might note that in fact the geometries of aeroengines, due to manifold design needs, usually do not show exact geometric symmetries. The problem of locating imperfect bifurcations and detached, isolated solution branches might thus be of substantial importance whenever the strength of the symmetry breaking becomes marked. The paper is set up as follows: at first, the concept of whole engine modelling, as understood in the community, is described briefly, and the underlying assumptions are stated. Then, the harmonic balance method and the contact models employed in the study are presented. Then, we show how the bifurcation analysis in the context of imperfect bifurcations and detached, isolated solution branches is conducted. An example based on a whole engine model of an aeroengine illustrates how the methods proposed may be applied. A detailed analysis of the bifurcations is provided, and the results obtained by HBM are compared to equivalent transient simulations. Finally, conclusions are given. 2 Whole engine modelling and numerical solution procedure In this study, we are most interested in the motion of the rotor centreline. The stator elements are considered in a stationary coordinate frame, and centrifugal forces and gyroscopic effects have to be included. Considering an aeroengine as an integral whole and discretising the geometry through finite elements (FE) in a straightforward manner may easily result in a model of several million of degrees of freedom when the complex geometries involved are to be captured. To yield smaller but still accurate enough models, the individual components are represented by the finite element method and component mode synthesis is used to reduce the size of each component FE model. Assembling the reduced models produces an overall whole engine model of the order of a few thousands of degrees of freedom, which can be computed by HBM. The equation of motion for the reduced whole engine model (WEM) with nonlinear contact interfaces has the following form: [M] u¨ + ([C] − [G]Ω ) u˙ + ([K] − [K]Ω ) u = p(t ) − f nl(u, u˙ ), (1) where u is the vector of generalised displacements of the reduced model. [M] , [C] and [K] are the matrices of mass, damping and stiffness, respectively. [K]Ω and [G]Ω are the matrices of centrifugal stiffening and gyroscopic effects. p(t ) is the vector of periodic excitation forces, and f nl is the vector of nonlinear forces caused by the contacts between rotor and stator elements. A periodic response of the system will be considered, with the response having the same frequency as the excitation, ω, which is given by the rotational speed of the rotor, due to an arbitrary unbalance mass p(t ) = −mdisω2 cos(ωt ). A truncated Fourier series is used to approximate the generalised displacements, u(t ) = U0 + Un,c cos(mnωt ) + Un,s sin(mnωt ), nh n=1 (2) (3) (4) where U0, Un,c and Un,s are vectors of the coefficients of the multi-harmonic representation, with mn denoting the harmonics retained in the Fourier series, and nh the truncation parameter. Inserting (2) into (1) leads to a frequency domain formulation of the problem. In frequency domain, the system can be reduced to the nonlinear degrees of freedom [ 19 ], which again permits to reduce the size of the system. The method is based on a accurate calculation of the forced response function matrix [A] (ω) = [A] + i ω [C] − ω2 [M] , [A] (ω) ≈ [A]0 + [A˜ ](ω), where [A˜ ](ω) is the modal synthesis of [A] truncated to the mode Nm and [A]0 is the static correction to provide the exact value of [A] at the frequency ω0. The residual of the nonlinear system at each frequency ω reduced to the nonlinear degrees of freedom for the j th harmonic is: R j (Qnl) = Qnl j − [A]nll (m j ω)Plj − [A]nlnl (m j ω)Pnjl + [A]nlnl (m j ω)Fnjl(Qnl) = 0 j = 1, . . . , nh. where R j is the residual function for the j th harmonic, Qnl is the vector of unknowns at the nonlinear nodes (Fourier coefficients of the displacements), Plj and Pnl j are the external forces applied to the linear and nonlinear nodes. The nonlinear forces Fnl include contact forces and gyroscopic effects. System (4) is solved using a continuation method with the frequency ω as the continuation parameter, coupled with a Newton– Raphson solver. The implementation of the contact forces follows [ 2 ]. A rotor in contact with a stator is considered. 0z is the axis of the rotor. As presented in Fig. 1, the displacements of the rotor and stator centrelines in the plane x y are (urx (t ), ury (t )) and (usx (t ), usy (t )) so that the distance between the rotor and the stator is u(t ) = ur − us + e, where e is the eccentricity of the rotor. Contact occurs when δ = u − g is positive, where g defines the clearance between the rotor and the stator. For a perfect circular shapes, the clearance is defined by g = Rs − Rr with Rs and Rr are the radius for the stator and for the rotor, respectively. The contact forces f applied to the rotor are decomposed into normal and tangential components at the contact point, f = fnn + ftt, (5) where n and t are the normal and tangential vector at the point where the contact occurs and are related to the distance u by: n = u 1 u = u ux u y and t = 1 u −u y ux . (6) The magnitude of the normal force fn is related to the normal penetration δ of the rotor into the stator. If a Hertzian contact is considered, the relation is force oriented opposite to the sliding direction, and ft = ξ μ fn, (8) (9) (10) where μ is the coefficient of friction and ξ = ±1 is a sign function, which depends on the rotor rotation and whirl directions. The expression of the contact forces in the absolute frame is: f (t ) = fx f y = fn(δ) u(t ) 1 ξ μ −ξ μ 1 ux u y . The torque Mf due to contact interactions and applied to the centreline of the structure is given by Mf = Rs ft = ξ μ Rs fn(δ). The expressions for the contact forces are not available as explicit analytic functions, so it is not possible to get explicit analytical expressions for the frequency domain representations of the contact force. An alternating frequency/time (AFT) procedure is used to get the frequency domain representation of the contact forces. The unknowns of the problem are the displacement coefficients of the centre of the rotor and the stator in the absolute frame. At each iteration of the nonlinear solution process (explained subsequently), the harmonic coefficients of the displacements are thus known and permit to get the penetration δ in time domain, using inverse discrete Fourier transform. Once the forces f (t ) and the torque Mf (t ) have been calculated, they are projected into the frequency domain via DFT (discrete Fourier transform) or FFT (fast Fourier transform). The Newton nonlinear solver needs the Jacobian matrix, which means that the derivatives [Jnl ] of the harmonic nonlinear forces F˜ nl with respect to the variables (displacements and frequency) have to be calculated, which is again achieved using the alternating frequency/time procedure [ 20 ]. The system (4) is then solved and an arc-length continuation method is employed to avoid problems with turning points. Solutions (Q, ω) are parameterised by the curvilinear parameter s defined by: fn = kcδ 23 . The tangential force is caused by friction. A Coulombtype sliding friction law is assumed, with the friction The size of the system is increased as ω becomes an unknown. A new equation has to be added. In this (7) s = QT Q + ω2. (11) study, the pseudo-arc-length method is applied and the nonlinear system to be solved is: R(Qk , ωk ) = 0 N (Qk , ωk , s) = gradk Y k − sk = 0, (12b) where R is the residual function, N is the scalar function defined by the pseudo-arc-length method, Y = (Q, ω)T is the vector of unknowns and gradk is the gradient of R at the previously calculated point Yk−1 defined by: grad = where [D]Q R is the Jacobian matrix of R and [D]ω R the gradient of R with regards to ω. A prediction–correction approach is used to build the response curve. After convergence at one frequency, a new point is predicted using the gradient and corrected solving the nonlinear system (12b). The corrector is based on a Newton–Raphson approach which is used together with a line search method, where the Newton direction is defined as the solution of a linear algebraic problem, where i is the iteration of the Newton solver. This system is not directly solved, but a block elimination technique is performed since the matrix in (14) is a bordered matrix [ 21 ]. This technique permits to factorize only the Jacobian matrix DQR, which is require by the later bifurcation analysis. In the code, a LU decomposition of DQR is performed using the LAPACK routines. One might note that the system (14) can be also solved using an iterative method based on Arnoldi iteration [ 22,23 ], which is, however, not followed here. 3 Analysis of singular points and of imperfect bifurcations Points on the frequency response function curve obtained by the above-described continuation method can be classified in regular and different types of singular points [22]: – A regular point of R(X, ω) = 0 is one for which the implicit function theorem works: RX = 0 or Rω R = 0. There is only one curve through the point. – A turning point corresponds to a singularity of the Jacobian in which zero is a simple eigenvalue of the Jacobian matrix DX R and Dω R ∈/ Im(DX R). The further properties mean that the dimension of the kernel of DY R (with Y = (X, ω)) is one. This kernel is spanned by the eigenvector Φ1 : Ker(DY R) = span [Φ1]. – A simple singular bifurcation point X0, ω0 has zero as a simple eigenvalue of DX R and Dω ∈ Im(DX R). A simple bifurcation point is characterised by a two-dimensional kernel of DY R, where one dimension corresponds to the eigenvector Φ2 = (v0, 1) and is solution of DX Rv0 + Dω R = 0, under the transversal condition Φ∗, DXX RΦ1Φ1 2 − Φ∗, DXX RΦ1Φ2 Φ∗, DXX RΦ2Φ2 > 0. Problem (12b) has exactly two different stationary solution branches intersecting at X0, ω0. – A Hopf bifurcation occurs when DX R has exactly one pair of eigenvalues ±ω0, ω0 > 0 on the imaginary axis. During the construction of the response curve (frequency response function in this case), singular points can occur and the method must be able to deal with them. In the case of a turning point, the continuation method passes the singular point without any problem, since merely a change of direction versus the frequency parameter occurs. When a simple singular point with a bifurcation is present, continuation methods without special consideration of it can experience convergence problems, since the direction of the gradient before and after the singular point may change. To overcome this issue, the singular point is usually first detected and then calculated precisely, and finally, the individual bifurcating branches arising are calculated. After this operation, each branch emanating from the singular point can be followed. Different approaches exist to detect a simple bifurcation point during the continuation procedure [ 22 ]. Here we make use of an approach proposed by Allgower [ 24 ]. It is based on the sign of the cosine between two consecutive gradients. If the cosine is negative, i.e., the two gradients have opposite direction, a simple bifurcation point exists. (15) (16) bifurcation points Fig. 2 Duffing system u¨ + 0.05u˙ + u + u3 = cos(ωt) [ 16 ] with a simple bifurcation point (a), and a modified duffing system u¨ + 0.05u˙ + u + u3 + 0.0001u2 = cos(ωt) with imperfect or perturbed bifurcations, leading to isolated solutions for the perturbed case (b). (Color figure online) To calculate the singular points precisely, a number of approaches have been developped [ 21,22,24 ]. Unfortunately, they cannot be applied directly to our problem, since the number of degrees of freedom is still too large to allow for efficient numerical implementation of them [ 22 ]. Instead, we developed the following approximate approach, which consists of four steps: 1. Approximation of the bifurcation point. Once a simple bifurcation has been detected, the bifurcation point is approximated using the secant method to find the root of the determinant H (X, ω) = DX R. 2. The vectors (Φ1, Φ2) which span the null space (kernel) of H (Xb, ωb) are approximated, by calculating the two first eigenvectors of H (Xb, ωb)∗ H (Xb, ωb). The size of this matrix is (neq + 1) ∗ (neq + 1). The span vector Ψ of the kernel of the adjoint gradient H (Xb, ωb)∗ is approximated calculating the first eigenvector of H (Xb, ωb) H (Xb, ωb)∗, whose size is neq ∗ neq 3. The algebraic bifurcation equation (ABE) defined by the following expression [ 22 ], (17) (18) Φ∗, DY Y RbΦ1Φ1 α2 c11 + 2 Φ∗, DY Y RbΦ1Φ2 αβ c12 + Φ∗, DY Y RbΦ2Φ2 β2 = 0, c22 is approximated using (Φ1, Φ2), Ψ and a numerical calculation of the coefficients: ci j = Ψ , H (Φi , Φ j ) i, j = 1, 2. A numerical finite difference method is used to approximate the second derivative ∂i ∂ j g(0, 0) = H (Φi , Φ j ) where g is defined by: g(α, β) = Φ∗ H (Yb + αΦ1 + βΦ2). (19) The following difference formulae are used: ∂12g(0, 0) = ∂22g(0, 0) = g( , 0) − 2g(0, 0) + g(− , 0) 2 , g(0, ) − 2g(0, 0) + g(0, − ) 2 , ∂1∂2g(0, 0) = g( , ) + g(− , − )4−2g( , − ) − g(− , ) , ∂2∂1g(0, 0) = ∂1∂2g(0, 0). (20) The mesh size can be chosen to be on the order of the cubic root of relative machine precision. The algebraic bifurcation equation (17) is solved and gives the coefficients αi , βi i = 1, 2. 4. Finally, the tangents τ1 and τ2 are approximated using these coefficients (τ1 = α1Φ1 + β1Φ2 and τ2 = α2Φ1 + β2Φ2). This continuation strategy works very well with proper bifurcations. Figure 2a shows as an example the application of the approach to a Duffing oscillator, as studied, e.g., in [ 16 ]. For a symmetric restoring forces, a bifurcation occurs at a frequency of about 0.69 Hz. Our approach succeeded in following the solution branches. When the symmetry of the restoring forces is broken by adding a quadratic term with a small coefficient (10−4), a perturbed, or imperfect, bifurcation results, see Fig. 2b. In that case, there are no crossing solution branches any more, and one can also not speak any more about a clear bifurcation point. Still, however, the resulting solution branches can be followed continuously, with one of the branches (branch 2) being isolated from the main branch (branch 1). From a path-following perspective, perturbing a system can be a way to deal with symmetry breaking bifurcations. Allgower calls this approach a branch switching by perturbation [ 24 ]. In the systems, we are looking at, imperfect bifurcations, indicating a lack of symmetries, do seem to be inherent to the nature of the system, however. We will see in the numerical application to a whole engine model with several rubbing elements that small even terms of the nonlinear forces often exist and create perturbated, imperfect bifurcation. With the above-described continuation technique, such imperfect bifurcations neither can be detected, nor can the isolated branches be detected or followed. This can lead to miss a whole solution branches. Moreover, in the numerical continuation scheme, there is with rubbing elements no rubbing elements 0.6 0.8 the risk of entering and being trapped in an isolated branch, which prevents building the complete frequency response function. A numerical strategy was therefore designed to deal with perturbated, imperfect bifurcations and their resulting isolated branches within continuation. When approaching an imperfect bifurcation during continuation, first the test on the cosine functions suggests the presence of a bifurcation. If the secant method to find the bifurcation point does not converge, it means that the bifurcation is perturbated. To determine the solutions branches, the approaches to calculate the continuous branch and to calculate the isolated branches Fig. 6 Frequency response of transversal force (y axis) for the second rubbing element. (Color figure online) 2 e c ro0 F − y −2 −4 −6 −80 differ. In order not to enter in the isolated branch in following the continuous branch, simply the step size s of the continuation is reduced and a flag is switched on. The flag permits to cancel the attempt to calculate a bifurcation point when the sign of the cosine between two following gradients is negative. The continuous or main branch is therefore the straightforward part in the technique: on the main branch, the step size s has plainly to be reduced to smoothly stay on the branch and overpass the imperfect bifurcation. To determine the isolated solutions, we exploit our knowledge on the existence of these other solutions. When the proximity of an imperfect bifurcation is detected as described above, the tangent at the calculated point which is in the isolated branch is saved in a file. The saved point and the tangent are used to restart the calculation and build the isolated branch. In the vicinity of the imperfect bifurcation, this approach was always sufficient to obtain convergence onto the isolated branch, which can then be continued again. A flow chart of the proposed modified continuation method for the capture of imperfect bifurcations is shown in Fig. 3. To check stability of the periodic solutions, a number of different techniques are already implemented in the available tool set and have been applied: sign of the Jacobian determinant, Floquet analysis, Hill’s criteria [ 16,25 ] and calculation of the monodromy matrix [26]. Since the techniques are all well known [ 27 ], we will not present them here. Whenever in the bifurcation 4 bifurcation 5 bifurcation 1 bifurcation 2 bifurcation 3 bifurcation 6 0.2 Frequency 0.05 0.1 0.15 following an indication on the stability of solutions is given, appropriate calculations on stability characteristics have been conducted to demonstrate either linear stability or instability. 4 Bifurcation analysis for an example system The rotor–stator contact elements and solution methods described above have been implemented in the FORTRAN 90-based computer code package FORSE developed at Imperial College London. The model subjected to the analysis is an idealisation of a three shaft aeroengine with four rubbing elements, see Fig. 4. The model consists of 1311 degrees of freedom after component mode synthesis. It was further reduced in the frequency domain to 28 degrees of freedom taking the contact into account in the rubbing elements (4 ∗ 3 ∗ 2 dofs) and gyroscopic effects (4 dofs) into account. The coefficient of friction is set to 0.1 for the four rubbing elements. The rotating system is excited by an unbalance mass at fixed eccentricity. All the resulting values (displacements, frequency, forces …) in this study have been normalised to a fixed but arbitrary value. Different values of truncation parameters for the Fourier series terms selected have been tested in the calculations, and from the results eleven harmonics have turned out to yield satisfactory results, resulting in a good compromise between accuracy and simulation time. Later on, we will show and discuss only results with full implementation of rubbing elements with corresponding normal contact force and tangential friction force models. The response functions have been calculated for a system without any rubbing elements and with rubbing elements. Corresponding frequency response functions for the centreline node of the fan are depicted in Fig. 5. Most clearly, once contacts occur in the rubbing elements, the shape of the response curve of the strongly nonlinear system strongly differs from the results without rubbing elements. Figure 6 depicts a resulting force response of the second contact element. The maximum (continuous line) and minimum (dashed line) of the force during one cycle are shown versus the frequency of excitation. For the bifurcated branches, only the maximum response values are drawn for clarity. In the considered frequency range, six bifurcations are present where five of them are imperfect. They have been numbered in ascending order with respect to their appearances in frequency. Figure 7 shows the maximum and minimum displacement of the centreline of the rotor at each rubbing element versus frequency. The initial gap values between the rotor and the stator are also shown to give an impression when contact occurs at the respective contact element. Note, however, that due to the compliance of the structure, contact does not necessarily coincide with displacements of this original gap value. Nevertheless, it can be seen that contact element 2 is the first to come into contact when the frequency increases. Figure 8 shows the contact conditions during a vibration cycle (ωt ) for the four rubbing elements over fre(a) 0.01 ten 0.008 lcem0.006 a isp 0.004 icD0.002 ronm 0 (b) a thd0H−−00..000042 e z la−0.006 i m ro−0.008 N−0.01 quency. The maps are plotted for the main branch solutions only, since the bifurcated branches do not show significantly different behaviour. The second rubbing element (Fig. 8b), which is the first to come in contact, as already been seen before. It has an intermittent contact at first, but then, the contact becomes permanent over the whole cycle for higher frequencies. The first contact element shows a very similar behaviour. The last two rubbing elements have a more complicated map with several clusters of intermittent contacts. We are now coming back to a fuller discussion of the characteristics of the six bifurcations that were detected. Figure 9 shows the bifurcation structures of all six bifurcations. As pointed out already, bifurcation number one (9a) is a simple perfect bifurcation, while bifurcations two to six are all imperfect to a lesser or stronger extent. In particular, bifurcations three to six demonstrate a remarkably complex structure. In all of the imperfect bifurcations, except for the sixth, isolated solution branches arise. Only the last bifurcation creates a continuous curve, with all the branches joining into one single curve. Figure 10 shows for bifurcations 1 and 4 how the contact forces behave. While for bifurcation 1 the solutions of the new branch differ only by a temporal phase shift of half a period and are apart from that identical, for bifurcation 4 the two new branches are actually distinct ones. The bifurcations at higher forcing frequencies are all imperfect. To determine the solutions, the continuation strategy presented above was successfully applied. Interestingly, the imperfect bifurcations are all related to a symmetry breaking. For example, there is the appearance of non-vanishing mean displacements. In frequency domain, this means that the 0th harmonic, and in consequence also higher-order even harmonics, does not vanish any more. This is shown in Fig. 11. Obviously, this effect has substantial implications for the amplitude of the resulting forces. For completeness, Fig. 12 shows two orbit plots for the centreline displacements. Two characteristic frequencies for bifurcations number five and six, 0.30 and 0.3725, have been selected. It can be seen that bifurcation number five leads to orbits with a mirror-plane symmetry, while the orbits resulting from bifurcation number six do not show this symmetry in the solution. For verification of our approach, the results of the frequency domain analysis have been compared to time 10000 9000 e8000 c ro7000 f la6000 m ro5000 n ing4000 b ub3000 R2000 1000 00 sol f1 (stable) sol f2 (unstable) sol f3 (stable) sol 1 (stable) sol 2 (unstable) sol 3 (stable) 1 2 Fig. 10 Symmetry breaking during the first bifurcation for the second rubbing element (a–c) corresponding to points (b–d) in Fig. 9a, respectively, and fourth bifurcation for the third rubbing element (d–f) corresponding to points (b–d) in Fig. 9d, respectively, a before the low-frequency bifurcation point. b Multiplicity of solutions. The bifurcated branches differ by a half-period phase shift only. c After the high-frequency bifurcation point. d Contact force versus angular position, e contact force versus angular position, d contact force versus angular position. (Color figure online) culated in the transient simulation to reach the steady state. Figure 14 shows that in the time domain simulation, where the forcing frequency is slowly, ideally adiabatically changed, one bifurcated branch is obtained during acceleration of the rotor, while the other one is reached during deceleration. Overall it turns out that the bifurcated branches calculated by HBM are excellently verified by means of the transient simulations. 5 Summary and conclusion A method has been demonstrated to calculate the frequency response functions of a whole aeroengine model based on harmonic balance method. A new strategy has been presented to detect and trace simple bifurcations and imperfect bifurcation with isolated branches. The new approach has allowed us to study the frequency response behaviour of an idealised whole engine model of a commercial aeroengine. The differMain Branch Bifurcation 2 Bifurcation 3 Bifurcation 4 Bifurcation 5 0.2 0.25 Frequency 0.3 0.35 Fig. 11 Mean or 0th harmonic displacement of the fourth contact element versus frequency. (Color figure online) domain simulation [ 28 ]. Figure 13 shows the excellent agreement between the two approaches. Some minor differences are observed in the zone of the sixth bifurcation due to the fact that not enough cycles were calFig. 12 Relative orbit of the centreline for the fourth rubbing element at a freq = 0.30, bifurcation 5, and b freq = 0.3725, bifurcation 6. (Color figure online) Fig. 13 Rubbing element 4: comparison between transient simulation (green curve) and frequency domain approach (lateral forces). (Color figure online) (a) 4 3 2 1 3 F0 O D −1 −2 −3 0 DOF 2 2 4 0 DOF 2 5 ent bifurcations occurring in the frequency response have been studied. The presence of several contact elements in the model leads to a comparatively complex bifurcation sequence in the dynamics response. The new method permits to build the full frequency response, while it is not possible with usual bifurcation analysis as the continuation algorithm enters in an isolated loop. The findings have been compared to the outcome of a time marching simulation, and excellent agreement was observed. The harmonic balance method permits to get the response in several minutes when almost day computation is necessary with time domain integration. The method was tested on a system containing 4 rubbing elements and 11 retained harmonics for the simulation which is already a large system compared to similar studies. It is possible that the method will have difficulty with very large system containing a lot of rubbing elements. The value of the determinant will increase exponentially with the number of rubbing elements. Every time a rubbing element is added, the determinant is more or less multiplied by the contact stiffness of the added rubbing element. In order to overcome this issue, the system of equation has to be rescaled, but it can be a difficult task for an industrial model. If the isolated branches are too much separated from the main branch, the continuation method would not detect them. The proposed method was not designed to find all possible solutions but to be able to build frequency response in the presence of imperfect bifurcations. Future work will comprise the implementation of a bladed rotor-to-stator contact model, which will allow to see the effects of rubbing on the dynamics of individual blades. Acknowledgments The authors are grateful to Innovate UK and Rolls-Royce Plc. for providing the financial support for this work and for giving permission to publish it. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 1. Sinha , S.K. : Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads . J. Sound Vib . 332 ( 9 ), 2253 - 2283 ( 2013 ). ISSN 0022460X. doi:10 .1016/j.jsv. 2012 . 11 .033 2. Petro , E.P. : Multiharmonic analysis of nonlinear whole engine dynamics with bladed disc-casing rubbing contacts . In: ASME Turbo Expo 2012: Turbine Technical Conference and Exposition , pp. 1181 - 1191 . American Society of Mechanical Engineers ( 2012 ). doi: 10 .1115/GT2012-68474 3. Aidanpää , J.-O.: Dynamic similarities of rotors with rubbing blades . In: 5th CHAOS 2012 International Conference , pp. 687 - 694 ( 2012 ) 4. Sinha , S.K. : Dynamic characteristics of a flexible bladedrotor with Coulomb damping due to tip-rub . J. Sound Vib . 273 ( 4-5 ), 875 - 919 ( 2004 ). ISSN 0022460X. doi:10 .1016/ S0022 -460X( 03 ) 00647 - 3 5. Lesaffre , N. , Sinou , J.-J. and Thouverez , F. : Contact analysis of a flexible bladed-rotor . Eur. J. Mech. A Solids 26 ( 3 ), 541 - 557 ( 2007 ). ISSN 09977538. doi: 10 .1016/j.euromechsol. 2006 . 11 .002 6. Gruin , M. , Thouverez , F. , Blanc , L. and Jean , P. : Nonlinear transverse steady-state periodic forced vibration of 2-dof discrete systems with cubic nonlinearities . Revue européenne de mécanique numérique 20(1-4) , 207 - 225 ( 2011 ). ISSN 17797179. doi: 10 .3166/ejcm.20. 207 - 225 7. Parent , M.- O. , Thouverez , F. and Chevillot , F. : Whole engine interaction in a bladed rotor-to-stator contact . In: ASME Turbo Expo, Number GT2014-25253. Dusseldorf ( 2014 ) 8. Choy , F.K. , Padovan , J.: Non-linear transient analysis of rotor-casing rub events . J. Sound Vib . 113 ( 3 ), 529 - 545 ( 1987 ). doi: 10 .1016/ S0022 -460X( 87 ) 80135 - 9 9. Choy , F.K. , Padovan , J.O.E. : Rub interactions of flexible casing rotor systems . Trans. ASME 111 ( 4 ), 652 - 658 ( 1989 ) 10. Legrand , M. , Batailly , A. , Magnain , B. , Cartraud , P. and Pierre , C. : Full three-dimensional investigation of structural contact interactions in turbomachines . J. Sound Vib . 331 ( 11 ), 2578 - 2601 ( 2012 ). ISSN 0022460X. doi:10 .1016/ j.jsv. 2012 . 01 .017 11. Ahrens , J. , Ulbrich , H. , Ahaus , G.: Measurement of contact forces during blade rubbing . In: Seventh International Conference on Vibration in Rotating Machinery , pp. 259 - 268 . IMechE, Nottingham ( 2000 ) 12. Kascak , A.F. : Effects of Different Rub Models on Simulated Rotor Dynamics Effects of Different Rub Models on Simulated Rotor Dynamics . Technical Report 2220 , AVSCOM 83-C-8, NASA ( 1984 ) 13. Kim , Y.B. , Noah , S.T.: Bifurcation analysis for a modified Jeffcott rotor with bearing clearances . Nonlinear Dyn . 1 ( 3 ), 221 - 241 ( 1990 ). ISSN 0924-090X. doi:10 .1007/ BF01858295 14. Peletan , L. , Baguet , S. , Jacquet-Richardet , G. , Torkhani , M. : Use and limitations of the harmonic balance method for rub-impact phenomena in rotor-stator dynamics . In: ASME Turbo Expo, Number GT2012-69450 . Copenhagen ( 2012 ). doi: 10 .1115/GT2012-69450 15. Kim , Y.B. , Noah , S.T.: Quasi-periodic response and stability analysis for a non-linear Jeffcott rotor . J. Sound Vib . 190 ( 2 ), 239 - 253 ( 1996 ) 16. Lazarus , A. , Thomas , O.: A harmonic-based method for computing the stability of periodic solutions of dynamical systems . Comptes Rendus Mécanique 338 ( 9 ), 510 - 517 ( 2010 ). ISSN 1631-0721 . doi: 10 .1016/j.crme. 2010 . 07 .020 17. Grolet , A. , Thouverez , F. : Free and forced vibration analysis of a nonlinear system with cyclic symmetry: application to a simplified model . J. Sound Vib . 331 ( 12 ), 2911 - 2928 ( 2012 ). ISSN 0022460X. doi:10 .1016/j.jsv. 2012 . 02 .008 18. Sarrouy , E. , Grolet , A. , Thouverez , F. : Global and bifurcation analysis of a structure with cyclic symmetry . Int. J. Non-linear Mech . 46 ( 5 ), 727 - 737 ( 2011 ). ISSN 00207462. doi: 10 .1016/j.ijnonlinmec. 2011 . 02 .005 19. Petrov , E.P.: A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces . J. Eng. Gas Turbines Power 133 ( 10 ) ( 2011 ). doi:10.1115/1 . 4002810 20. Salles , L. , Blanc , L. , Thouverez , F. , Gouskov , A.M. , Jean , P. : Dynamic analysis of a bladed disk with friction and frettingwear in blade attachments . In: ASME Turbo Expo 2009 , pp. 465 - 476 . American Society of Mechanical Engineers ( 2009 ) 21. Govaerts , W.J.F. : Numerical Methods for Bifurcations of Dynamical Equilibria , vol. 66 . SIAM, Philadelphia ( 2000 ) 22. Mei , Z. : Numerical Bifurcation Analysis for ReactionDiffusion Equations . Springer, Berlin ( 2000 ) 23. Salles , L. : Vibration non-linéaire des structures avec frottement:application aux modèles possédant de nombreux ddls non-linéaires . In: CSMA 2013-2011e Colloque National en Calcul des Structures , Giens, France. CSMA ( 2013 ) 24. Allgower , E.L. , Georg , K. : Numerical Continuation Methods , vol. 13 . Springer, Berlin ( 1990 ) 25. Von Groll , G. , Ewins , D.J.: The harmonic balance method with arc-length continuation in rotor/stator contact problems . J. Sound Vib . 241 ( 2 ), 223 - 233 ( 2001 ). ISSN 0022- 460X. doi:10 .1006/jsvi. 2000 .3298 26. Cardona , A. , Lerusse , A. , Géradin , M. : Fast fourier nonlinear vibration analysis . Comput. Mech . 22 ( 2 ), 128 - 142 ( 1998 ) ISSN 0178-7675 . doi: 10 .1007/s004660050347 27. Peletan , L. , Baguet , S. , Torkhani , M. , Jacquet-Richardet , G. : A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics . Nonlinear Dyn . 72 ( 3 ), 671 - 682 ( 2013 ). ISSN 0924-090X. doi:10.1007/s11071-012-0744-0 28. Williams , R.J.: Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching . In: ASME Turbo Expo , pp. 1 - 10 . Vancouver ( 2011 )

This is a preview of a remote PDF:

Loïc Salles, Bernard Staples, Norbert Hoffmann, Christoph Schwingshackl. Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions, Nonlinear Dynamics, 2016, DOI: 10.1007/s11071-016-3003-y