Generalized Quantitative Stability Analysis of Time-Dependent Comprehensive Rotorcraft Systems

Abstract

Rotorcraft stability is an inherently multidisciplinary area, including aerodynamics of rotor and fuselage, structural dynamics of flexible structures, actuator dynamics, control, and stability augmentation systems. The related engineering models can be formulated with increasing complexity due to the asymmetric nature of rotorcraft and the airflow on the rotors in forward flight conditions. As a result, linear time-invariant (LTI) models are drastic simplifications of the real problem, which can significantly affect the evaluation of the stability. This usually reveals itself in form of periodic governing equations and is solved using Floquet’s method. However, in more general cases, the resulting models could be non-periodic, as well, which requires a more versatile approach. Lyapunov Characteristic Exponents (LCEs), as a quantitative method, can represent a solution to this problem. LCEs generalize the stability solutions of the linear models, i.e., eigenvalues of LTI systems and Floquet multipliers of linear time-periodic (LTP) systems, to the case of non-linear, time-dependent systems. Motivated by the need for a generic tool for rotorcraft stability analysis, this work investigates the use of LCEs and their sensitivity in the stability analysis of time-dependent, comprehensive rotorcraft models. The stability of a rotorcraft modeled using mid-fidelity tools is considered to illustrate the equivalence of LCEs and Floquet’s characteristic coefficients for linear time-periodic problems.

1. Introduction

Rotorcraft use rotating blades to generate the loads required for flight, i.e., lift, propulsion, and control, as needed to take-off and land vertically, as well as hover [1]. They can fulfill many roles that fixed-wing aircraft cannot perform as effectively, such as search and rescue, if they can at all. To provide a sufficient level of thrust at hover and low airspeeds, a relative motion must be maintained between the rotating wings and the air. Rotors should have a long enough span to enable the efficient lifting, propulsion, and control of the vehicle [2]. These three fundamental functions significantly affect each other, and lead to cross-couplings, which result in more complex flight physics compared to fixed-wing aircraft [3].

Stability assessment is an important task in helicopter design, especially considering that from an aeromechanical viewpoint they are inherently unstable [4]. A specific form of complexity arises from the rotation of the blades in forward flight, which induces periodic aerodynamic loading on the rotor blades. Additionally, blades in the advancing azimuthal range add their rotational velocity with a component of the forward flight speed of the vehicle, and exhibit substantial compressibility effects at high forward speeds, with the blade tip approaching Mach 1. Conversely, blades in the other half of the azimuth move in a direction opposite to that of the vehicle, resulting in low relative speed, the need for high angles of attack to compensate for the reduction of lift, proximity to stalled flow, and even reverse flow conditions near the blade root. As a consequence of extended stall and reverse flow regions, high levels of drag result, with a significant drop in the lift. In addition to the time-dependence of aerodynamic origin, similar effects also occur if an unbalance is present in the rotating elements. The latter can usually happen as a result of deviations from nominal mass distribution, or aerodynamic shape distortion [5]. Another source could be a malfunction of one of the repetitive mechanical elements, such as lead-lag dampers [6]. Whatever the reason, the loss of symmetry in the rotation field results in

• alternating loads that excite the fuselage, and
• alternating system properties leading to a time-dependent system.

While the former represents a forced vibration problem, the latter results in a time-variant problem, whose stability analysis is more difficult.

In dynamical systems, the study of nearby solutions of an equilibrium point is referred to as stability [7]. In other words, the system response after a perturbation of the system’s state is introduced about an equilibrium condition [8]. There are several methods to assess the stability of a system. A conventional one is perturbing the system through typical input or disturbance sources and observing the free-decay characteristic. This could be performed by conducting experiments on a manufactured product or through computer simulations. Alternatively, spectral methods exist, which quantify the decaying behaviour of independent components of the solution, which is referred to as the system’s exponential multipliers or characteristic exponents. Once the stability spectrum that quantifies stability is obtained, it is straightforward to interpret the results and track the changes in stability characteristics with the evolution of the design. Besides, the spectral methods are also suitable for the analytical estimation of the sensitivity of stability estimates. This work focuses on spectral methods of estimating stability.

The validity of rotorcraft stability estimation is defined by the capability of mathematical and numerical tools used in evaluating the stability and the simplifications made during the process. For the linear, time-invariant (LTI) case, stability can be assessed using eigenanalysis of the system’s matrix. A more difficult case is that of linear, time-periodic (LTP) systems, a special class of linear, time-variant (LTV) problems, which can be considered to address the intrinsic periodicity of idealized rotorcraft, especially originating from the main rotor’s aerodynamics in non-axial flow conditions, or dissimilarities in the rotor components, as in ground resonance with degraded dampers. A widespread tool, in this case, is Floquet’s theory, which represents a natural extension of eigenvalue analysis for LTP systems [9]. The eigenvalue decomposition for LTI systems and Floquet analysis of LTP problems may successfully evaluate stability in several flight conditions; however, rotorcraft dynamics is generally described by non-linear, generically time-dependent equations [10]. Therefore, these methods imply simplifications when the system is neither time-invariant nor periodic. For these reasons, capturing the effects of time-dependence in rotorcraft stability could provide more insight into increasingly complex rotorcraft stability problems.

To achieve a generalized framework for the stability analysis of rotorcraft, a technique that can estimate the stability of time-dependent systems is required. This method is further needed to require no simplifications regarding time-dependence nor any assumption of periodicity. Additionally, the technique should give the same results and interpretation when applied to relatively simple models of time-invariant and periodic systems. Lyapunov Characteristic Exponents (LCE) satisfy all these requirements and can represent a universal framework for rotorcraft aeroelastic stability. LCEs are the indicators of the stability properties of non-linear, time-dependent differential equations [11]. Therefore, no simplifications are required concerning the time-dependence of the system, paving the way for a more versatile and accurate assessment of rotorcraft stability.

