Development of a Fully-Coupled Harmonic Balance Method and a Refined Energy Method for the Computation of Flutter-Induced Limit Cycle Oscillations of Bladed Disks with Nonlinear Friction Contacts
Christian Berthold, Johann Gross, Christian Frey, Malte Krack
DDevelopment of a Fully-Coupled Harmonic Balance Method and a RefinedEnergy Method for the Computation of Flutter-Induced Limit Cycle Oscillationsof Bladed Disks with Nonlinear Friction Contacts
Christian Berthold b , Johann Gross a , Christian Frey b , Malte Krack a a Institute of Aircraft Propulsion Systems, University of Stuttgart, Pfaffenwaldring 6, 70569 Stuttgart, [email protected], [email protected] b Institute of Propulsion Technology, German Aerospace Center, Linder Höhe, 51147 Cologne, [email protected], [email protected]
Abstract
Flutter stability is a dominant design constraint of modern gas and steam turbines. To further increase the feasibledesign space, flutter-tolerant designs are currently explored, which may undergo Limit Cycle Oscillations (LCOs) ofacceptable, yet not vanishing, level. Bounded self-excited oscillations are a priori a nonlinear phenomenon, and canthus only be explained by nonlinear interactions such as dry stick-slip friction in mechanical joints. The currentlyavailable simulation methods for blade flutter account for nonlinear interactions, at most, in only one domain, thestructure or the fluid, and assume the behavior in the other domain as linear. In this work, we develop a fully-couplednonlinear frequency domain method which is capable of resolving nonlinear flow and structural effects. We demonstratethe computational performance of this method for a state-of-the-art aeroelastic model of a shrouded turbine bladerow. Besides simulating limit cycles, we predict, for the first time, the phenomenon of nonlinear instability, i.e. , asituation where the equilibrium point is locally stable, but for sufficiently strong perturbation (caused e.g. by animpact), the dry frictional dissipation cannot bound the flutter vibrations. This implies that linearized theory doesnot necessary lead to a conservative design of turbine blades. We show that this phenomenon is due to the nonlinearcontact interactions at the tip shrouds, which cause a change of the vibrational deflection shape and frequency, whichin turn leads to a loss of aeroelastic stability. Finally, we extend the well-known energy method to capture these effects,and conclude that it provides a good approximation and is useful for initializing the fully-coupled solver.Keywords: fluid-structure interaction, flutter, limit cycle oscillation, aircraft engine, friction damping, harmonicbalance
Nomenclature Δ 𝑊 Work per cyclei imaginary unit, √− (cid:61) Imaginary part 𝝌 Finite volume mesh location 𝝍 Modal deflection shape 𝒆 𝑗 Unit vector pointing in the 𝑗 -th coordinate di-rection 𝒇 Vector of forces 𝑮 Aerodynamic influence coefficient matrix 𝑲 Stiffness matrix 𝑴 Mass matrix 𝑹 Residual vector 𝒖 Vector of structural coordinates 𝒙 spatial coordinate 𝜔 Angular oscillation frequency (cid:60)
Real partˆ 𝑢 Sensor
Displacement amplitude of sensor node in cir-cumferential direction 𝐷 Modal damping ratio
Preprint submitted to JFS February 12, 2021 a r X i v : . [ c s . C E ] F e b Aerodynamic influence coefficient 𝐻 Harmonic truncation order 𝐿 Span
Blade span 𝑇 Time period 𝑡 time 𝑢 DisplacementFSI Fluid-Structure InteractionHB Harmonic BalanceLCO Limit Cycle OscillationMAC Modal assurance criterionND Nodal diameter
Superscripts ˆ Fourier coefficientH HermitianT Transposea Aerodynamic quantitycb Refers to the Craig-Bampton modelc Friction contactfem Refers to the Finite Element modelnl Refers to quantities of the refined energy methodstick Refers to fully sticking contact conditionss Solid quantity
1. Introduction
Blade vibrations cause high-cycle fatigue and therefore put at risk the safe operation of gas and steam turbines.Besides resonances with rotor-synchronous forcing, flutter is one of the most important causes for blade vibrations.Flutter denotes a dynamically unstable interaction between an elastic structure and a surrounding flow. In turbomachin-ery, flutter is mainly caused by the aerodynamic interference among the vibrating blades within a row (cascade effect)[1]. The current design paradigm is to completely avoid flutter-induced vibrations. This is achieved by linear vibrationtheory: The normal modes of vibration of the linearized structural model are determined, and for each relevant mode,the aerodynamic damping is computed assuming infinitesimally small vibrations. If the sum of the (possibly negative)aerodynamic damping and the positive structural damping (estimated from experience) of each mode is positive, thedesign is accepted; otherwise it is viewed as aeroelastically unstable and discarded. Flutter avoidance has become adominant design constraint with regard to both compressor [2–4] and turbine blades [5, 6]. It is currently explored toovercome the strict requirement of flutter avoidance by flutter tolerance, where a dynamic instability is permitted aslong as the flutter-induced vibrations remain bounded at an acceptably low level. This change in design paradigm hasthe potential to further improve aerodynamic performance, and thus to reduce emissions and save resources.In the presence of a dynamic instability, linear theory only predicts the unbounded exponential growth of the vibrations(divergence). A saturation and, thus, bounded self-excited oscillations are a priori a nonlinear phenomenon. Nonlin-ear force-deformation relations, in principle, can have various physical origins. For large vibrations, geometric andmaterial nonlinearity may become important on the structural side. Aerodynamic forces may depend nonlinearly onthe vibration amplitude, e.g. , in the case of separating and transonic flows. Finally, the contact interactions occurringat mechanical joints (blade attachment in disk, underplatform friction dampers, interlocked tip shrouds, etc.) give riseto nonlinear contact forces. This work focusses on nonlinear contact forces.In the steady state, flutter-induced vibrations typically take the form of a periodic oscillation, denoted as Limit CycleOscillation (LCO) [7, 8]. This is also the case when multiple modes of vibration are aeroelastically unstable [8].Besides LCOs, quasi-periodic oscillations (or Limit Torus Oscillations) may also occur. Limit Torus Oscillationsare caused by the interaction between multiple unstable modes of vibration, if their frequencies satisfy a certaincombination resonance condition, and nonlinearity is sufficiently strong [9]. The methods presented in this work focuson LCOs, the generalization to Limit Torus Oscillations is left for future research. Over a period of an LCO, the work Δ 𝑊 f < Δ 𝑊 s > Δ 𝑊 f + Δ 𝑊 s = Δ 𝑊 s is dominated by dry frictional dissipationin the joints. An LCO is characterized by its angular frequency 𝜔 , deflection shape 𝝍 , amplitude ˆ 𝑢 and potentialhigher-harmonic content. The accurate prediction of these properties relies on a precise modeling of aerodynamic andstructural forces. For modeling of the structural forces within frictionally coupled bladed disks, we refer to a recent2verview [10]. For modeling of the aerodynamic forces, several simplifying assumptions are commonly made, asdiscussed in the following.Firstly, it is very common to assume the aerodynamic forces as linear in the structural displacement amplitude ˆ 𝑢 . Inthe amplitude range relevant for high-cycle fatigue, this assumption has been supported by experimental [11, 12] andnumerical [13–15] evidence.Secondly, it is often assumed that 𝜔 and 𝝍 correspond to a certain normal vibration mode of the structure in vacuum,obtained under linearized contact boundary conditions in joints (i.e. sticking contact in preloaded interfaces). Conse-quently, these properties can be efficiently computed using finite elements and linear modal analysis. The associatedaerodynamic work per cycle Δ 𝑊 f can then be determined for each mode of interest using computational fluid dynamics.This leads to the conventional energy method , where frequency and shape of the vibration are assumed as descriedabove, and the amplitude is determined from the requirement that supplied and dissipated works per cycle cancel eachother ( Δ 𝑊 f + Δ 𝑊 s = 𝝍 and frequency 𝜔 of the LCO. However, the nonlinear contact boundary conditions can easily account for morethan 10% frequency shift and significantly alter the deflection shape [19], as compared to the linearized case. Giventhat the reduced frequency and the deflection shape are the driving determinants for flutter, this strategy has a strictlylimited range of validity. Indeed, the sign of the aerodynamic work can be affected by the boundary conditions at thecontact interfaces. For instance, Leyes et al. [20] associate a blade failure in an aircraft engine to a flutter instabilitycaused by the altered boundary conditions at the tip shrouds due to fretting wear. Similarly, Wu et al. [21] demonstratea test case that was stable for sticking and unstable for frictionless sliding tip shroud interfaces.To account for a potential change of the deflection shape 𝝍 and frequency 𝜔 with the amplitude ˆ 𝑢 , the linear aero-dynamic forces can be modeled in terms of influence coefficients in the frequency domain [5, 8, 22, 23]. These areformulated with respect to a selected set of modes. To predict these coefficients, a harmonic deformation in each modeis imposed (again, with infinitesimal amplitude) and the resulting aerodynamic force, more specifically its fundamen-tal Fourier coefficient, is determined using computational fluid dynamics. The force induced by unit deformationin mode 𝑗 , acting on mode 𝑖 is the influence coefficient 𝐺 𝑖 𝑗 . This way, a modal aerodynamic influence coefficientmatrix (MAIM) 𝑮 = [ 𝐺 𝑖 𝑗 ] is set up. The MAIM is a complex-valued matrix; the real parts of the diagonal entriesare associated with the aerodynamic stiffness and the imaginary parts with the aerodynamic damping. The MAIMgenerally depends on the frequency 𝜔 of the imposed vibration. If the actual vibration frequency is accurately known,the frequency dependence can be neglected [24]. For cases where the relevant frequency range is large, a piecewiselinear interpolation was proposed in [25]. Just as with the conventional energy method, an advantage of the methodof influence coefficients is the decoupling of computational fluid dynamics and nonlinear structural dynamics: Oncethe MAIM is set up, one has a closed-form expression for the aerodynamic forces. This can subsequently be usedto compute the LCO by considering the balance with the nonlinear structural forces. As the computational effort forsetting up the MAIM grows quickly with the number of modes, the appropriate choice of a modal basis for describingstructural vibration is a crucial aspect.An alternative to using influence coefficients is the multi-physical dynamic co-simulation of both structural and fluidmechanics, i.e. , a fully-coupled Fluid-Structure Interaction computation (see e.g. [26, 27]). An advantage of thisapproach is that it is straight-forward to account for nonlinear effects in both domains. However, the available methodsrely on numerical time integration and lead to computational efforts prohibitive for the design turbomachinery compo-nents [26, 28]. To overcome the high computational effort of numerical time integration, nonlinear frequency domainmethods such as the Harmonic Balance (HB) method have been developed which are now popular in both structuraldynamics and computational fluid dynamics. Compared to numerical time integration, the computational effort istypically reduced by two to four orders of magnitude [29, 30]. The reasons for the computational advantages of HBare: (a) HB can efficiently deal with stiff differential equations. (b) HB avoids the costly simulation of the usuallylong transient. (c) A few harmonics often suffice to resolve the periodic oscillation. (d) Assuming symmetric behaviorwithin the blade row, the computational domains can be reduced to a reference sector (only one blade or one passage)with appropriate cyclic symmetry boundary conditions. First attempts to co-simulate cascade flutter in the time orfrequency domain are limited to linear [31–34] or single-degree-of-freedom nonlinear [35, 36] structural models.As described above, the current approach to design against flutter relies on linear vibration theory. On the one hand,this strategy is too conservative, as this excludes the design space where flutter would lead to vibrations of tolerablelevel. On the other hand, the current approach is not conservative enough, as the case of nonlinear instability is ignored;3.e., a situation where the static equilibrium is aeroelastically stable, but an instability occurs for (perhaps rather smallbut) finite vibration amplitudes due to nonlinear effects. There is thus a need for simulation methods that predict thenonlinear Fluid-Structure Interaction. The purpose of this work is to develop efficient frequency-domain methods forcomputing flutter-induced LCOs of bladed disks, which are capable of accounting for both structural and aerodynamicnonlinear effects. In Section 2, the model of the considered benchmark is presented, along with the frequency-domainformulation of the governing equations in either domain. Then, the conventional energy method is refined to take intoaccount structural nonlinearities and their effect on the aerodynamic damping (Section 3). The fully-coupled solver ispresented in Section 4. The computational performance of the methods is then assessed for the case of a stable LCOand an unstable LCO (stability limit) in Section 5. The article ends with the conclusions in Section 6.
2. Benchmark and modeling
A model of a low pressure turbine bladed disk with interlocked tip shrouds is considered. The structural finiteelement model of one sector is shown in Fig. 1; the whole bladed disk contains 50 sectors. The structural modeling,briefly described in the following, represents the current state of the art, and is described in more detail in [10].Geometric and material properties are listed in Tab. 1. Mistuning is known for its potential to significantly distort themodal deflection shapes and for its tendency to stabilize flutter [11, 24]. On the other hand, for the case of stronginter-sector coupling, here via tip shrouds, the effect of small mistuning is known to be negligible [37]. In this work,the system is assumed as rotationally periodic; i.e. , each sector has identical properties. It is further assumed that theperiodic response of the system inherits the symmetry suggested by the system; i.e. , the response is assumed to takethe form of a traveling wave, where each sector undergoes the same oscillation, however, with a constant time shiftbetween neighboring sectors. This permits to reduce the problem domain to only a reference sector with appropriatetime-shift boundary conditions at the cyclic sector boundaries. cyclicboundary sensornodecontact areafixedboundary
Figure 1: Structural finite element model of the sector of a bladed disk (kindly provided by MTU Aero Engines AG)
Nonlinear dry frictional interactions are taken into account between the tip shrouds of neighboring blades. Inaccordance with the cyclic symmetry assumption, the contact is defined under consideration of the time shift betweenthe interfaces on either side of the sector. In the tangential plane, dry friction is modeled in terms of the commonelastic Coulomb law, which can be written in differential form asd 𝒑 t = (cid:40) 𝑘 t d 𝒈 t : || 𝒑 t + 𝑘 t d 𝒈 t || ≤ 𝜇𝑝 n 𝑘 t d 𝒈 t | | 𝒑 t + 𝑘 t d 𝒈 t | | : || 𝒑 t + 𝑘 t d 𝒈 t || > 𝜇𝑝 n . (1)4umber of blades 50Blade span 86 mmBlade chord 22 mmYoungs modulus blade + disk 110 GPaPoisson ratio blade + disk 0 . . × kg m − Youngs modulus shroud 210 GPaPoisson ratio shroud 0 . . × kg m − Radius hub 134 mmRadius tip 221 mm
Table 1: Geometric and material properties of the benchmark model.
Herein, 𝒑 t and 𝑝 n are tangential traction vector and normal pressure, 𝒈 t is the tangential gap vector, and 𝑘 t is thetangential stiffness per area value. A constant friction coefficient 𝜇 is considered for both static and dynamic friction.The normal pressure 𝑝 n is assumed as time-constant and uniformly distributed within the contact interface. Hence,the contact is always closed in the normal direction; liftoff is prohibited. The specific values of 𝜇𝑝 n and 𝑘 t are definedlater. The contact law in Eq. (1) is spatially discretized using six C2D3 surface elements with a total of 8 nodes (eachside), conform with the underlying solid elements. The elements and nodes of the opposing sides coincide in the cyclicsense. The nodes are used as integration points. Dry frictional damping in the mechanical joints is the only source ofstructural damping considered in this work; i.e. , material damping is neglected.For convenience, a coordinate transform to the relative tangential displacements at each contact node pair is carried out(two orthogonal directions per node pair). Assuming small vibrations around the equilibrium position, the structuralelastic and inertia forces are linear in the structural displacement. To reduce the mathematical model order, theCraig-Bampton method is applied [38]. Accordingly, the vector of nodal displacements of the finite element model, 𝒖 fem is approximated as a linear combination, 𝒖 fem = 𝑻 cb 𝒖 s with 𝑻 cb = (cid:20) 𝑰 𝚿 𝚽 (cid:21) , (2)of a set of component modes (the columns of the matrix 𝑻 cb ) and associated coefficients (the elements of the vector 𝒖 s ).The set of component modes is formed by static constraint modes and fixed-interface normal modes. Here, the staticconstraint modes are associated with the relative tangential displacements of each contact node pair. 𝚿 representsthe static displacement of the remaining nodal degrees of freedom in response to a unit displacement of each relativetangential coordinate. The fixed-interface normal modes are obtained by fixing all relative tangential coordinates(rigidly sticking contact); the deflection shapes associated with the lowest-frequency normal modes are collected ascolumns in the matrix 𝚽 . Based on the results of a convergence study, the 5 lowest-frequency normal modes wereretained. Consequently, the dimension of 𝒖 s is 𝑁 s = × + =
21. Substituting the approximation in Eq. (2) intothe equation of motion and requiring that the residual is orthogonal with respect to the component modes, yields areduced order model of dimension 21. The mass and stiffness matrices, 𝑴 cb and 𝑲 cb , in the reduced order model arerelated to the full-order finite element matrices via 𝑴 cb = 𝑻 cb H 𝑴 fem 𝑻 cb and 𝑲 cb = 𝑻 cb H 𝑲 fem 𝑻 cb , and (aerodynamicand contact) forces are projected as 𝒇 cb = 𝑻 cb H 𝒇 fem . The superscript H denotes the Hermitian (or complex-conjugatetranspose). It is emphasized that the matrix 𝑻 cb depends on the structural properties only. The structural vibration isnow described by the reduced vector of generalized coordinates 𝒖 s , which is the solution of the coupled fluid-structureproblem.As mentioned before, we focus on LCOs; i.e. , time-periodic responses. The periodic solution is approximated with atruncated Fourier series, 𝒖 s ( 𝑡 ) ≈ (cid:60) (cid:40) 𝐻 s ∑︁ 𝑘 = e i 𝑘 𝜔𝑡 ˆ 𝒖 s 𝑘 (cid:41) , (3)5here 𝑡 denotes time, i = √− 𝐻 s is the harmonic truncation order, 𝜔 is the fundamentaloscillation frequency and ˆ 𝒖 s 𝑘 are the Fourier coefficients. Applying the Harmonic Balance technique [39] to thereduced order model, the algebraic equation system, 𝑹 s ( ˆ 𝒖 s , 𝜔, ˆ 𝒇 a ) : = 𝑺 ( 𝜔 ) ˆ 𝒖 s + ˆ 𝒇 c ( ˆ 𝒖 s ) − ˆ 𝒇 a = , (4)is obtained. The Fourier coefficients are here stacked as ˆ 𝒖 s = [ ˆ 𝒖 s0 ; ˆ 𝒖 s1 ; . . . ; ˆ 𝒖 s 𝐻 ] , where ; denotes vertical concatenation. 𝑺 ( 𝜔 ) is the block diagonal dynamic stiffness matrix, which describes the linear-elastic and linear inertia forces andthus depends on 𝑴 cb and 𝑲 cb . ˆ 𝒇 c and ˆ 𝒇 a are the vectors of Fourier coefficients of the contact and aerodynamic forces,stacked analogous to ˆ 𝒖 s . ˆ 𝒇 c ( ˆ 𝒖 s ) is computed by sampling 𝒖 s at equidistant time instants, to determine the steadyhysteresis cycle of the contact forces in discrete time, and applying the discrete Fourier transform to determine theFourier coefficients ˆ 𝒇 c . This procedure is well-known as alternating frequency-time scheme. The compressible unsteady Reynolds-averaged-Navier-Stokes equations are formulated with respect to the con-servative variables 𝒖 a = [ 𝜌, 𝜌𝑢, 𝜌𝑣, 𝜌𝑤, 𝜌𝐸 ] 𝑇 (fluid density, 𝑥 -, 𝑦 -, 𝑧 - component of linear momentum, energy pervolume). The equations are closed with the 𝑘 - 𝜔 turbulence model [40], the ideal gas law and the Sutherland law forthe molecular viscosity. The problem domain is discretized with a finite volume mesh containing 𝑁 cell = , 𝒖 a ( 𝑡 ) , which are approximated analogously to Eq. (3), but withtruncation order 𝐻 a . For each cell 𝑗 , harmonic balance yields the equation [41]: 𝑹 a 𝑘, 𝑗 ( ˆ 𝒖 a , 𝜔, ˆ 𝝌 ) : = i 𝑘𝜔 ( (cid:100) 𝑉 𝒖 a ) 𝑘, 𝑗 + (cid:98) 𝑹 𝑘, 𝑗 ( ˆ 𝒖 a , ˆ 𝝌 , 𝜔 ) = 𝑘 = , . . . , 𝐻 𝑗 = , . . . , 𝑁 cell . (5)Herein, ˆ 𝝌 represents the mesh location, which is dynamically deformed under consideration of a Arbitrary Lagrangian-Eulerian formulation. In accordance with the assumption of a traveling-wave-type response, only a single passage ofthe blade row is considered in the fluid model. At the cyclic sector boundaries, phase lag boundary conditions areimposed. At the inlet and outlet, two-dimensional non-reflecting boundary conditions are considered [42]. The fluid mesh must deform compatibly to the structural deformation at the blade surface. To this end, the structuralnodal displacements 𝒖 fem are mapped onto the fluid mesh boundary Γ Blade Surface , 𝒖 ( 𝒙 ) interpolated = Φ ( 𝒖 fem , 𝒙 ) = Φ ( 𝑻 cb 𝒖 s , 𝒙 ) ∀ 𝒙 ∈ Γ Blade Surface , (6)using a bilinear spatial interpolation Φ with respect to the spatial coordinate 𝒙 . The mesh deformation is treated as alinear superposition, 𝝌 ( 𝒖 s ) = 𝝌 + [ 𝝌 𝝌 ... 𝝌 𝑁 ] 𝒖 s , (7)of the mesh deformations 𝝌 ℓ with 𝑁 s ≥ ℓ ≥ ℓ -th generalized structuralcoordinate, and the coordinates 𝒖 s used as weights. 𝝌 is the location of the undeformed mesh. 𝝌 ℓ with ℓ ≥ ∇ · ( 𝜇 ( 𝒙 )∇ 𝝌 ℓ ( 𝒙 )) = ∀ 𝒙 ∈ Ω Fluid Domain , (8) 𝝌 ℓ ( 𝒙 ) = Φ ( 𝑻 cb 𝒆 ℓ , 𝒙 ) ∀ 𝒙 ∈ Γ Blade Surface , (9)where the vector 𝒆 ℓ is the ℓ -th unit vector (containing only zeros, except for the ℓ -th element which is equal to one), ∇ denotes the gradient, · denotes the inner product, and 𝜇 ( 𝒙 ) corresponds to the equivalent elastic modulus. In eachcell, 𝜇 ( 𝒙 ) is made proportional to the inverse of the cell volume which makes sure that small cells, e.g. in a boundarylayer, are not sheared excessively. Eqs. (8)-(9) are efficiently solved with the generalized minimal residual method.Having an explicit expression for the mesh deformation, obtained in a pre-processing step, substantially reduces thecomputation effort compared to computing the mesh deformation simultaneously with the solution of the governingequations of the fluid. The described procedure preserved good mesh quality throughout the amplitude range consideredin this work. By integration of the fluid pressure and traction over the blade surface, nodal forces consistent withthe underlying finite elements are obtained. Finally, these are projected onto the component modes to obtain thegeneralized aerodynamic forces 𝒇 a . 6 .3. Damping and Excitation To understand the occurrence of limit cycles and to assess their stability, it is useful to analyze the work performedby aerodynamic forces on the blade surface and the work dissipated by dry friction, per cycle of oscillation. Thereforethe damping ratio (or loss factor) is used, which is defined as [44], 𝐷 = Δ 𝑊 𝜋𝐸 p , max , (10)where 𝐸 p , max is the maximum potential energy and Δ 𝑊 the work per cycle. The aerodynamic work per cycle isobtained by integrating the power due to aerodynamic forces over one oscillation period, Δ 𝑊 a = ∫ 𝜋 / 𝜔 (cid:164) 𝒖 s T 𝒇 a d 𝑡 ,where (cid:164) (cid:3) denotes derivative with respect to time 𝑡 . The work dissipated by dry friction is obtained analogously by Δ 𝑊 s = ∫ 𝜋 / 𝜔 (cid:164) 𝒖 s T 𝒇 c d 𝑡 .A negative value of the aerodynamic damping 𝐷 a indicates an aerodynamic self-excitation. Elastic sticking and frictionless sliding contact conditions are linear limit cases of the elastic dry friction nonlinearityreached for sufficiently small and (asymptotically) for very large vibration amplitudes, respectively. For very smallvibration amplitudes the norm of the frictional traction will not exceed 𝜇𝑝 n . In this case the contact elements willbehave as linear springs (which corresponds to elastic sticking). For the case of very large amplitudes the contactelements slide most of the time during the vibration cycle. The norm of the frictional traction is limited by the finitevalue of 𝜇𝑝 n , and thus become negligible compared to the stresses within the body, which further increase linearlywith the vibration amplitude. Hence, the behavior for large vibrations asymptotically approaches the behavior forfrictionless sliding.The normal modes of a rotationally periodic structure have the form of sinusoidal waves, standing or traveling aroundthe circumference. The wave number is characterized by the number of nodal diameters (ND), which are diametrallines of material points with zero deflection. ND = = − − Fig. 3 shows an instantaneous snapshot of the cyclic symmetric flow field around the blade. The contours areplotted at constant radius. The lower contour shows the position of a shock and the upper contour visualizes thegeneration of entropy due to the flow separation in the latter part of the suction side. The flow separation coversroughly the upper two third of the blade span.In the following, the vibrational amplitude denoted by ˆ 𝑢 Sensor is the amplitude of the (circumferential) 𝑌 -componentof the sensor node depicted in Fig. 1. This coordinate is chosen since the investigated mode with ND -4 deformsconsiderably in this direction at the selected sensor node. For this number nodal diameters, the tip shroud displacementis relatively large, such that the displacement of the sensor node (located near the tip shroud) properly reflects theglobal vibration level of the blade. In the following plots ˆ 𝑢 Sensor is normalized by the blade span 𝐿 Span .To analyze the amplitude-dependence of the aerodynamic forces, the aerodynamic damping is computed for differentvibration amplitudes of the mode shape 𝝍 stick and natural frequency 𝜔 stick obtained with sticking contact conditions.7 tableunstable -20 -10 0 10 20 ND el. stickingsliding (a) -8 -7 -6 -5 -4 -3 ND -0.02-0.0100.010.020.030.04 el. stickingsliding (b) Figure 2: Aerodynamic damping ratio vs. ND, configuration 2: (a) Overview; (b) Detail.Figure 3: Instantaneous flow field around the blade (mach and entropy contours).
In Fig. 4 a nonlinear behavior of the aerodynamic response can be identified, however the change in damping due tothe amplitude is not significant compared to the impact of the ND or the contact conditions. For amplitudes larger thanˆ 𝑢 Sensor / 𝐿 Span = × − the fluid solver failed to converge. However, the computed range covers the LCO amplitudeswhich will later be analyzed. The modes associated with several other NDs showed similar behavior under variousflow conditions for the given test case.It should be emphasized that a linear dependence of the aerodynamic forces on the vibration amplitude should notbe confused with a linear dependence on the conservative flow variables. The highly unsteady flow phenomenacan certainly be captured only with nonlinear computational fluid dynamics. Even if there is no significant directamplitude dependence (for fixed deflection shape and frequency, 𝝍 stick , 𝜔 stick ), the nonlinear contact conditions willinduce an amplitude dependent natural frequency 𝜔 nl ( ˆ 𝑢 ) and deflection shape 𝝍 nl ( ˆ 𝑢 ) , which in turn generates an indirect amplitude dependence of the aerodynamic damping and stiffness.
3. Energy methods
The basic idea behind energy methods is to first determine the frequency 𝜔 and the deflection shape 𝝍 , and in asecond step determine the amplitude ˆ 𝑢 by the requirement that dissipated and supplied energies per cycle cancel eachother. The result of such an analysis is illustrated in Fig. 5(a). The right plot depicts Δ 𝑊 s and − Δ 𝑊 a as function ofthe amplitude ˆ 𝑢 ; at the intersection points, thus, Δ 𝑊 a + Δ 𝑊 s = Δ 𝑊 a and Δ 𝑊 s on the amplitude ˆ 𝑢 , one can alsoinfer the stability of a given limit cycle: A positive slope at the zero crossing of Δ 𝑊 a + Δ 𝑊 s means that for increasedamplitude, the dissipated work exceeds the supplied one. Thus, the system returns to the limit cycle in the presenceof a small perturbation; the limit cycle is stable. If, in contrast, the slope is negative, an arbitrarily small perturbation8 -4 -3 -5 Figure 4: Aerodynamic damping ratio vs. amplitude for ND -4 mode under sticking contact conditions, configuration 2. (a) (b)
Figure 5: Example of a stable LCO and a stability limit: (a) Damping ratio vs. amplitude; (b) Work vs. amplitude. would drive the amplitude away from this point; the limit cycle is unstable.In the example shown in Fig. 5, two limit cycles can be identified. This is a typical situation under aerodynamicself-excitation and dry frictional dissipation. Assuming that the aerodynamic forces depend linearly on the structuraldisplacements, and the mode shape and frequency are fixed, the aerodynamic damping ratio is constant and Δ 𝑊 a grows quadratically with the amplitude. For small vibration amplitudes, the contact is sticking, and consequently Δ 𝑊 s is zero. For sufficiently large amplitudes, sliding friction occurs. For constant normal load, the area enclosed in theforce-displacement hysteresis cycle, which corresponds to Δ 𝑊 s , grows linearly beyond this point. The left intersectioncorresponds to a stable, the right to an unstable limit cycle. As there is no further limit cycle at higher amplitudes (forthe given model), the unstable limit cycle can also be viewed as a stability limit .The available energy methods differ by the simplifying assumptions made for the determination of 𝜔 and 𝝍 . Traditionally, the normal modes are computed for the structure in vacuum (neglecting the aerodynamic forces)and under linearized contact conditions (here: sticking shrouds). More specifically, 𝝍 stick and 𝜔 stick are obtainedfrom an eigenvalue analysis based on the finite element structural model, ( 𝑲 fem + 𝑲 contact − ( 𝜔 stick ) 𝑴 fem ) 𝝍 stick = .For a given mode of interest, subsequently, the work dissipated by structural forces, Δ 𝑊 s and that provided byaerodynamic forces, Δ 𝑊 a , are determined, as functions of the modal amplitude ˆ 𝑢 . In this step, a harmonic displacement 𝒖 fem ( 𝑡 ) = (cid:60) (cid:110) 𝝍 stick ˆ 𝑢𝑒 i 𝜔 stick 𝑡 (cid:111) with amplitude ˆ 𝑢 is imposed. Without loss of generality, a positive real-valued amplitudeˆ 𝑢 is assumed in the following. 9or blade vibration problems, it is common to assume that structural dissipation is dominated by the dry frictionalcontact forces in the mechanical joints, 𝒇 c , fem . Indeed, material damping is neglected in this study. Under imposedharmonic displacement 𝒖 fem ( 𝑡 ) = 𝒖 fem ( 𝑡 + 𝑇 ) with period 𝑇 = 𝜋 / 𝜔 stick , the steady-state contact forces will also beperiodic. The work dissipated by these forces per cycle of vibration is then Δ 𝑊 s ( ˆ 𝑢 ) = ∫ ( 𝑇 ) (cid:16) (cid:164) 𝒖 fem (cid:17) T 𝒇 c , fem d 𝑡 = ∫ ( 𝑇 ) (cid:60) (cid:110) i 𝜔 stick 𝝍 stick ˆ 𝑢 e i 𝜔 stick 𝑡 (cid:111) T (cid:60) (cid:110) ˆ 𝒇 c , fem1 e i 𝜔 stick 𝑡 (cid:111) d 𝑡 = ∫ ( 𝜋 ) (cid:60) (cid:8) i 𝝍 stick ˆ 𝑢 e i 𝜏 (cid:9) T (cid:60) (cid:110) ˆ 𝒇 c , fem1 e i 𝜏 (cid:111) d 𝜏 = ∫ ( 𝜋 ) (cid:18) i 𝝍 stick ˆ 𝑢 e i 𝜏 + i 𝝍 stick ˆ 𝑢 e − i 𝜏 (cid:19) T (cid:18) ˆ 𝒇 c , fem1 e i 𝜏 + ˆ 𝒇 c , fem1 e − i 𝜏 (cid:19) d 𝜏 = (cid:18) − i (cid:16) 𝝍 stick (cid:17) H ˆ 𝒇 c , fem1 / + i (cid:16) 𝝍 stick (cid:17) T ˆ 𝒇 c , fem1 / (cid:19) ˆ 𝑢 ∫ ( 𝜋 ) d 𝜏 = (cid:60) (cid:26) − i 𝜋 (cid:16) 𝝍 stick (cid:17) H ˆ 𝒇 c , fem1 (cid:27) ˆ 𝑢 . (11)Herein, 𝜏 = 𝜔 stick 𝑡 , (cid:3) denotes the complex conjugate and ˆ 𝒇 c , fem1 is the complex fundamental Fourier coefficient theperiodic contact forces. It should be remarked that all other harmonics of the periodic contact forces are orthogonal tothe purely harmonic velocity. Since the relative displacement is imposed, the dissipated work can be easily determinedseparately for each contact element, and then the sum is taken over all elements. For given modal deflection shape Δ 𝑊 s is only a function of the modal amplitude ˆ 𝑢 .For the imposed harmonic structural displacement, the flow field is simulated using computational fluid dynamics, seee.g. [45, 46]. The subsequent spatial integration over the blade surface yields the nodal aerodynamic force vector 𝒇 a , fem . The aerodynamic work per cycle is determined in analogy to the work dissipated by dry friction, Δ 𝑊 a ( ˆ 𝑢 ) = (cid:60) (cid:26) − i 𝜋 (cid:16) 𝝍 stick (cid:17) H ˆ 𝒇 a , fem1 (cid:27) ˆ 𝑢 = (cid:60) (cid:8) − i 𝜋 ˆ 𝑓 a (cid:9) ˆ 𝑢 = 𝜋 (cid:61) { 𝐺 } ˆ 𝑢 . (12)In the last step of this equation, it is assumed that the vibration-induced aerodynamic force is linear with respect to themodal amplitude, so that the complex amplitude of the modal aerodynamic force can be expressed as ˆ 𝑓 a = 𝐺 ˆ 𝑢 , withthe complex modal aerodynamic influence coefficient 𝐺 . Like the work dissipated by dry friction, the aerodynamicwork per cycle is a function of ˆ 𝑢 only for given modal deflection shape and natural frequency. In the following, we propose a refinement of the conventional energy method, which takes into account the nonlinearamplitude-dependence of mode shape, 𝝍 nl ( ˆ 𝑢 ) , and natural frequency, 𝜔 nl ( ˆ 𝑢 ) . To determine these, a nonlinear modalanalysis is carried out in accordance with the Extended Periodic Motion Concept [47]. Again, Harmonic Balance isused to compute the periodic oscillations. This yields the algebraic equation system, 𝑺 𝑘 ( 𝜔 nl ) ˆ 𝒖 s 𝑘 + ˆ 𝒇 c 𝑘 ( ˆ 𝒖 s ) = 𝐷 s 𝜔 nl ( i 𝑘𝜔 nl ) 𝑴 cb 𝑘 ˆ 𝒖 s 𝑘 𝑘 = , , . . . , 𝐻. (13)The term on the right hand side represents the mass-proportional excitation term in accordance with the definition in[47]. Here, the effect of the flow is still neglected (as in the conventional energy method), but the nonlinear boundary10onditions at the contact interfaces are considered. The nonlinear modal analysis starts with a small amplitude in thelinear regime in which only sticking (but no sliding) friction occurs. By numerical path continuation, the amplitude-dependent mode shapes 𝝍 nl ( ˆ 𝑢 ) and natural frequency 𝜔 nl ( ˆ 𝑢 ) are obtained. The mode shape is mass normalized,i.e. , ˆ 𝒖 s1 = 𝝍 nl ˆ 𝑢 with (cid:0) 𝝍 nl (cid:1) H 𝑴 cb 𝝍 nl as in the linear case. The modal deflection shape in the nodal displacementcoordinates of the initial finite element is recovered according to ˆ 𝒖 fem = 𝑻 cb 𝝍 nl ( ˆ 𝑢 ) ˆ 𝑢 . The work dissipated by dryfriction immediately follows from Eq. (10).To account for the change of the mode shape, the aerodynamic force ˆ 𝒇 a is approximated using the aerodynamicinfluence coefficient matrix 𝑮 with respect to the component modes,ˆ 𝒇 a ( ˆ 𝑢 ) = 𝑮 (cid:16) 𝜔 nl ( ˆ 𝑢 ) (cid:17) 𝝍 nl ( ˆ 𝑢 ) ˆ 𝑢 . (14) 𝑮 is set up column-wise, 𝑮 = [ 𝑮 . . . 𝑮 𝑁 ] , where 𝑮 𝑖 = 𝑻 cb H ˆ 𝒇 a , fem1 (cid:0) 𝑻 cb 𝒆 𝑖 (cid:1) for 𝑖 = , . . . , 𝑁 . Thus, the fundamentalFourier coefficient of the aerodynamic forces must be determined, obtained under a harmonic structural displacementin the form of every component mode (index 𝑖 ). To account for the change of the natural frequency, 𝑮 is linearlyinterpolated as 𝑮 ( 𝜔 ) = 𝑮 slip + ( 𝑮 stick − 𝑮 slip )( 𝜔 − 𝜔 slip )/( 𝜔 stick − 𝜔 slip ) . Herein, 𝜔 slip and 𝜔 stick are the naturalfrequencies under the frictionless sliding and sticking linear limit cases, and 𝑮 slip , 𝑮 stick are obtained by imposing thestructural displacement with the corresponding frequency. The aerodynamic work per cycle is then given by Δ 𝑊 a ( ˆ 𝑢 ) = (cid:60) (cid:110) − i 𝜋 (cid:16) 𝝍 nl ( ˆ 𝑢 ) (cid:17) H 𝑮 (cid:16) 𝜔 nl ( ˆ 𝑢 ) (cid:17) 𝝍 nl ( ˆ 𝑢 ) (cid:111) ˆ 𝑢 . (15)An alternative to using the influence coefficient matrix in Eq. (14), the aerodynamic forces can be computed byimposing the deflection shape 𝝍 nl ( ˆ 𝑢 ) and frequency 𝜔 nl ( ˆ 𝑢 ) directly, for specific values of the amplitude ˆ 𝑢 . Inpreliminary investigations, it was found that the results of both alternatives are in very good agreement in the relevantamplitude and frequency range, in spite of the relatively large frequency shift by about 19%. It should be remarkedthat the situation might be different for another benchmark. The analysis of 𝑮 ( 𝜔 ) can provide deeper insight intothe aerodynamic interaction of certain modes of vibration. The direct evaluation of ˆ 𝒇 a1 using 𝝍 nl ( ˆ 𝑢 ) and 𝜔 nl ( ˆ 𝑢 ) can indeed be more computationally feasible, as discussed later, and is not a priori limited to linear amplitude- andfrequency-dependence of the aerodynamic forces.
4. Fully-coupled solver
The algorithm of the fully-coupled solver is illustrated in Fig. 6. A serial coupling strategy is followed, wherethe structural and the fluid subproblems are solved in an alternating way. An initial guess for the frequency 𝜔 ( ) and the structural deformation ˆ 𝒖 s ( ) (in terms of Fourier coefficients of the generalized coordinates) is required, asdefined later. The outer iteration loop starts by updating the flow mesh deformation in accordance with the periodicstructural displacement, as described in Section 2. Next, the harmonic balance residual of the fluid domain (Eq. (5))is minimized by iteratively improving the flow variables ˆ 𝒖 a . In this step, the frequency, 𝜔 , and the Fourier coefficientsof the structural deformation, ˆ 𝒖 s , are held constant. Once the residual norm falls below a specified tolerance, theaerodynamic force is determined. Subsequently, the harmonic balance residual of the structural domain (Eq. (4)) isminimized by iteratively improving both the frequency 𝜔 and the structural variables ˆ 𝒖 s . In this step, the aerodynamicforce, ˆ 𝒇 a , is considered as input. As the residuals in both domains are always reduced below the specified tolerance,the current estimate is considered as solution of the coupled problem if it does not change between two iterations(according to a given tolerance).An important benefit of the serial coupling strategy is that it is not intrusive, in the sense that existing codes foreither domain can be used, and individual single-domain solvers can be employed that are optimized for the individualsubproblems. For the structural model, a Newton-like solver with analytical gradients is used, as implemented in theopen source Harmonic Balance tool NLvib [39]. The fluid model is implemented in the tool TRACE [48, 49], wherethe Harmonic Balance equations are solved with a pseudo time stepping scheme.In contrast to energy methods, it is straight-forward to account for all nonlinear effects that can be modeled with theunderlying fluid and structure solvers. Hence, the fully-coupled solver is more versatile and has the potential for amuch higher prediction accuracy. It should be emphasized that the vibrational frequency is not a priori known and11 tructure iterate on ˆ 𝒖 s , 𝜔 until (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝑹 s (cid:16) ˆ 𝒖 s , 𝜔, ˆ 𝒇 a ( 𝑚 ) (cid:17)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < 𝜖 𝑠 (cid:12)(cid:12) 𝜔 ( 𝑚 ) − 𝜔 ( 𝑚 − ) (cid:12)(cid:12) < 𝜖 𝜔 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ 𝒖 s ( 𝑚 ) − ˆ 𝒖 s ( 𝑚 − ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < 𝜖 𝑢 flow mesh ˆ 𝝌 ( 𝑚 ) = [ 𝝌 𝝌 . . . 𝝌 𝑁 ] ˆ 𝒖 s ( 𝑚 ) 𝑚 ← 𝑚 + flow iterate on ˆ 𝒖 a until (cid:12)(cid:12)(cid:12)(cid:12) 𝑹 a (cid:0) ˆ 𝒖 a , 𝜔 ( 𝑚 − ) , ˆ 𝝌 ( 𝑚 − ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) < 𝜖 𝑎 aerodynamic force ˆ 𝒇 a ( 𝑚 ) = ˆ 𝒇 a (cid:16) ˆ 𝒖 a ( 𝑚 − ) , 𝜔 ( 𝑚 − ) (cid:17) solvedinitialization ˆ 𝒖 s ( ) , 𝜔 ( ) no yes Figure 6: The algorithm scheme for the FSI frequency domain solver. considered as part of the solution. The more complicated situation of a not locked-in vibration of the coupled system(as it can occur in case of vortex shedding or acoustic resonance), however, must still be investigated and is not part ofthe paper.A key to the computational efficiency of the coupled solver is that the data exchanged between the respective single-domain solvers does not require significant amount of memory. Only Fourier coefficients of generalized coordinatesand forces plus the oscillation frequency are exchanged, so only ( 𝐻 + ) 𝑁 + 𝐻 = min ( 𝐻 s , 𝐻 a ) real numbersmust be exchanged. In this work, with 𝐻 s = 𝐻 a = 𝑁 =
21, 64 real numbers are exchanged between thesingle-domain solvers. This contrasts the common multi-physics simulation tools, which exchange displacements andforces at all surface nodes and possibly many time levels. In the following subsections, a few crucial aspects for thecomputational efficiency and robustness of the full-coupled solver are addressed.
As the system is autonomous, a given periodic solution can be time-shifted arbitrarily and will still remain asolution. In other words, the phase of the solution is arbitrary in the autonomous case. To achieve local uniqueness, itis common practice to impose a phase constraint. Taking into account that the oscillation frequency is an additionalunknown in the autonomous case, this phase constraint is also necessary to make the number of equations equal tothe number of unknowns. Typically, the real or imaginary part of the fundamental Fourier coefficient of a specificcoordinate is commonly set to zero. While this constrains the phase, it does not exclude the trivial solution (staticequilibrium, ˆ 𝒖 s = ). Without appropriate countermeasures, therefore, the structural solver has the tendency to runinto the static equilibrium. To overcome this, we propose a different phase constraint: We set the real part of thefundamental Fourier coefficient of a selected generalized coordinate to a suitable nonzero value (Fig. 7). Since thesolver can still adjust the imaginary part, this does not constrain the physical solution (except if the prescribed real partexceeds the magnitude, which can be easily checked upon solving). We selected the generalized coordinate associatedwith the dominant fixed-interface mode. This is expected to provide great numerical performance, as this coordinatecharacterizes the motion of the structure very well and is not likely to suddenly flip or significantly change its phase12eIm ˆ 𝑢 s 𝑖 Figure 7: Non-trivial phase constraint to avoid trivial solution. relative to the rest of the coordinates during the solution process. It is important to emphasize that the describedprocedure has no impact on the physical solution but is actually necessary to obtain a well posed problem.
Within the structural subproblem, the aerodynamic force ˆ 𝒇 a is modeled as linear with respect to the vibrationamplitude, but always strictly consistent with the force ˆ 𝒇 a ( 𝑚 ) output by the previous run of the fluid solver. This ensuresthat the converged solution is independent of the assumed linearization. However, the linearization greatly speeds upthe overall convergence of the coupled solver as the assumed linearity is at least locally a reasonable approximation. Asmentioned in the introduction, linearity with respect to the vibration amplitude is well-supported by experimental andnumerical evidence, including the results presented in this study (Fig. 4). Two variants of the described linearizationare presented in the following. If the fully-populated influence coefficient matrix is available the aerodynamic force ˆ 𝒇 a in Eq. (4) can be split asˆ 𝒇 a1 ( ˆ 𝒖 s , 𝜔 ) = ˆ 𝒇 a1 , ( 𝑚 ) + 𝑮 ( 𝜔 ) ˆ 𝒖 s1 − 𝑮 (cid:0) 𝜔 ( 𝑚 − ) (cid:1) ˆ 𝒖 s1 , ( 𝑚 − ) ˆ 𝒇 a 𝑘 = ˆ 𝒇 a 𝑘, ( 𝑚 ) ∀ 𝑘 ≠ . (16)Note that ˆ 𝒇 a ( 𝑚 ) is obtained by imposing ˆ 𝒖 s ( 𝑚 − ) and 𝜔 ( 𝑚 − ) in the previous run of the fluid solver (Fig. 6). 𝑮 ( 𝜔 ) ˆ 𝒖 s canbe viewed as the part of the aerodynamic forces explainable by the influence coefficient matrix. A simplified variant of the previous approach is to neglect the frequency dependance of 𝑮 . This seems justifiedif the change in frequency due to nonlinear effects is small. The simplification halves the computational effort inthe pre-processing, since 𝑮 does not have to be computed at two frequencies in order to obtain an interpolation withrespect to 𝜔 . The disadvantage of the above method is that the computation of a fully-populated influence coefficient matrix isrelatively expensive. Another possibility is to linearize the aerodynamic force, in an ad-hoc way, asˆ 𝒇 a = ˆ 𝒇 a ( 𝑚 ) 𝒆 T 𝑖 𝒖 s1 𝒆 T 𝑖 𝒖 s1 , ( 𝑚 − ) . (17)Recall that 𝒆 𝑖 is the unit vector whose 𝑖 -th element is 1 and the rest are zeros. 𝒆 𝑇𝑖 ˆ 𝒖 s1 thus extracts the 𝑖 -th elementof ˆ 𝒖 s1 , which is, again, the fundamental Fourier coefficient of the dominant fixed-interface mode. This method willperform well if the change of the mode shape and thus the change of ˆ 𝒇 a between two iterations of the coupled solveris moderate. 13 .3. Initialization A good initialization of the fully-coupled solver is important. We used the results from the refined energy methodˆ 𝒖 s ( ) = 𝝍 nl ( ˆ 𝑢 lco ) ˆ 𝑢 lco , 𝜔 ( ) = 𝜔 nl ( ˆ 𝑢 lco ) , (18)where ˆ 𝑢 lco is the modal amplitude of the limit cycle predicted by the refined energy method. The phase of ˆ 𝒖 s ( ) ismodified in order to meet the phase constraint illustrated in Fig. 7.
5. Results
If multiple ND modes are aerodynamically self-excited, it is expected that each of these gives rise to a limitcycle [8]. The larger the aerodynamic self-excitation ( − 𝐷 a ) of a given mode, the larger is the corresponding basinof attraction, i.e. , the set of initial conditions from which the trajectories approach the respective limit cycle. Asa representative example, we analyze the fundamental mode with ND −
4, which also has the highest aerodynamicself-excitation in the frictionless sliding limit case shown in Fig. 2.configuration: 1 2inlet total pressure 1 . × Pa 1 . × Pastatic pressure at outflow 9 . × Pa 7 . × Pacontact area 4 .
55 mm .
55 mm contact stiffness 𝑘 t .
19 N mm − .
95 N mm − limit frictional traction 𝜇𝑝 n .
22 N mm − . × − N mm − 𝐷 a (stick) 4 . × − . × − 𝐷 a (slip) * − . × − natural frequency (stick) 437 Hz 520 Hznatural frequency (slip) 418 Hz 418 Hz Table 2: Parameters of the two benchmark configurations and corresponding linear modal properties. *The aerodynamic damping 𝐷 a (slip) ofconfiguration 1 (used for the stable limit cycle) was not computed since it is not relevant. Two configurations are studied in the following, which differ in terms of contact properties and static outflow pressure,as listed in Tab. 2. During a long service life, the normal preload in the interlocked shrouds can decrease significantlydue to creep and wear. This is the motivation to set different normal preloads in the two configurations. In this sense,configuration 1 corresponds to the initially built setting , while configuration 2 corresponds to a setting after a longservice life . First the expected and well-known stable limit cycle is analyzed for configuration 1. As discussed in Section2, the ND − 𝐻 s = 𝐻 f = Fig. 8 shows the structural and aerodynamic damping ratios as function of the vibration amplitude. For sufficientlysmall vibration amplitudes, the contact is always sticking, and thus the structure behaves linear and there is no frictionaldissipation ( 𝐷 s = | ˆ 𝑢 Sensor | ≈ . 𝐿 Span ),dissipative sliding friction starts to occur, which leads to an increase of 𝐷 s . For the given contact stiffness and area,a maximum structural damping ratio of about 4% is reached. The aerodynamic damping ratio 𝐷 a is negative, hasa relatively small magnitude and does not considerably depend on the amplitude, as compared to 𝐷 s . Recall that14 FSI Solver (a)
FSI Solver (b)
Figure 8: Results of the conventional and refined energy methods, and the fully-coupled (FSI) solver, configuration 1: (a) Global view; (b) Detailedview at LCO solution. a limit cycle requires 𝐷 s + 𝐷 a =
0. Hence, a limit cycle occurs at the intersections in Fig. 8(a) where 𝐷 s = − 𝐷 a .Because of the relatively small magnitude of 𝐷 a , the aerodynamic self-excitation is cancelled by friction dampingalready just shortly beyond the onset of sliding friction. Consequently, the oscillation frequency at the limit cycle isvery close to that of the sticking limit case. Indeed, the difference in frequency is only 0 .
003 % here. Therefore, thefrequency-dependence in the refined energy methods was neglected in this case. In Tab. 3 the results of the energymethods and the fully-coupled solver are summarized in terms of amplitude and frequency of the stable limit cycle.It can be concluded that all methods are in very good agreement for this particular case. For completeness, also thechange of the vibrational deflection shape is assessed by computing the
𝑀 𝐴𝐶 (modal assurance criterion) defined asthe correlation measure
𝑀 𝐴𝐶 ( 𝝍 , 𝝍 ) = | 𝝍 H1 𝝍 | ( 𝝍 H1 𝝍 )( 𝝍 H2 𝝍 ) . (19)The 𝑀 𝐴𝐶 values shown in Table 3 are all very close to unity. Thus, it is concluded that the deflection shape at thestable limit cycle is practically identical in all methods and agrees well with the modal deflection shape obtained forsticking contact conditions.
The decay of the residuals in either domain is shown in Fig. 9(a) for the first four iterations of the fully-coupledsolver. The residuals decrease with rates typical for the decoupled subproblems. The convergence of the limit cycleamplitude and frequency is depicted in Fig. 9(b). Recall that ˆ 𝑢 s 𝑖 denotes the fundamental Fourier coefficient of thedominant fixed-interface mode. The results almost do not change after the first iteration and the full-coupled solverreaches convergence in accordance with the specified tolerances already after 5 iterations. Both the linearizationinvolving the frequency-constant influence coefficient matrix (with the natural frequency of the sticking contact) andthe dominant-mode linearization were assessed and found to lead to numerically identical results in the same numberof iterations. The simulation took 19 . .
952 Hz 436 .
941 Hz 437 .
032 Hzˆ 𝑢 Sensor / 𝐿 Span . × − . × − . × − 𝑀 𝐴𝐶 ( 𝝍 , 𝝍 stick ) .
999 999 95 0 .
999 996 6
𝑀 𝐴𝐶 ( 𝝍 , 𝝍 nl ) - 1 0 .
999 997 4
Table 3: Summary of stable limit cycle data. Pseudo time step, Newton solver step -8 -7 -6 -5 -4 -3 -20 -10 coupled FSI solver iteration }}}} 𝑚 = 𝑚 = 𝑚 = 𝑚 = (a) (b) Figure 9: Behavior of the fully-coupled solver, configuration 1: (a) Successive decay of the residuals in fluid and structural domain; (b) convergenceof the dominant amplitude and frequency of the limit cycle.
The structural and aerodynamic damping ratios are depicted as function of the vibration amplitude in Fig. 10.In contrast to the case of the stable limit cycle, the unstable limit cycle does not occur close to the linear case withsticking contact conditions. Consequently, the oscillation frequency and the deflection shape differ considerably fromthe modal properties of the sticking limit case. The conventional energy method ignores the amplitude-dependenceof the oscillation frequency and deflection shape (assuming 𝝍 stick , 𝜔 stick for all ˆ 𝑢 ), and thus the aerodynamic dampingratio is invariant according to this method. As the damping ratio is positive in the sticking limit case, the conventionalenergy method does not predict a stability limit, but predicts stable behavior for all amplitudes.In Tab. 4 the results of the energy methods and the fully-coupled solver are summarized in terms of amplitude andfrequency of the stability limit. The results of the refined energy method and the fully-coupled solver agree very well.The deviations of frequency, ˆ 𝑢 Sensor are 0 .
02 %, 0 .
88 % respectively. The
𝑀 𝐴𝐶 values of about 0 .
91 with respectto the sticking limit case show that the vibrational deflection shape is significantly affected by the nonlinear contactinteractions in the tip shrouds. This leads to a considerable influence on the aerodynamic response, cf. − 𝐷 a ( 𝜔 nl , 𝝍 nl ) and − 𝐷 a ( 𝜔 nl , 𝝍 stick ) in Fig. 11. The vibrational deflection shape predicted by the coupled solver and the refined energyConventional Energy Method Refined Energy Method full-coupled solverFrequency 519 .
661 Hz 418 .
049 Hz 418 .
671 Hzˆ 𝑢 Sensor / 𝐿 Span - 1 . × − . × − 𝑀 𝐴𝐶 ( 𝝍 , 𝝍 stick ) .
910 69 0 .
910 73
𝑀 𝐴𝐶 ( 𝝍 , 𝝍 nl ) - 1 0 .
999 996
Table 4: Summary of stability limit data. method agree very well, as indicated by a
𝑀 𝐴𝐶 value close to unity (right column, bottom row). This means thatthe vibrational deflection shape is mainly driven by structural forces, and can thus be well-represented by a nonlinearmode. Hence, the fluid does not add substantial stiffness (causing a significant frequency shift) or even change thedeflection shape. The aerodynamic self-excitation merely determines the vibration level.16 -4 -2 -50510152025 FSI Solver (a) -4 -2 -0.04-0.0200.020.040.060.08 FSI Solver (b)
Figure 10: Results of the conventional and refined energy methods, and the fully-coupled (FSI) solver, configuration 2: (a) Global view; (b) Detailedview.
To further study the influence of the amplitude-dependent modal deflection shape 𝝍 nl and frequency 𝜔 nl , therefined energy method is applied with either the frequency or modal deflection shape held constant, denoted by − 𝐷 a ( 𝜔 stick , 𝝍 nl ) and − 𝐷 a ( 𝜔 nl , 𝝍 stick ) , respectively. Both results are plotted in Fig. 11. For larger amplitudes, there isa noticeable difference to − 𝐷 a ( 𝜔 nl , 𝝍 nl ) which shows that both effects are important to determine the accurate flowresponse for this case. Here, the frequency dependency is mainly responsible for the aerodynamic behavior changingfrom damping to self-excitation. This appears plausible since a decrease in vibrational frequency leads to a decreasein reduced frequency which can trigger aerodynamic instability [1]. -4 -2 -0.04-0.0200.020.040.060.08 FSI Solver
Figure 11: Results of the refined energy method (ND -4) with constant frequency or constant modal deflection shape.
The convergence behavior of the fully-coupled solver is illustrated in Fig. 12(a) and Fig. 12(b). Again, thesolver converges quickly, but more iterations are needed for the same tolerances compared to the previous benchmarkconfiguration. In this configuration, the type of linearization used for the representation of the aerodynamic forceswithin the structural subproblem (cf. Section 4.2) was found to have an important influence. The dominating modelinearization performed well, whereas the frequency-dependent influence-coefficient-matrix linearization did not leadto convergence in a reasonable number of iterations. The behavior could be improved by introducing a relaxationfactor of 0 . 𝒖 s ( 𝑚 ) , see the results in Fig. A.1. The coupled solver required 83 . Pseudo time step, Newton solver step -8 -7 -6 -5 -4 -3 -30 -20 -10 coupled FSI solver iteration }}}} 𝑚 = 𝑚 = 𝑚 = 𝑚 = (a) (b) Figure 12: Behavior of the fully-coupled solver, configuration 2: (a) Successive decay of the residuals in fluid and structural domain; (b) convergenceof the dominant amplitude and frequency of the limit cycle. the dominating-mode linearization and the influence-coefficient-matrix linearization (with the described relaxation)respectively.
The computational costs for solving the structural subproblem was practically negligible in the considered bench-mark. For instance, the nonlinear modal analysis took only about 30 seconds. For the refined energy method, thusthe computational effort is driven by the computation of the influence coefficient matrix 𝑮 . This effort scales withthe number of generalized coordinates, 𝑁 , of the reduced structural model (here 𝑁 = 𝑮 ,one flow simulation is required which took approximately 5 .
25 hours on 9 cores. This sums up to 110 .
25 hourson 9 cores. For the linear interpolation with respect to frequency, 𝑮 must be computed for two frequencies, whichdoubles the computational effort. The advantage of this approach is that one obtains a closed-form expression forthe aerodynamic forces, which can be useful to analyze the aerodynamic damping (and stiffness) for a large numberof different vibrational deflection shapes and frequencies (in the considered frequency range). An alternative to thisapproach is to directly compute the aerodynamic forces for a given structural vibration, and to determine the resultingaerodynamic damping. The effort for such a computation is about the same as that for computing a single column of 𝑮 . Thus, this alternative becomes computationally more attractive than using the influence coefficient matrix once thenumber of amplitude levels for which Δ 𝑊 a is of interest falls below the number of generalized coordinates (or twicethis number in the case of the two-point interpolation with respect to frequency).
6. Conclusions
In this work, the conventional energy method has been refined and a fully-coupled frequency-domain method weredeveloped and assessed. These methods are useful to compute flutter-induced limit cycle oscillations of turbomachinerybladed disks. In contrast to the conventional energy method, the new methods account for the amplitude-dependenceof the oscillation frequency and the vibrational deflection shape. The refined energy method is particularly suited toseparate the different effects of contact and aerodynamic forces, which allows deep understanding of the physics ofnonlinear aeroelastic interactions. As the fully-coupled solver follows a serial coupling strategy, available codes canbe used which are optimized for the respective structural and fluid subproblem. This also permits to take into accountvarious nonlinear effects in either domain in a straight-forward way. Even if the aerodynamic forces remain linearin the vibration amplitude, the fully-coupled solver might be an attractive alternative. The reason for this is that the18nfluence coefficients need to be computed for a large number of component modes in order to accurately describevariation of the vibrational deflection shape within the refined energy method.As benchmark, a state-of-the-art aeroelastic model of a bladed disk subjected to nonlinear contact interactions at theinterlocked tip shrouds was considered. The fully-coupled solver showed excellent convergence behavior for reasonablecomputational effort. The good agreement between fully-coupled solver and refined energy method confirms that theaerodynamic forces can, in a wide range, be well approximated as amplitude-linear and using a linear interpolation ofthe modal aerodynamic influence coefficients with respect to frequency. In the first configuration, a stable limit cycleunder relatively low aerodynamic self-excitation was analyzed. The oscillation frequency and vibrational deflectionshape were thus close to the linear case under sticking contact conditions here. Hence, the conventional energy methodyields a reasonable approximation and is therefore in good agreement with the newly developed methods. The secondconfiguration demonstrated, for the first time, the phenomenon of nonlinear blade flutter instability, i.e. , a situationwhere the static equilibrium (sticking contact conditions) is stable, but the vibrations start to grow unboundedly oncethe initial vibration amplitude exceeds a certain threshold. Neither linear theory nor the conventional energy methodare able to predict this behavior, while the newly developed methods are both able to determine the critical amplitude(stability limit).An open question is under which circumstances a considerable amplitude-nonlinearity of the aerodynamic forces isto be expected and how the fully-coupled solver performs in this case. It would also be interesting to assess theperformance of the methods for a case of more severe aerodynamic influence, where the aerodynamic stiffness is alsoexpected to affect the oscillation frequency and perhaps the vibrational deflection shape. Finally, the full-coupledsolver could be improved with regard to its initialization, as this currently involves the relatively costly computationof the aerodynamic influence coefficient matrix (with respect to the generalized coordinates of the reduced structuralmodel).
Acknowledgements
The research project was financially supported by the German Research Foundation (DFG) (funding no. 382141955)and the FVV (Research Association for Combustion Engines eV) (FVV funding no. 6013081). The research projectwas performed by the Institute of Propulsion Technology (IAT) at the German Aerospace Center under the directionof Prof. Reinhard Mönig and by the Institute of Aircraft Propulsion Systems (ILA) at University of Stuttgart under thedirection of Jun.-Prof. Malte Krack. The project was conducted by an expert group under the direction of Dr. AndreasHartung (MTU Aero Engines AG). The authors greatfully acknowledge the support received from the FVV (ResearchAssociation for Combustion Engines eV) and the DFG (German Research Foundation).19 ppendix A. Convergence Details
Figure A.1: Convergence of the dominant amplitude and frequency of the limit cycle for the sbtability limit with full-influence-matrix linearization.
References [1] A. Srinivasan, Flutter and resonant vibration characteristics of engine blades, Journal of Engineering for Gas Turbines and Power 119 (1997)– (Oct. 1997).[2] S. Leichtfuss, F. Holzinger, C. Brandstetter, F. Wartzek, H. P. Schiffer, Aeroelastic investigation of a transonic research compressor, in:Proceedings of ASME Turbo Expo 2013: Turbine Technical Conference and Exposition, no. 55270, 2013, pp. V07BT33A006– (2013). doi:10.1115/GT2013-94730 .URL http://dx.doi.org/10.1115/GT2013-94730 [3] M. May, Sensitivity analysis with respect to flutter-free design of compressor blades, in: ASME Turbo Expo, no. GT2010-23557, Glasgow,UK, 2010 (Jun. 2010).[4] H. Schönenborn, T. Breuer, Aeroelasticity at reversed flow conditions—part ii: application to compressor surge, Journal of turbomachinery134 (6) (2012) 061031 (2012).[5] J. J. Waite, R. E. Kielb, Physical understanding and sensitivities of low pressure turbine flutter, Journal of Engineering for Gas Turbines andPower 137 (1) (2014) 012502–012502 (Aug. 2014). doi:10.1115/1.4028207 .URL http://dx.doi.org/10.1115/1.4028207 [6] J. J. Waite, R. E. Kielb, Shock structure, mode shape, and geometric considerations for low-pressure turbine flutter suppression, in: TurboExpo: Power for Land, Sea, and Air, Vol. 49842, American Society of Mechanical Engineers, 2016, p. V07BT34A009 (2016).[7] M. Lassalle, C.M. Firrone, A parametric study of limit cycle oscillation of a bladed disk caused by flutter and friction at the blade root joints,Journal of Fluids and Structures 76 (Supplement C) (2018) 349–366 (2018). doi:10.1016/j.jfluidstructs.2017.10.004 .[8] C. Martel, R. Corral, R. Ivaturi, Flutter amplitude saturation by nonlinear friction forces: Reduced model verification, Journal of Turboma-chinery 137 (4) (2014) 041004 (2014). doi:10.1115/1.4028443 .[9] J. Gross, M. Krack, Multi-wave vibration caused by flutter instability and nonlinear tip shroud friction, Journal of Engineering for Gas Turbinesand Power (2019).[10] M. Krack, L. Salles, F. Thouverez, Vibration prediction of bladed disks coupled by friction joints, Archives of Computational Methods inEngineering 24 (3) (2017) 589–636 (2017).[11] M. Nowinski, J. Panovsky, Flutter mechanisms in low pressure turbine blades, Journal of Engineering for Gas Turbines and Power 122 (1)(1999) 82–88 (1999). doi:10.1115/1.483179 .[12] C. E. Seeley, C. Wakelam, X. Zhang, D. Hofer, W.-M. Ren, Investigations of flutter and aero damping of a turbine blade: Part 1 - experimentalcharacterization, in: ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, no. 49842, 2016, pp. V07BT34A024–(2016). doi:10.1115/GT2016-57930 .URL http://dx.doi.org/10.1115/GT2016-57930 [13] R. Corral, J. M. Gallardo, Verification of the vibration amplitude prediction of self-excited lpt rotor blades using a fully coupled time-domainnon-linear method and experimental data, in: Turbo Expo: Power for Land, Sea, and Air, Vol. 43154, 2008, pp. 835–847 (2008).[14] W.-M. Ren, C. E. Seeley, X. Zhang, B. E. Mitchell, H. Ju, Investigations of flutter and aero damping of a turbine blade: Part 2 - numericalsimulations, in: ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, 2016, pp. V07BT34A025–V07BT34A025(2016). doi:10.1115/GT2016-57935 .URL http://dx.doi.org/10.1115/GT2016-57935 [15] T. R. Müller, D. M. Vogt, K. Vogel, B. A. Phillipsen, Influence of intrarow interaction on the aerodynamic damping of an axial turbine stage,in: Turbo Expo: Power for Land, Sea, and Air, 2018 (2018). doi:10.1115/GT2018-76777 .[16] A. Sinha, J. H. Griffin, Friction damping of flutter in gas turbine engine airfoils, AIAA Journal of Aircraft 20 (4) (1983) 372–376 (1983).[17] A. Sinha, J. H. Griffin, Effects of friction dampers on aerodynamically unstable rotor stages, AIAA Journal 23 (2) (1985) 262–270 (1985). doi:10.2514/3.8904 .[18] A. Sinha, J. H. Griffin, Stability of limit cycles in frictionally damped and aerodynamically unstable rotor stages, Journal of Sound andVibration 103 (3) (1985) 341–356 (1985). doi:10.1016/0022-460X(85)90427-4 .
19] A. Hartung, J. Gross, T. Heinze, L. Panning-von Scheidt, M. Krack, H.-P. Hackenberg, Rig and engine validation of the non-linear forcedresponse analysis performed by the tool oragl, J. Eng. Gas Turbines Power, accepted (2018).[20] R. A. Leyes, W. A. Fleming, The history of North American small gas turbine aircraft engines, Aiaa, 1999 (1999).[21] X. Wu, M. Vahdati, C. Schipani, M. Imregun, Analysis of low-pressure turbine flutter for different shroud interfaces, in: ASME Turbo Expo2007: Power for Land, Sea, and Air, American Society of Mechanical Engineers, 2007, pp. 629–636 (2007).[22] E. P. Petrov, Analysis of flutter-induced limit cycle oscillations in gas-turbine structures with friction, gap, and other nonlinear contactinterfaces, Journal of Turbomachinery 134 (6) (2012) 061018/1—061018/13 (2012). doi:10.1115/1.4006292 .[23] C. Fuhrer, D. M. Vogt, On the impact of simulation approaches on the predicted aerodynamic damping of a low pressure steam turbine rotor,in: ASME TurboExpo 2017, 2017 (2017). doi:10.1115/GT2017-63401 .URL http://dx.doi.org/10.1115/GT2017-63401 [24] C. Martel, R. Corral, J. M. Llorens, Stability increase of aerodynamically unstable rotors using intentional mistuning, Journal of Turboma-chinery 130 (1) (2007). doi:10.1115/1.2720503 .[25] K. D’Souza, C. Jung, Epureanu B. I., Analyzing mistuned multi-stage turbomachinery rotors with aerodynamic effects, Journal of Fluids andStructures 42 (0) (2013) 388–400 (2013). doi:10.1016/j.jfluidstructs.2013.07.007 .[26] C. Berthold, C. Frey, G. Dhondt, H. Schönenborn, Fully coupled aeroelastic simulations of limit cycle oscillations in the time domain, in:Proceedings of the 15th ISUAAAT, 2018 (September 2018).URL https://elib.dlr.de/127137/ [27] H.-S. Im, G.-C. Zha, Simulation of non-synchronous blade vibration of an axial compressor using a fully coupled fluid/structure interaction,no. 44731, 2012, pp. 1395–1407 (2012). doi:10.1115/GT2012-68150 .URL http://dx.doi.org/10.1115/GT2012-68150 [28] R. Corral, J. M. Gallardo, Nonlinear dynamics of bladed disks with multiple unstable modes: Aiaa journal, AIAA Journal 52 (6) (2014)1124–1132 (2014). doi:10.2514/1.J051812 .[29] H.-P. Kersken, C. Frey, G. Ashcroft, H. Schönenborn, Flutter analysis of an embedded blade row with a harmonic balance solver, in:Proceedings of 12th European Conference on Turbomachinery Fluid dynamics & Thermodynamics, 2017 (2017).[30] M. Krack, J. Gross, The harmonic balance method and its application to nonlinear vibrations: Introduction and current state of the art,submitted to Mechanical Systems and Signal Processing (2018).[31] K. Ekici, H. Huang, An assessment of frequency-domain and time-domain techniques for turbomachinery aeromechanics, in: 30th AIAAApplied Aerodynamics Conference 2012, 2012, pp. 1807–1823 (Jun. 2012).URL [32] D. Su, W. Zhang, Z. Ye, A reduced order model for uncoupled and coupled cascade flutter analysis, Journal of Fluids and Structures 61 (2016)410–430 (2016).[33] A. Cadel, G. N. Boum, F. Thouverez, A. Dugeai, M.-O. Parent (Eds.), Computing fluid structure interaction coupling time spectral method(TSM) and harmonic balance method (HBM), Vol. 7B-2017, 2017. doi:10.1115/GT2017-64260 .[34] Y. Gong, W. Zhang, Efficient aeroelastic solution based on time-spectral fluid–structure interaction method, AIAA Journal (2019) 3014–3025(2019).[35] C. R. Berthold, Development of a coupled fluid-structure simulation method in the frequency domain, Master’s thesis, Faculty of Mechanical,Maritime and Materials Engineering (3mE), Delft University of Technology, p&E report number: 2760 (2016).[36] C. Berthold, C. Frey, H. Schönenborn, Coupled fluid structure simulation method in the frequency domain for turbomachinery applications,in: ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, American Society of Mechanical Engineers, 2018, pp.V07CT36A014–V07CT36A014 (2018).[37] S. T. Wei, C. Pierre, Localization phenomena in mistuned assemblies with cyclic symmetry. part i: Free vibrations, Journal of Vibration,Acoustics, Stress, and Reliability in Design 110 (4) (1988) 429–438 (1988). doi:10.1115/1.3269547 .[38] R. R. Craig, M. C. Bampton, Coupling of substructures for dynamic analysis, AIAA Journal 6 (7) (1968) 1313–1319 (1968).[39] M. Krack, J. Gross, Harmonic Balance for Nonlinear Vibration Problems, Springer, 2019 (2019). doi:10.1007/978-3-030-14023-6 .[40] D. C. Wilcox, Reassessment of the scale-determining equation for advanced turbulence models, AIAA J. 26 (11) (1988) 1299–1310 (November1988).[41] G. Ashcroft, C. Frey, H.-P. Kersken (Eds.), On the development of a harmonic balance method for aeroelastic analysis, 2014.[42] D. Schlüß, C. Frey, Time domain flutter simulations of a steam turbine stage using sptectral 2d non-reflecting boundary conditions, in: 15thInternational Symposium on Unsteady Aerodynamics Aeroacoustics and Aeroelasticity Turbomachines, 2018 (06 2018).[43] C. Voigt, C. Frey, H.-P. Kersken, Development of a generic surface mapping algorithm for fluid-structure-interaction simulations in turbo-machinery, in: J. C. F. Pereira, A. Sequeira, J. M. C. Pereira (Eds.), V European Conference on Computational Fluid Dynamics ECCOMASCFD 2010, 2010 (Jun. 2010).URL http://elib.dlr.de/64893/ [44] D. Ottl, Modelle und mathematische beschreibung, versuchstechniken zur messung von dämpfungskenngrößen, in: VDI-Richtlinie 3830:Werkstoff- und Bauteildämpfung, 2007, pp. 1–16 (2007).[45] H.-P. Kersken, C. Frey, C. Voigt, G. Ashcroft, Time-linearized and time-accurate 3D RANS methods for aeroelastic analysis in turbomachinery,J. Turbomach. 134 (5) (2012) 051024 (2012). doi:10.1115/1.4004749 .[46] G. Kahl, Aeroelastic effects of mistuning and coupling in turbomachinery bladings, Ph.D. thesis, École polytechnique fédérale de Lausanne(2002). doi:10.5075/epfl-thesis-2629 .[47] M. Krack, Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept, Computers and Structures 154(2015) 59–71 (2015). doi:10.1016/j.compstruc.2015.03.008 .[48] C. Frey, G. Ashcroft, H.-P. Kersken, C. Voigt, A harmonic balance technique for multistage turbomachinery applications, in: ASME TurboExpo 2014: Turbine Technical Conference and Exposition, no. 45615, 2014, p. V02BT39A005 (2014). doi:10.1115/gt2014-25230 .URL http://dx.doi.org/10.1115/GT2014-25230
49] G. Ashcroft, C. Frey, H.-P. Kersken, On the development of a harmonic balance method for aeroelastic analysis, in: 6th European Conferenceon Computational Fluid Dynamics (ECFD VI), 2014 (2014).49] G. Ashcroft, C. Frey, H.-P. Kersken, On the development of a harmonic balance method for aeroelastic analysis, in: 6th European Conferenceon Computational Fluid Dynamics (ECFD VI), 2014 (2014).