Cassoni, G.

Cocco, A.

Masarati, P.

Tamer, A.

Zanoni, A.

English

Archivio istituzionale della ricerca - Politecnico di Milano

Paper 124Tiltrotor Whirl-Flutter Stability Investigation using LyapunovCharacteristic Exponents and Multibody DynamicsGianni Cassoni†, Alessandro Cocco†,∗, Aykut Tamer?, Andrea Zanoni†, Pierangelo Masarati††Politecnico di Milano, Milano - Italy?Imperial College London, London - UK* alessandro.cocco@polimi.itAbstractThis work discusses the use of Lyapunov Characteristic Exponents to assess the stability of nonlinear, time-dependentmechanical systems. Specific attention is dedicated to methods capable of estimating the largest exponent without requir-ing the Jacobian matrix of the problem, which can be applied to time histories resulting from existing multibody solvers.Tiltrotor whirl-flutter stability is analyzed. With respect to the available literature, the proposed method does not require thesystem to be strictly periodic, no linearization is required about a reference steady solution, and characteristic nonlinearaspects of stationary solutions like limit cycle oscillations are correctly identified and pointed out. A limitation lies in theability to correctly identify the stability but no information is inferred about the related characteristic frequencies/periods, ifany.IntroductionStability assessment is a fundamental aspect of the anal-ysis and design of dynamical systems. The foundationsof modern stability theory lie in Aleksandr M. Lyapunov’swork [1].Although for nonlinear problems of the general formẋ = f(x, t)(1)stability is a local property of a specific solution, x(t), re-sulting from a specific set of initial conditions, x(t0) = x0,for Linear, Time-Invariant (LTI) problems, like those result-ing from linearization about a steady reference condition inthe formẋ = Ax,(2)it becomes a characteristic of the entire system. It can beevaluated from the inspection of the real part of the eigen-values of matrix A.In many applications associated with rotorcraft dynam-ics, and specifically in helicopter rotor aeromechanics, prob-lems often need to be formulated as time-periodic, althoughoften still linear, through linearization about a periodic or-bit, resulting in Linear, Time-Periodic (LTP) problems. Inthis case, stability assessment may exploit the periodic-ity of the motion through the well-known Floquet-Lyapunovapproach. Asymptotic stability is associated with contrac-tion, over a period T , of the so-called “monodromy matrix”,which corresponds to the State Transition Matrix (STM) ofthe problem over one period. This approach is relativelypopular in the rotorcraft aeromechanics community [2–5].In general, however, when problems are non-linear andsubjected to non-(strictly) periodic time dependence, asmay occur in many transient-related problems, the abilityto evaluate the stability of reference trajectories can be ex-tremely useful. The use of Lyapunov Characteristic Expo-nents (LCE) or, in short, Lyapunov Exponents (LE), hasbeen recently proposed in the field of rotorcraft aerome-chanics [6–8], also with a focus on their sensitivity study[9–12], and other aerospace-related applications [10].Lyapunov Characteristic ExponentsThe LCEs indicate the rate of expansion or contraction ofperturbations of a generic solution of the nonlinear differen-tial problem of Eq. (1) along independent directions in thestate space. As such, they describe the stability of the ref-erence solution with respect to such directions.Consider a solution x(t) of Eq. (1) for t ≥ t0 (some au-thors refer to it as the ‘fiducial trajectory’), and a solutionix(t) of the problemiẋ = f/x∣∣x(t),t ix, ix(t0) = ix0(3)for a perturbation ix0 of arbitrary magnitude and direction.LCEs are defined asλi = limt→∞1tlog‖ix(t)‖ .(4)Each λi is calculated from one of n linearly independentix0, the equivalent of the principal directions of a LTI prob-lem. Since Eq. (3) is linear time-dependent, its solutionPresented at the 48th European Rotorcraft Forum, Winterthur, Switzerland, 6-8 September, 2022This work is licensed under the Creative Commons Attribution International License (CC BY). Copyright © 2022 by author(s).Page 1 of 6can be expressed through the State Transition Matrix (STM)Y(t, t0),ix = Y(t, t0) ix0(5)The STM Y(t, t0) is the solution of Ẏ = f/x∣∣x(t),t Y, withY(t0, t0) = I, namely the matrix that represents the evolu-tion of the problem’s state from time t0 to t. Then Eq. (4) isequivalent toλi = limt→∞1tRe(log(eig(Y(t, t0)))) .(6)When all LCEs are negative, the solution is exponentiallystable. When at least one LCE is positive, the solution isunstable, or leads to a chaotic attractor. When the largestLCE is zero, or the largest LCEs are zero, a limit cycle os-cillation (LCO) is expected; i.e., there exists one direction,or multiple independent directions, in the state space alongwhich the solution neither expands nor contracts. In case ofmultiple largest LCEs equal to zero, a higher order periodicor quasi-periodic attractor exists, e.g. a torus.Note the analogy with the LTI case, since Y(t, t0)LTI≡eA(t−t0) and thus λiLTI= Re(eig(A)), and with the LTP one, inwhichλiLTP=1TRe(log(eig(Y(t0 +T, t0)))) .(7)In this sense, one may consider the LCEs as sort of theeigenvalues of matrix f/x, averaged over time.LCEs are often called the ‘spectrum’ of the associatedproblem, much like the spectrum of LTI problems is repre-sented by the eigenvalues of matrix A = f/x.In most cases, the definition of Eq. (6) cannot be used inpractice, because usually at least some of the elements ofthe STM contract to zero, in case of asymptotic stability ofthe solution, or expand to infinity, in case of instability, result-ing in either under- or overflowing, or a mix of both. As pro-posed in the work of Benettin et al. [13], to estimate LCEsin practice one needs to exploit the re-orthogonalization ofthe local directions of evolution of the solution. Alternativeapproaches based on well known orthogonal decomposi-tions (the Singular Value Decomposition (SVD) and the QRfactorization, respectively) have been proposed [13–15].Jacobian-Less Methods: Max LCEfrom Time SeriesThe estimation of LCEs from problems of interest in the ro-torcraft field using those methods has been recently dis-cussed in [6,8]. However, those methods require the knowl-edge of the Jacobian matrix of the problem, and the abilityto integrate the associated linear, time dependent problem.The largest LCE, or Maximum LCE (MLCE) can also beestimated directly from a time series. This makes the ap-proach rather interesting because it can be used in conjunc-tion with experimental data, or with the results of numericalsimulations performed using complex multibody models. Al-though the full spectrum cannot be obtained, estimation ofthe MLCE allows to identify the LCE that belongs to theleast damped principal direction of the system, which is themost critical stability indicator. Among the algorithms pro-posed in the literature (see for example [16]), in this workthe one proposed by Rosenstein et al. [17], described be-low, is used.The trajectory matrix, X, is constructed from an N-pointtime series xi, i = 1, ...,N, using the time delay method.Each row of matrix X is a phase-space vector, namelyX =[X1 X2 . . . Xm](8)withXk =[x1+(k−1)J x2+(k−1)J . . . xM+(k−1)J]T(9)with k = 1, ...,m. Thus, X ∈ RM×m, with m, M, J, and Nrelated byM = N− (m−1)J(10)where m is the embedding dimension, N the length of thetime series, and J represents the so-called reconstructiondelay. The embedding dimension is usually estimated inaccordance with Takens’ theorem, i.e., m > 2n.After constructing the trajectory matrix, the algorithm lo-cates the nearest neighbor, X ĵ, of each point on the trajec-tory. It is found by searching for the point that minimizesthe distance from each particular reference point, X j. Thedistance is expressed asd j(0) = minX ĵ∥∥∥X j−X ĵ∥∥∥(11)where d j(0) is the initial distance from the jth point to itsnearest neighbor, and ‖·‖ denotes the Euclidean norm.An additional constraint is that nearest neighbors havea temporal separation greater than the mean period (T̄ , thereciprocal of the mean frequency of the power spectrum,although it can be expected that any comparable estimate,e.g., using the median frequency of the magnitude spec-trum, yields equivalent results) of the time series,∣∣ j− ĵ∣∣> T̄(12)Thanks to this, each pair of neighbors can be considered asnearby initial conditions for different trajectories.The largest Lyapunov exponent is then estimated as themean rate of separation of the nearest neighbors. The jthpair of nearest neighbors diverge approximately at a rategiven by the largest Lyapunov exponent:d j(i)≈C jeλ1(i∆t)(13)where C j is the initial separation. By taking the logarithm ofboth sideslogd j(i)≈ logC j +λ1(i∆t)(14)Page 2 of 6which represents a set of approximately parallel lines, forj = 1, ...,M, each with a slope roughly proportional to λ1.The largest Lyapunov exponent is calculated using a least-squares fit to the “average” line defined byy(i) =1∆t〈logd j(i)〉(15)where 〈·〉 denotes the average over all values of j.Numerical ResultsThe use of the Jacobian-less method by Rosenstein etal. was proposed in [18] to address the rather importantrotorcraft-related problem of ground resonance with rela-tively simple models, for validation purposes, including amodel formulated using MBDyn1, a free general-purposemultibody solver [19]. In such case, this method is of par-ticular interest, since it does not require the capability of thesolver to expose the Jacobian matrix of the problem, andcan thus be applied to results obtained from any solver.The problem analyzed in the following is the aeroelasticsimulation of a tiltrotor semi-span model (Fig. 1) in prox-imity of whirl-flutter conditions, using a rather sophisticatedaeroservoelastic model [20].Figure 1: Tiltrotor semispan model.The entire model can be divided into two main sub-components: the wing-pylon assembly and the propro-tor. The wing-pylon model has been developed using ageometrically exact, composite ready beam finite elementmodel [21, 22] to reproduce the fundamental frequenciesand mode shapes of an equivalent finite element stickmodel, which was tuned to match the full-scale aircraft dy-namics at the rotor. The proprotor is a three-bladed stiff-in-plane rotor with gimballed hub (Fig. 2). It consists of thepitch control chain, three blades, and the yoke, all modeledusing the previously mentioned beam elements.Figure 2: Proprotor multibody model.In the original formulation, the algorithm proposed byRosenstein et al. is applied to a single time series, namely ascalar signal. In complex multivariable problems, the choiceof a suitable signal can be critical, along with some param-eters of the method, related to the duration of the signaland other factors discussed in the original paper. To over-come this critical aspect, a Proper Orthogonal Decompo-sition of the multiple signals resulting from the multibodydynamics simulation is performed, and the transformed sig-nals corresponding to the largest singular values are con-sidered. Figure 3 compares the results obtained with theproposed analysis to that resulting from the Periodic Oper-ational Modal Analysis (POMA) proposed in [23–25]. TheMLCEs estimated from the first three principal components(PC) are considered, starting from the top.The first PC (top plot) shows results compatible withthose of the first wing chordwise mode. This occurs sincethe excitation that caused the transient was dominated by acollective pitch perturbation, which directly excites the wingchordwise mode.The second PC (mid plot) shows results compatiblewith those of the first wing bending mode, the most lightlydamped one, except for the analysis at the highest speed,which is dominated by whirl-flutter, with a substantial partic-ipation of wing torsion.1https://www.mbdyn.org/Page 3 of 6a) MLCE from first Principal Component (mainly first wing chordwise mode)b) MLCE from second Principal Component (mainly first wing bending mode)c) MLCE from third Principal Component (mainly first wing torsion mode)Figure 3: MLCE compared with real part of eigenvalues from POMA.Page 4 of 6Indeed, the wing torsional mode dominates the third PC(bottom plot), which shows greater damping that the othermodes in the speed range below whirl-flutter.ConclusionsThis paper discussed the use of estimated Lyapunov Char-acteristic Exponents to evaluate the stability of a complexmechanical system, namely a tiltrotor semi-span modelclose to whirl-flutter conditions. The problem was formu-lated without a priori linearization or coordinate transfor-mation to reduce the azimuthal dependence of the equa-tions. Nonlinear phenomena, like limit-cycle oscillations,can be detected. By resorting to Jacobian-less methods,the Maximum LCE, i.e., the most critical one regarding sta-bility, can be estimated from time series computed usinggeneral-purpose multibody formulations, without the needto modify existing solvers.Copyright StatementThe authors confirm that they, and/or their company or or-ganization, hold copyright on all of the original material in-cluded in this paper. The authors also confirm that theyhave obtained permission, from the copyright holder of anythird party material included in this paper, to publish it aspart of their paper. The authors confirm that they give per-mission, or have obtained permission from the copyrightholder of this paper, for the publication and distribution ofthis paper as part of the ERF proceedings or as individualoffprints from the proceedings and for inclusion in a freelyaccessible web-based repository.References[1] Aleksandr Mikhailovich Lyapunov. The General Problem ofthe Stability of Motion. Taylor & Francis, London - Washing-ton, DC, 1992. Translated and edited by A. T. Fuller.[2] David A. Peters and Kurt H. Hohenemser. Application ofthe Floquet transition matrix to problems of lifting rotor stabil-ity. Journal of the American Helicopter Society, 16(2):25–33,1971. doi:10.4050/JAHS.16.25.[3] James C. Biggers. Some approximations to the flapping sta-bility of helicopter rotors. Journal of the American HelicopterSociety, 19(4):24–33, 1974. doi:10.4050/JAHS.19.24.[4] P. Friedmann and J. Silverthorn. Aeroelastic stability ofperiodic systems with application to rotor blade flutter. InAIAA/ASME/SAE 15th Structures, Structural Dynamics andMaterials Conference, Las Vegas, Nevada, April 17–19 1974.[5] C. Eugene Hammond. An application of Floquet the-ory to prediction of mechanical instability. Journalof the American Helicopter Society, 19(4):14–23, 1974.doi:10.4050/JAHS.19.14.[6] Aykut Tamer and Pierangelo Masarati. Helicopter rotoraeroelastic stability evaluation using Lyapunov exponents.In 40th European Rotorcraft Forum, Southampton, UK,September 2–5 2014.[7] Aykut Tamer and Pierangelo Masarati. Do we really needto study rotorcraft as linear periodic systems? In AHS 71stAnnual Forum, Virginia Beach, VA, USA, May 5–7 2015.[8] Aykut Tamer and Pierangelo Masarati. Stability of nonlinear,time-dependent rotorcraft systems using Lyapunov charac-teristic exponents. Journal of the American Helicopter Soci-ety, 61(2), April 2016. doi:10.4050/JAHS.61.022003.[9] Pierangelo Masarati and Aykut Tamer. Sensitivity of trajec-tory stability estimated by Lyapunov characteristic exponents.Aerospace Science and Technology, 47:501–510, November2015. doi:10.1016/j.ast.2015.10.015.[10] Aykut Tamer and Pierangelo Masarati. Aeroelastic stabil-ity estimation of control surfaces with freeplay nonlinear-ity. l’Aerotecnica Missili e Spazio, 96(3):154–164, 2017.doi:10.1007/BF03404750.[11] Aykut Tamer and Pierangelo Masarati. Sensitivity ofLyapunov exponents in design optimization of nonlineardampers. J. of Computational and Nonlinear Dynamics,14(2):021002, February 2019. doi:10.1115/1.4041827.[12] Aykut Tamer and Pierangelo Masarati. Generalized quanti-tative stability analysis of time-dependent comprehensive ro-torcraft systems. Aerospace, 9(1), 2022.[13] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, andJean-Marie Strelcyn. Lyapunov characteristic exponents forsmooth dynamical systems and for Hamiltonian systems; amethod for computing all of them. part 1: Theory. Mecca-nica, 15(1):9–20, March 1980. doi:10.1007/BF02128236.[14] Karlheinz Geist, Ulrich Parlitz, and Werner Lauterborn. Com-parison of different methods for computing Lyapunov expo-nents. Progress of Theoretical Physics, 83(5):875–893, May1990. doi:10.1143/PTP.83.875.[15] Luca Dieci and Erik S. Van Vleck. Lyapunov spec-tral intervals: Theory and computation. SIAM Jour-nal on Numerical Analysis, 40(2):516–542, 2002.doi:10.1137/S0036142901392304.[16] Ch. Skokos. The Lyapunov Characteristic Exponents andTheir Computation, pages 63–135. Springer Berlin Hei-delberg, Berlin, Heidelberg, 2010. doi:10.1007/978-3-642-04458-8 2.[17] Michael T. Rosenstein, James J. Collins, and Carlo J. DeLuca. A practical method for calculating largest Lyapunovexponents from small data sets. Physica D: Nonlinear Phe-nomena, 65(1–2):117–134, May 1993. doi:10.1016/0167-2789(93)90009-P.[18] Gianni Cassoni, Andrea Zanoni, Aykut Tamer, andP. Masarati. Stability of rotorcraft ground resonance by es-timating Lyapunov characteristic exponents from multibodydynamics. In ASME IDETC/CIE 2022, St. Louis, Missouri,USA, August 14–17 2022. doi:10.1115/DETC2022-88995.[19] Pierangelo Masarati, Marco Morandini, and Paolo Man-tegazza. An efficient formulation for general-purpose multi-body/multiphysics analysis. J. of Computational and Nonlin-ear Dynamics, 9(4):041001, 2014. doi:10.1115/1.4025628.[20] Alessandro Cocco, Alberto Savino, and Pierangelo Masarati.Flexible multibody model of a complete tiltrotor for aeroser-voelastic analysis. In ASME IDETC/CIE 2022, St. Louis,Missouri, USA, August 14–17 2022. doi:10.1115/DETC2022-88973.Page 5 of 6[21] Gian Luca Ghiringhelli, Pierangelo Masarati, and PaoloMantegazza. A multi-body implementation of finite volumeC0 beams. AIAA Journal, 38(1):131–138, January 2000.doi:10.2514/2.933.[22] Olivier A. Bauchau, Peter Betsch, Alberto Cardona, Jo-hannes Gerstmayr, Ben Jonker, Pierangelo Masarati, andValentin Sonneville. Validation of flexible multibody dynam-ics beam formulations using benchmark problems. MultibodySystem Dynamics, 37(1):29–48, 2016. doi:10.1007/s11044-016-9514-y.[23] Alessandro Cocco, Pierangelo Masarati, Stefan van ’t Hoff,and Bart Timmerman. Tiltrotor whirl flutter analysis in supportof NGCTR aeroelastic wind tunnel model design. In 47th Eu-ropean Rotorcraft Forum, Glasgow, UK (virtual), September7–10 2021.[24] Alessandro Cocco, Alberto Savino, Andrea Zanoni, andPierangelo Masarati. Comprehensive simulation of a com-plete tiltrotor with pilot-in-the-loop for whirl-flutter stabilityanalysis. In 48th European Rotorcraft Forum, Winterthur,Switzerland, September 6–8 2022.[25] Alessandro Cocco, Pierangelo Masarati, Stefan van ’tHoff, and Bart Timmerman. Numerical whirl-flutter anal-ysis of a tiltrotor semi-span wind tunnel model. CEASAeronautical Journal, available online September 2022.doi:10.1007/s13272-022-00605-2.Page 6 of 6

Tiltrotor Whirl-Flutter Stability Investigation Using Lyapunov Characteristic Exponents and Multibody Dynamics

https://re.public.polimi.it/bitstream/11311/1220940/1/CASSG02-22.pdf

Tiltrotor Whirl-Flutter Stability Investigation Using Lyapunov Characteristic Exponents and Multibody Dynamics

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Politecnico di Milano