Concerning their application to rotorcraft aeromechanics, LCEs and their sensitivity were studied for the first time for specific non-linear engineering systems and subcomponents in Refs. [12,13]. This work extends their use to time-dependent high fidelity rotorcraft problems as a quantitative stability assessment method. The rest of the paper is organized as follows. Section 2 briefly introduces Lyapunov Characteristic Exponents and their practical computation. Two time-dependent example problems are utilized in Section 3 to illustrate the approach, where the LCE stability indicators are compared with the conventional Floquet multiplier. Conclusions are finally drawn in Section 4.

2. Method

This section illustrates the theory and estimation of Lyapunov Characteristic Exponents (or Lyapunov Exponents, in short) and their sensitivity.

2.1. Lyapunov Characteristic Exponents

An LTI dynamical system can be represented in state-space form as:

$$\begin{array}{cccc}\hfill \dot{\mathbf{x}}& =\mathbf{A}\mathbf{x},\hfill & \hfill \mathbf{x}\left(t_0\right)& ={\mathbf{x}}_0,\hfill \end{array}$$

(1)

where vector $\mathbf{x}\in {\mathbb{R}}^n$ contains the n states of the system, and matrix $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is referred to as the state-space matrix. When matrix $\mathbf{A}$ is not singular, $\mathbf{x}=\mathbf{0}$ (the trivial solution) is the unique equilibrium condition. In this case, the stability of the LTI system in Equation (1) is estimated by the eigenvalues of the state-space matrix $\mathbf{A}$ [14]. Three fundamental outcomes may originate from this consideration.

• First, the time response of the system converges to zero, which happens when all eigenvalues have a negative real part.
• Second, the system responds in a divergent manner when some eigenvalues have a positive real part.
• Third, a periodic orbit can be the resulting motion in the presence of purely imaginary eigenvalues.

A combination of these fundamental response types can occur, as well [15]. Since only one equilibrium point exists, which is trivial, the stability properties of all possible solutions of $\mathbf{x}$ are equivalent to those of the trivial solution. Although the time response amplitude of the system depends on the initial conditions [16], the stability characteristics, i.e., being divergent, convergent, or periodic, do not depend on them.

For LTP systems, the same remarks are valid, as well. The only difference is the need to consider, instead of the state-space matrix $\mathbf{A}$, the monodromy matrix $\mathbf{H}$, which represents the state transition matrix (STM) over one period (often referred to as Floquet transition matrix), i.e., the result of integrating the problem

$$\begin{array}{cc}\hfill \dot{\mathbf{H}}& =\mathbf{A}\left(t\right)\mathbf{H}\hfill \end{array}$$

(2)

over a period from $t=0$ to T, with initial conditions $\mathbf{H}\left(0\right)=\mathbf{I}$. For an LTI system, the monodromy matrix over the same period T corresponds to

$$\begin{array}{cc}\hfill \mathbf{H}& ={\mathrm{e}}^{\mathbf{A}T}.\hfill \end{array}$$

(3)

While LTI and LTP systems can usually provide sufficiently accurate representations of rotorcraft problems for several applications, the relation between the system states $\mathbf{x}$ and its rate of change in time $\dot{\mathbf{x}}$ cannot be limited to being time-invariant or time-periodic. In other words:

$$\begin{array}{cccc}\hfill \dot{\mathbf{x}}& =\mathbf{f}(\mathbf{x},t),\hfill & \hfill \mathbf{x}\left(t_0\right)& ={\mathbf{x}}_0\hfill \end{array}$$

(4)

is a more generic representation of rotorcraft dynamics, where the non-linear function $\mathbf{f}(\mathbf{x},t)$ governs the system dynamics and can include relations not necessarily time-invariant nor periodic.

To address the more general time-dependence problem, this work suggests the use of Lyapunov Characteristic Exponents (LCEs). LCEs indicate the stability properties of the solutions of differential equations [11,14]. Mathematically speaking, LCEs define the spectrum of the corresponding initial value, or Cauchy, problem. Moreover, Lyapunov’s theory applies also to non-linear time-dependent systems. The stability of trajectories in the state-space can be estimated while integrating their time response. Therefore, Lyapunov Exponents represent a generalization of the stability properties that are familiar in current engineering practice.

By definition, the trajectory of a system is a solution to a differential equation of the form of Equation (4). Special cases occur if the model is linear, i.e., $\mathbf{f}(\mathbf{x},t)=\mathbf{A}\left(t\right)\mathbf{x}\left(t\right)$, or time-periodic, i.e., linear with $\mathbf{A}(t+T)=\mathbf{A}\left(t\right)$, for a given period T. Time-invariant problems arise when $\mathbf{f}\left(\mathbf{x}\right)$ does not explicitly depend on time t. The simplest form is obtained when the problem is simultaneously linear and time-invariant, i.e., $\mathbf{f}\left(\mathbf{x}\right)=\mathbf{A}\mathbf{x}$.

Consider the system with $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$, where the previously introduced vector $\mathbf{x}\in {\mathbb{R}}^n$ contains the states, $t\in \mathbb{R}$ is the time, and $\mathbf{f}\in {\mathbb{R}}^{n+1}\to {\mathbb{R}}^n$ is a non-linear time-dependent function. For a solution $\mathbf{x}\left(t\right)$ originating from initial conditions $\mathbf{x}\left(0\right)={\mathbf{x}}_0$, the LCEs ${\lambda }_i$ are defined as:

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{t\to \infty }{lim}\frac{1}{t}log\frac{{\left|\right|}_i\mathbf{x}\left(t\right)\left|\right|}{{\left|\right|}_i\mathbf{x}\left(t_0\right)\left|\right|},\hfill \end{array}$$

(5)

where ${}_i\mathbf{x}\left(t\right)$ represents the growth of the ith axis of the ellipsoid that evolves from an initially infinitesimal n-dimension sphere. The evolution happens according to the map ${\mathbf{f}}_{/\mathbf{x}}$, tangent to $\mathbf{f}$ along the solution trajectory $\mathbf{x}\left(t\right)$. In other words, ${}_i\mathbf{x}\left(t\right)$ is the solution of the linear, time-dependent problem

$$\begin{array}{cc}\hfill {}_i\dot{\mathbf{x}}\left(t\right)& ={\mathbf{f}}_{/\mathbf{x}}(\mathbf{x}\left(t\right),t){}_i\mathbf{x}\left(t\right),\hfill \end{array}$$

(6)

with initial conditions ${}_i\mathbf{x}\left(t_0\right)={}_i{\mathbf{x}}_0$.

Equation (5) can be simplified by splitting the log function of Equation (5) as:

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{t\to \infty }{lim}\frac{1}{t}log{\left|\right|}_i\mathbf{x}\left(t\right)\left|\right|-\underset{t\to \infty }{lim}\frac{1}{t}log{\left|\right|}_i\mathbf{x}\left(t_0\right)\left|\right|.\hfill \end{array}$$

(7)

The second limit involves a constant, ${\left|\right|}_i\mathbf{x}\left(t_0\right)\left|\right|$, divided by time t, which, in the limit, approaches $+\infty$, obviously tending to zero. The first limit, thus, defines the ith LCE as:

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{t\to \infty }{lim}\frac{1}{t}{log\left|\right|}_i\mathbf{x}\left(t\right)\left|\right|.\hfill \end{array}$$

(8)

Due to the presence of a limit for $t\to \infty$ in the definition, the practical evaluation of LCEs can only be numerically obtained for a sufficiently long time t. For this reason, the term “LCEs” in the following will rather refer to their estimation for a large enough value of t to reach practical convergence.

2.2. Application to LTP Problems

For an LTV problem, the STM is the solution of

$$\begin{array}{cc}\hfill \dot{\mathbf{Y}}& =\mathbf{A}\left(t\right)\mathbf{Y},\hfill \end{array}$$

(9)

integrated from $\tau$ to t with initial conditions, $\mathbf{Y}\left(\tau \right)=\mathbf{I}$. When replaced in the LCEs definition, it yields

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{t\to \infty }{lim}\frac{1}{t}log\left|\right|\mathbf{Y}(t,0){}_i{\mathbf{x}}_0\left|\right|.\hfill \end{array}$$

(10)

When the problem is LTP, the time t can be nondimensionalized with respect to the problem’s period T, namely $t=T\psi$, yielding

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{\psi \to \infty }{lim}\frac{1}{T\psi }log\left|\right|\mathbf{Y}(T\psi ,0){}_i{\mathbf{x}}_0\left|\right|.\hfill \end{array}$$

(11)

Considering now the limit for discrete values of time, for $\psi \in {\mathbb{N}}^{+}$ to $+\infty$, one obtains $\mathbf{Y}(T\psi ,0)={\mathbf{Y}}^{\psi }(T,0)$, and, recalling the monodromy matrix resulting from the solution of Equation (2), $\mathbf{H}=\mathbf{Y}(T,0)$, one obtains

$$\begin{array}{c}\hfill {\lambda }_i=\underset{\psi \to \infty }{lim}\frac{1}{T\psi }log\left|\right|{\mathbf{H}}^{\psi }{}_i{\mathbf{x}}_0\left|\right|,\end{array}$$

(12)

and, again, expressing the monodromy matrix in spectral form, namely $\mathbf{H}=\mathbf{X}{\mathrm{e}}^{\mathbf{B}T}{\mathbf{X}}^{-1}$, and assuming that ${}_i{\mathbf{x}}_0$ is its ith eigenvector, namely the ith column of matrix $\mathbf{X}$, one obtains

$$\begin{array}{c}\hfill {\lambda }_i=\underset{\psi \to \infty }{lim}\frac{1}{T\psi }log\left|\right|{}_i{\mathbf{x}}_0{\mathrm{e}}^{{\beta }_{i}T\psi }\left|\right|=\mathrm{Re}\left({\beta }_i\right),\end{array}$$

(13)

i.e., the LCEs correspond to the real part of the logarithm of the eigenvalues of the monodromy matrix, divided by the period, T.

2.3. Numerical Estimation of LCEs

Starting from the fundamental work of Benettin et al. [11], several methods have been formulated for the estimation of all LCEs, or at least the largest ones. They can be broadly partitioned into two classes: continuous and discrete. Noteworthy continuous methods include those based on the QR and SVD factorizations, whereas, among the discrete ones, that based on QR is worth mentioning, and will be considered in the rest of this work.

LCE evaluation using continuous formulas suffers from the computational difficulty of dealing with matrices whose coefficients either diverge or go to zero rapidly. As a result, practical formulations have been developed. For example, the discrete QR method is a common technique in LCE calculation. In this method, the LCEs estimates are incrementally updated with the diagonal elements of the matrix $\mathbf{R}$, which is extracted using the QR decomposition of the STM evaluated between two successive time steps.

For a given matrix $\mathbf{Y}(t,t_{j-1})$ from time $t_{j-1}$ to the consecutive time t, define the state transition of the problem $\dot{\mathbf{Y}}={\mathbf{f}}_{/\mathbf{x}}(\mathbf{x}\left(t\right),t)\mathbf{Y}$ with $\mathbf{Y}\left(t_{j-1}\right)=\mathbf{I}$. Setting ${\mathbf{Y}}_j=\mathbf{Y}(t_j,t_{j-1})$ and considering the QR decomposition of ${\mathbf{Y}}_j{\mathbf{Q}}_{j-1}$ implies that ${\mathbf{Q}}_j{\mathbf{R}}_j={\mathbf{Y}}_j{\mathbf{Q}}_{j-1}$. When the ${\mathbf{R}}_{{\Pi }_j}$ matrix is defined as ${\mathbf{R}}_{{\Pi }_j}={\Pi }_{k=0}^j{\mathbf{R}}_{j-k}$, it can be shown that:

$$\begin{array}{cc}\hfill {\mathbf{Y}}_j{\mathbf{Q}}_{j-1}{\mathbf{R}}_{{\Pi }_{j-1}}& ={\mathbf{Q}}_j{\mathbf{R}}_j{\mathbf{R}}_{{\Pi }_{j-1}}={\mathbf{Q}}_j{\mathbf{R}}_{{\Pi }_j}.\hfill \end{array}$$

(14)

Using the above relation, ${\mathbf{Y}}_j{\mathbf{Q}}_{j-1}{\mathbf{R}}_{{\Pi }_{j-1}}$ ${\mathbf{R}}_{{\Pi }_j}$ can be constructed using only incremental QR decompositions over ${\mathbf{Y}}_j{\mathbf{Q}}_{j-1}$, which allows limited contraction/expansion over one time-step. Now, the LCEs can be estimated from the ${\mathbf{R}}_{{\Pi }_j}$ matrix as

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{j\to \infty }{lim}\frac{1}{t_j}log{r}_{ii}\left(t_j\right),\hfill \end{array}$$

(15)

where $j\in \mathbb{N}$ and $r_{ii}\left(t_j\right)$ are the diagonal elements of matrix $\mathbf{R}\left(t_j\right)={\mathbf{R}}_{{\Pi }_j}$. Since the product of two upper triangular matrices $\mathbf{C}=\mathbf{A}\mathbf{B}$ is also an upper triangular matrix, the diagonal elements are $c_{ii}=a_{ii}{b}_{ii}$. Thus, the $log\left(c_{ii}\right)$ can be computed incrementally, namely $log\left(a_{ii}{b}_{ii}\right)=log\left(a_{ii}\right)+log\left(b_{ii}\right)$. This helps prevent overflow/underflow in numerical computations. Furthermore,

$$\begin{array}{cc}\hfill r_{ii}\left(t_j\right)& ={\Pi }_{k=0}^j{r}_{{(j-k)}_{ii}};\hfill \end{array}$$

(16)

thus,

$$\begin{array}{cc}\hfill log\left(r_{ii}\left(t_j\right)\right)& =\sum _{k=0}^{j}log\left(r_{k_{ii}}\right),\hfill \end{array}$$

(17)

which leads to

$$\begin{array}{cc}\hfill {\lambda }_i& =\underset{j\to \infty }{lim}\frac{1}{t_j}\sum _{k=0}^{j}log\left(r_{k_{ii}}\right).\hfill \end{array}$$

(18)

2.4. Computation of State Transition Matrix

The discrete QR decomposition method needs the calculation of the state transition matrix, which can be obtained by numerical integration. For a sufficiently small time-step, Equation (4) can be linearized:

$$\begin{array}{cc}\hfill \delta \dot{\mathbf{x}}& =\mathbf{A}\left(\mathbf{x}\right(t),t)\delta \mathbf{x},\hfill \end{array}$$

(19)

where $\mathbf{A}={\mathbf{f}}_{/\mathbf{x}}$ results from the partial derivative of the non-linear function $\mathbf{f}$ with respect to the state-space variables $\mathbf{x}$. On the other hand, $\mathbf{f}$ can be an arbitrary non-linear function of ${\mathbf{x}}_j$ (trajectory) and t (time). Therefore, to integrate the STM, the knowledge of the fiducial trajectory is required, i.e., the equations of motion need to be integrated first. The numerical integration of differential equations is a relatively straightforward task; as such, it is not discussed here in detail. However, the computation of the state transition matrix is presented, as its determination is essential for the procedure.

An efficient approach for the evaluation of the state transition matrix is Hsu’s method [17], which is also reported in Reference [18] and followed in this work. The method is applied to LTV problems of the form $\dot{\mathbf{x}}=\mathbf{A}\left(t\right)\mathbf{x}$ by considering a piecewise-constant approximation of matrix $\mathbf{A}\left(t\right)$, namely $\mathbf{A}(\mathbf{x},t)\approx \mathbf{A}(\widehat{\mathbf{x}},\widehat{t})=\widehat{\mathbf{A}}$ with $\widehat{t}\in [t_j,t_{j+1}]$, where $t_j\le t\le t_{j+1}$ and $\widehat{\mathbf{x}}=\mathbf{x}\left(\widehat{t}\right)$. The selection of $\widehat{t}$ may slightly alter the resulting STM. Finally, the state transition matrix is readily obtained as

$$\begin{array}{cc}\hfill \mathbf{Y}(t,t_j)& \approx {\mathrm{e}}^{\widehat{\mathbf{A}}(t-t_j)},\hfill \end{array}$$

(20)

where an approximation of the matrix exponential may be necessary to improve the computational efficiency, such as truncating the matrix power series.

2.5. Analytical Sensitivity of LCEs

In a dynamical system, especially in the design phase, the rate of change of stability estimates with respect to an evolving or uncertain parameter plays a significant role in assessing the stability characteristics. Such sensitivity is useful to gain insight into the dependence of stability indicators on system parameters or can be integrated into gradient-based (or gradient-aware) optimization procedures [19,20] and continuation algorithms [21], or into uncertainty evaluation problems. Rather than using finite-differences, analytical sensitivity is more favorable for several reasons, such as: to avoid issues related to sharp changes in sensitivity parameters, and to track the evolution of the stability indicators using continuation algorithms.

Consider a parameter p as a member of the set of bounded parameters $\mathbf{p}\in \mathcal{P}$. In addition, suppose that the system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t,\mathbf{p})$ depends on this parameter set. The sensitivity of the LCEs with respect to the parameter $p\subset \mathbf{p}$, using the form of Equation (15), can be stated as:

$$\begin{array}{cc}\hfill {\lambda }_{i/p}& =\underset{t\to \infty }{lim}\frac{1}{t}\frac{r_{ii/p}\left(t\right)}{r_{ii}\left(t\right)}.\hfill \end{array}$$

(21)

The ability to compute the sensitivity of matrix $\mathbf{R}$, namely ${\mathbf{R}}_{/p}$, is needed. Considering the definition of $\mathbf{R}\left(t_j\right)={\mathbf{R}}_{{\Pi }_j}$ as formulated within the discrete QR method, its sensitivity is:

$$\begin{array}{c}\hfill {\mathbf{R}}_{/p}\left(t_j\right)=\sum _{k=0}^j\left({\Pi }_{n=0}^{k-1}{\mathbf{R}}_{j-n}\right){\mathbf{R}}_{k/p}\left({\Pi }_{n=k+1}^j\right){\mathbf{R}}_{j-n}.\end{array}$$

(22)

Since $\mathbf{R}\left(t_j\right)={\mathbf{R}}_j\mathbf{R}\left(t_{j-1}\right)$, the ${\mathbf{R}}_{/p}$ can be computed incrementally according to:

$$\begin{array}{cc}\hfill {\mathbf{R}}_{/p}\left(t_j\right)& ={\mathbf{R}}_{j/p}\mathbf{R}\left(t_{j-1}\right)+{\mathbf{R}}_j{\mathbf{R}}_{/p}\left(t_{j-1}\right);\hfill \end{array}$$

(23)

thus, the perturbation is accumulated straightforwardly.

Finally, the form of Equation (18) can be used to formulate the sensitivity of LCEs:

$$\begin{array}{cc}\hfill {\lambda }_{i/p}& =\underset{j\to \infty }{lim}\frac{1}{t_j}\sum _{k=1}^j\frac{r_{k_{ii}/p}}{r_{k_{ii}}},\hfill \end{array}$$

(24)

which only needs the sensitivity of ${\mathbf{R}}_j$, which is explained next.

2.5.1. Sensitivity of QR Decomposition

The sensitivity of the elements of QR decomposition can be computed along the lines of the state transition matrix computation with differentiation of the QR decomposition, which is utilized to formulate the continuous QR method for LCE estimation [11,22]. Revisiting the QR decomposition of an arbitrary matrix $\mathbf{M}\in {\mathbb{R}}^{n\times n}$:

$$\begin{array}{cc}\hfill \mathbf{M}& =\mathbf{Q}\mathbf{R},\hfill \end{array}$$

(25)

with orthogonal $\mathbf{Q}\in {\mathbb{R}}^{n\times n}$ (${\mathbf{Q}}^T\mathbf{Q}=\mathbf{I}$), and upper triangular $\mathbf{R}\in {\mathbb{R}}^{n\times n}$ having positive diagonal elements. Differentiation of $\mathbf{M}$ with respect to a scalar parameter p gives:

$$\begin{array}{cc}\hfill {\mathbf{M}}_{/p}& ={\mathbf{Q}}_{/p}\mathbf{R}+\mathbf{Q}{\mathbf{R}}_{/p},\hfill \end{array}$$

(26)

which requires the derivative of the orthogonality condition ${\mathbf{Q}}^T\mathbf{Q}=\mathbf{I}$:

$$\begin{array}{c}\hfill {\left({\mathbf{Q}}^T\right)}_{/p}\mathbf{Q}+{\mathbf{Q}}^T{\mathbf{Q}}_{/p}={\left({\mathbf{Q}}^T{\mathbf{Q}}_{/p}\right)}^T+{\mathbf{Q}}^T{\mathbf{Q}}_{/p}=\mathbf{0}.\end{array}$$

(27)

The second relation implies that ${\mathbf{Q}}^T{\mathbf{Q}}_{/p}$ must be skew-symmetric; therefore, only $n(n-1)/2$ coefficients are independent (for example, those in the lower triangular part, excluding the diagonal). Finally, ${\mathbf{M}}_{/p}$ is pre-multiplied by ${\mathbf{Q}}^T$:

$$\begin{array}{cc}\hfill {\mathbf{R}}_{/p}& ={\mathbf{Q}}^T{\mathbf{M}}_{/p}-{\mathbf{Q}}^T{\mathbf{Q}}_{/p}\mathbf{R}.\hfill \end{array}$$

(28)

The ${\mathbf{R}}_{/p}$ matrix is upper triangular; hence, the whole problem can be re-cast in the following steps:

$$\begin{array}{cccc}\hfill \mathrm{compute}& \hfill & & \mathbf{W}={\mathbf{Q}}^T{\mathbf{Q}}_{/p}\hfill \end{array}$$

(29a)

$$\begin{array}{cccc}\hfill \mathrm{considering}& \hfill & & \mathrm{stril}\left({\mathbf{Q}}^T{\mathbf{M}}_{/p}-\mathbf{W}\mathbf{R}\right)=\mathrm{stril}\left(\mathbf{0}\right)\hfill \end{array}$$

(29b)

$$\begin{array}{cccc}\hfill \mathrm{constrained}\mathrm{by}& \hfill & & {\mathbf{W}}^T+\mathbf{W}=\mathbf{0},\hfill \\ \hfill \end{array}$$

(29c)

$$\begin{array}{cccc}\hfill \mathrm{compute}& \hfill & & {\mathbf{R}}_{/p}\hfill \end{array}$$

(29d)

$$\begin{array}{cccc}\hfill \mathrm{considering}& \hfill & & \mathrm{triu}\left({\mathbf{R}}_{/p}\right)=\mathrm{triu}\left({\mathbf{Q}}^T{\mathbf{M}}_{/p}-\mathbf{W}\mathbf{R}\right)\hfill \end{array}$$

(29e)

$$\begin{array}{cccc}\hfill \mathrm{and}& \hfill & & \mathrm{stril}\left({\mathbf{R}}_{/p}\right)=\mathrm{stril}\left(\mathbf{0}\right),\hfill \end{array}$$

(29f)

where the $\mathrm{triu}(·)$ extracts the upper triangular part, including the diagonal elements, and the $\mathrm{stril}(·)$ extracts the strictly lower triangular part, excluding the diagonal elements.

It should be noted that, since $\mathbf{R}$ is upper triangular, $\mathbf{W}$ is computed as:

$$\begin{array}{cccc}\hfill {\mathbf{W}}_L& =\mathrm{stril}\left({\mathbf{Q}}^T{\mathbf{M}}_{/p}{\mathbf{R}}^{-1}\right)\hfill & \hfill \mathbf{W}& ={\mathbf{W}}_L-{\mathbf{W}}_L^T,\hfill \end{array}$$

(30)

where the matrix ${\mathbf{R}}^{-1}$ can be computed by back-substitution. In fact, since $\mathbf{B}={\mathbf{Q}}^T{\mathbf{M}}_{/p}$, the generic coefficient of $\mathrm{stril}\left(\mathbf{W}\right)$ is

$$\begin{array}{cccccc}\hfill w_{ij}& =\frac{1}{r_{jj}}\left(b_{ij}-\sum _{k=1}^{j-1}{w}_{ik}{r}_{kj}\right)\hfill & \hfill j& =1,n-1\hfill & \hfill i& =j+1,n,\hfill \end{array}$$

(31)

leading to:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_{/p}& =\mathbf{Q}\mathbf{W}.\hfill \end{array}$$

(32)

The decomposition of ${\mathbf{Y}}_j{\mathbf{Q}}_{j-1}$ is required here; thus, the sensitivity of $\mathbf{M}={\mathbf{Y}}_j{\mathbf{Q}}_{j-1}={\mathbf{Q}}_j{\mathbf{R}}_j$ is computed, i.e.

$$\begin{array}{c}\hfill {\mathbf{M}}_{/p}={\mathbf{Q}}_{j/p}{\mathbf{R}}_j+{\mathbf{Q}}_j{\mathbf{R}}_{j/p}={\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{Y}}_j{\mathbf{Q}}_{(j-1)/p},\end{array}$$

(33)

where ${\mathbf{Q}}_{j-1}$ and ${\mathbf{Q}}_{(j-1)/p}$ are already known from the previous step. Next, Equation (33) is premultiplied by ${\mathbf{Q}}_j^T$ to get:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_j^T{\mathbf{Q}}_{j/p}{\mathbf{R}}_j+{\mathbf{R}}_{j/p}& ={\mathbf{Q}}_j^T{\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{Q}}_j^T{\mathbf{Y}}_j{\mathbf{Q}}_{(j-1)/p}.\hfill \end{array}$$

(34)

Then, the strictly lower triangular part of the equation is computed to obtain ${\mathbf{W}}_L$,

$$\begin{array}{cc}\hfill {\mathbf{W}}_L& =\mathrm{stril}\left(\left({\mathbf{Q}}_j^T{\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{Q}}_j^T{\mathbf{Y}}_j{\mathbf{Q}}_{(j-1)/p}\right){\mathbf{R}}_j^{-1}\right)\hfill \\ \hfill & =\mathrm{stril}\left(\left({\mathbf{Q}}_j^T{\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{R}}_j{\mathbf{W}}_{j-1}\right){\mathbf{R}}_j^{-1}\right).\hfill \end{array}$$

(35)

The strictly lower triangular part of ${\mathbf{W}}_j$ is ${\mathbf{W}}_j={\mathbf{Q}}_j^T{\mathbf{Q}}_{j/p}={\mathbf{W}}_L-{\mathbf{W}}_L^T$. Finally, the upper triangular part of Equation (34) is computed to attain ${\mathbf{R}}_{j/p}$,

$$\begin{array}{cc}\hfill {\mathbf{R}}_{j/p}& ={\mathbf{Q}}_j^T\left({\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{Y}}_j{\mathbf{Q}}_{(j-1)/p}\right)-{\mathbf{W}}_j{\mathbf{R}}_j\hfill \\ \hfill & ={\mathbf{Q}}_j^T{\mathbf{Y}}_{j/p}{\mathbf{Q}}_{j-1}+{\mathbf{R}}_j{\mathbf{W}}_{j-1}-{\mathbf{W}}_j{\mathbf{R}}_j,\hfill \end{array}$$

(36)

which depends on computation of the sensitivity of ${\mathbf{Y}}_j$, from $t_{j-1}$ to $t_j$, as explained next.

2.5.2. Sensitivity of the State Transition Matrix

The sensitivity of the state transition matrix $\mathbf{Y}$ at time $t_j$, namely ${\mathbf{Y}}_{/p}\left(t_j\right)$, is finally required to complete the computation of sensitivity of the LCEs. In principle, one can achieve this by integrating the sensitivity of the problem. $\dot{\mathbf{Y}}=\mathbf{A}\mathbf{Y}$. Then,

$$\begin{array}{cc}\hfill {\dot{\mathbf{Y}}}_{/p}& =\mathbf{A}{\mathbf{Y}}_{/p}+\left({\mathbf{A}}_{/\mathbf{x}}{\mathbf{x}}_{/p}+{\mathbf{A}}_{/p}\right)\mathbf{Y},\hfill \end{array}$$

(37)

i.e., a problem with the same matrix $\mathbf{A}$ of the original one, forced by a term $(\mathrm{d}\mathbf{A}/\mathrm{d}p)\mathbf{Y}$ that depends on the reference solution. The term ${\mathbf{A}}_{/\mathbf{x}}$ is the second-order derivative of function $\mathbf{f}$ with respect to the state $\mathbf{x}$. Therefore, the solution becomes independent for the change in the trajectory (${\mathbf{x}}_{/p}$) for LTP problems. Finally, the sensitivity of the state transition matrix can be computed using Hsu’s method [17] for the state transition matrix is readily obtained as $\mathbf{Y}(t,t_j)\approx {\mathrm{e}}^{\widehat{\mathbf{A}}(t-t_j)}$:

$$\begin{array}{cc}\hfill {\mathbf{Y}}_{/p}(t,t_j)& \approx {\widehat{\mathbf{A}}}_{/p}(t-t_j){\mathrm{e}}^{\widehat{\mathbf{A}}(t-t_j)}={\widehat{\mathbf{A}}}_{/p}(t-t_j)\mathbf{Y}(t,t_j).\hfill \end{array}$$

(38)

3. Numerical Results

This section presents the verification of the method using an analytical solution and later demonstrates the method over a complex multidisciplinary rotorcraft model.

3.1. Verification with an Analytical Solution

Before providing the detailed rotorcraft model and its stability analysis using LCEs, the proposed procedure was verified against analytical solution. Consider the homogeneous mass-damper system with a periodic damping term:

$$\begin{array}{c}\hfill m\ddot{x}+(c_0+c_p{cos}^{2}t)\dot{x}=0,\end{array}$$

(39)

where the system has a period of $\pi$ radians. Then, the system can be considered as a first order system with $q=\dot{x}$.

$$\begin{array}{c}\hfill \dot{q}=-\frac{(c_0+c_p{cos}^{2}t)}{m}q=0\end{array}$$

(40)

The monodromy matrix (in this case a scalar) can be obtained by integrating Equation (40) from an initial time $t=0$ to $t=T$. Let $m=1$;

$$\begin{array}{c}\hfill \frac{q_T}{q_0}=H={\mathrm{e}}^{-{\int }_0^T\left[c_0+c_p{cos}^{2}t\right]dt}={\mathrm{e}}^{-T{c}_0-\frac{T}{2}{c}_p}\end{array}$$

(41)

The eigenvalue of the system $\lambda$ can be obtained from the eigenvalue of the monodromy matrix:

$$\begin{array}{c}\hfill \lambda =\frac{1}{T}log\left(H\right)=\frac{1}{T}log\left({\mathrm{e}}^{-T{c}_0-\frac{T}{2}{c}_p}\right)=-c_0-\frac{1}{2}{c}_p\end{array}$$

(42)

If the coefficient of the periodic term, $c_p$, is considered as the parameter; the sensitivity of the eigenvalue, ${\lambda }_p$, can be analytically obtained as:

$$\begin{array}{c}\hfill {\lambda }_{/p}=-\frac{1}{2}.\end{array}$$

(43)

The analytical results are compared with the LCE numerical solution for $c_0=c_p=1$. Figure 1 presents the characteristic exponent obtained by the analytical and numerical (named Floquet), and Figure 2 gives the corresponding sensitivity of the characteristic exponents. Results show good correlation.

3.2. Complex LTP Rotorcraft Model

The nature of rotary-wing stability presents the interaction of structural dynamics, aerodynamics, and control systems to model rotorcraft components, such as rotors, fuselage, gearbox, etc. Each component may require a considerable number of degrees of freedom to describe the problem within the desired accuracy, even in its simplest and most compact form. As discussed earlier in this paper, LCEs estimation requires matrix operations at each time-step, including multiplications and orthogonalizations. Theoretically, this procedure should be carried on for an infinite limit ($t\to \infty$), and in practice till LCEs converge to a steady finite value, a criterion that may not be verifiable in a trivial manner. This may require significant computational cost, compared to other classical methods that are otherwise limited to LTI or LTP problems [12].

In problems described with a relatively large number of degrees of freedom, only extracting the most critical LCEs or the few largest ones closest to zero can help improve computation efficiency [23,24]. However, in a detailed rotorcraft model, many lightly damped LCEs may appear, especially when the aircraft’s rigid body motion and fuselage elastic degrees of freedoms are considered. It may be impossible to discern between lightly damped modes that are not going to jeopardize stability and modes that, despite being more damped in the reference configuration, will eventually turn unstable for some variations of the parameters of the problem. Therefore, within the scope of this work, even when applied to a detailed problem, the method needs to be verified for the whole spectrum.

3.2.1. Model Description

To illustrate the application of the proposed approach to a complex problem, an aeroelastic model of a rotorcraft was developed, based on data for the AS-330 PUMA helicopter as presented in Reference [25]. Its general characteristics are reported in

Table 1. AS-330 PUMA general characteristics [25].

Parameter	Value	Units
Helicopter
Gross Weight	7400	kg
Max Speed	140	kn
Main Rotor
Number of blades	4
Radius	7.49	m
Solidity	0.0913	(n.d.)
Lock number	8.70	(n.d.)
Speed	270	rpm
Flap frequency	1.03	/rev
Lag Frequency	0.26	/rev

. The solidity parameter is the ratio of the actual rotor area over the area of the disk,

$$\begin{array}{cc}\hfill \mathrm{Solidity}& =\frac{N_{b}c}{\pi R},\hfill \end{array}$$

(44)

where $N_b$ is the number of blades, and c is the chord, whereas the Lock number expresses the ratio between the aerodynamic and the inertia moments about the flap hinge that are associated with the flapping motion of the blades,

$$\begin{array}{cc}\hfill \mathrm{Lock}\mathrm{number}& =\frac{\rho ac{R}^4}{I_{\beta }},\hfill \end{array}$$

(45)

where $\rho$ is the air density, a the slope of the lift coefficient curve (often approximated as $a=2\pi$ according to thin airfoil theory), and $I_{\beta }$ the moment of inertia of the blade about the flap hinge. These parameters are used to nondimensionalize the dynamics and aeromechanics of the main rotor (see, for example, Reference [2]). The overall model comprises:

• all six rigid body modes (Fore/Aft, Lateral, Plunge, Roll, Pitch, and Yaw);
• flight mechanics derivatives of the fuselage determined by fuselage/wing-body, horizontal tail, and the vertical tail;
• ten elastic airframe modes captured at airframe and connection locations (such as airframe rotor connection); additionally, $1.5$% modal damping was also added;
• three bending modes of the rotor in multiblade coordinates, formulated using a linerazized finite volume approach [26];
• three main rotor servo-actuators modelled as transfer function from the force applied by controls ($f_c$) and requested displacement ($x_c$) to the servoactuator displacement ($x_s$), namely $x_s=H_x{x}_c+H_f{f}_c$;
• lead-lag dampers for the main rotor blades.

The above components are blended in MASST [27], an aeroservoelastic simulation platform with proven capabilities for rotorcraft comprehensive modeling [28,29]. The helicopter model is developed for hover flight, so the rotor model is isotropic; therefore, the model is linear time-invariant under normal conditions. Periodicity is introduced by considering the characteristic coefficient of the lead-lag damper of one blade as a parameter, varying from zero ($C_L=0$) to twice its nominal value ($C_L=2C_d$), while the other dampers maintain their nominal value ($C_L=C_d$). The resulting overall system of equations can be stated as:

$$\begin{array}{c}\hfill \mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}(t,T)\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}=\mathbf{0},\end{array}$$

(46)

where T is the period of the system, namely reciprocal of the rotor angular speed $T=1/\Omega$.

While the mass matrix $\mathbf{M}$ and stiffness matrix $\mathbf{K}$ are constant in time, the damping matrix $\mathbf{C}$ is periodic due to the asymetric distribution of lead-lag dampers. The resulting model is LTP with respect to the angular motion of the rotor. Stability and sensitivity analyses of the LTP model are performed using Lyapunov’s theory, as explained in a previous section. To illustrate the effectiveness of the proposed method, LCEs and their sensitivity are estimated and compared to that of the Floquet multipliers and their analytical sensitivities using the formulation in Ref. [30].

3.2.2. Analysis

Figure 3 presents the characteristic exponents resulting from both Floquet’s and Lyapunov’s analysis in the parameter range of interest. The LCEs show a good correlation with Floquet’s exponents. One can also observe that the higher modes (larger absolute valued LCEs) are not sensitive to the parameter change. However, some lightly damped modes that strongly interact with the in-plane motion of the rotor show higher sensitivity to the dissimilarity in the rotor blades’ dynamics. This can be better observed in Figure 4, which zooms on lightly damped modes.

Among the lightly damped modes, the one corresponding to regressive cyclic lag is the one that has the largest sensitivity to the parameter change. The corresponding LCEs and Floquet multipliers were separately plotted in Figure 5. It can be observed that the system does not lose the stability of the coupled airframe/lead-lag motion as one lead-lag damper reduces its effectiveness. Even when one damper is completely inoperative, the system is stable, which is most likely due to the presence of aerodynamic damping resulting from the coupled lag-flap motion.

The sensitivities corresponding to the LCEs and Floquet multipliers of Figure 5 are shown in Figure 6. Sensitivity estimates of the LCEs are compared with the corresponding sensitivity of the Floquet exponents and with analogous results obtained from numerical differentiation through finite differences. Results show that the sensitivity estimates obtained using LCEs are close enough to those obtained via finite differences, but not as close to those obtained through the sensitivity analysis of Floquet’s exponents.

4. Conclusions

This work presented the use of Lyapunov Characteristic Exponents (LCEs) for the estimation of the stability properties of time-dependent rotorcraft models. The theory of LCEs and the computation of their parameter sensitivity were presented with practical estimation methods. The proposed approach was applied to a comprehensive aeroelastic rotorcraft model. Time-variance was introduced by changing the properties of one of the lead-lag dampers, resulting in a time-periodic system. The resulting LCEs and their sensitivity were compared with those resulting from the Floquet method. The LCEs were in very good agreement with Floquet multipliers. Considering that the LCEs can also be generalized to non-linear problems, LCEs can represent a replacement and a generalization of the Floquet method. The analytical sensitivity estimation was not as successful, although found in relatively good agreement with finite-difference estimations. More efficient LCE sensitivity estimation algorithms are needed to improve that aspect of LCEs.