A Generalized State-Space Aeroservoelastic Model Based on Tangential Interpolation

Abstract

In this work, a new approach for the generation of a generalized state-space aeroservoelastic model based on tangential interpolation is presented. The resulting system of differential algebraic equations (DAE) is reduced to a set of ordinary differential equations (ODE) by residualization of the non-proper part of the transfer function matrix. The generalized state-space is of minimal order and allows for the application of the force summation method (FSM) for the aircraft loads recovery. Compared to the classical rational function approximation (RFA) approach, the presented method provides a minimal order realization with exact interpolation of the unsteady aerodynamic forces in tangential directions, avoiding any selection of poles (lag states). The new approach is applied first for the generation of an aerodynamic model for the bidimensional unsteady incompressible flow in the time domain. Next, an application on the generation of an aeroservoelastic model for loads evaluation of the flutter reduced order assessment (FERMAT) model under atmospheric disturbances is done, showing an excellent agreement with the reference model in the frequency domain. The proposed aeroservoelastic model of minimal order is suited for loads analysis and multivariable control design, and an application to a gust loads alleviation (GLA) strategy is shown.

1. Introduction

One of the major tasks for aircraft certification is the consideration of dynamic load cases, in particular those caused by gust and continuous turbulence encounters. The scenarios to be considered are specified by the regulatory agencies [1,2]. Typically, the sizing of the aircraft components is determined by different loadings corresponding to different excitations and both discrete gust and continuous turbulence are of major importance. In order to reduce the internal loads caused by these atmospheric disturbances without compromising the structural weight, different load alleviation schemes are applied. For modern transport aircraft configurations, an accurate model suited for the design of advance control laws covering the complete flight envelope is of most relevance, as the physical effects governed by the subsonic and transonic flow surrounding the aircraft structure must be precisely described. Thus, the need arises for an accurate but at the same time efficient model description over the complete flight envelope in order to design suitable control laws for an appropriate load reduction when encountering atmospheric disturbances. Moreover, modern transport aircraft are becoming more flexible and, as a result, the low frequency flexible modes tend to interact with the rigid-body behaviour of the structure. Thus, a model description containing the rigid and flexible effects should be taken into account for an effective control design.

Modern approaches of control theory require a description of the model in the time domain. However, the aerodynamic models describing the flow over the aircraft are provided within the frequency domain in tabular form as they cannot usually be obtained in an explicit form and, if they are, they contain transcendental functions which cannot be transformed in the time domain by a direct application of an inverse Laplace transform. Thus, there is a need for converting the frequency-domain description into time-domain models. For the application of the control theory techniques, it is moreover desirable that the size of the time-domain description is kept minimal.

Historically, the unsteady subsonic air forces for three-dimensional configurations have been determined for pure harmonic motion. For analytic functions in the frequency domain, the function in the whole complex plane defined by the Laplace variable s can be then deduced from the values over the imaginary axis $i\omega$ for pure harmonic motion [3]. The method of Roger [4] and Abel [5], based on a rational function approximation (RFA), represents the three-dimensional subsonic aerodynamic by means of a least-squares technique with a series of poles which represent the aerodynamic lags due to the presence of the wake. All elements of the transfer function matrix share the same poles, which in turn have to be chosen in advance. Note that the least-squares fit is ill-conditioned when increasing the number of real poles and is prone to numerical instabilities [6,7]. When applied at a half-generalized level in an effort not to increase the size of the state-space model, the least-squares fit for the gust excitation shows additional issues related to the time delay behaviour of the input. To overcome this problem, some authors divide the gust input in different zones along the flight direction and thereby increase the number of effective inputs [8,9]. Other authors apply an interpolation for the distributed aerodynamic loads, avoiding the delay term to be included in the frequency-domain description at the cost of drastically increasing the number of states [10,11]. In order to take into account the poles as free parameters to the least-squares fit, they can be further considered in an optimization problem [12,13], ensuring a minimum error of the least-squares for a fixed number of poles. Karpel presented the minimum-state method [14], where an iterative least-squares process is applied to the aeroelastic system, reducing the total number of augmented states. This iterative process leads to a significant increase of the computational effort which may be of up to three orders of magnitude [3].

Regarding the load alleviation strategies, the objective function to be minimized for a controller synthesis is a set of cut loads acting at specific aircraft locations. There are two common ways to compute the nodal loads acting over the airframe, namely the mode displacement method (MDM) and the force summation method (FSM) [15]. Several authors [15,16] showed that the FSM method has a superior convergence with respect to the number of structural modes retained for the projection of the aeroelastic equation. As discussed by Castrichini et al. [11], the MDM is not appropriate for cases where the loads are not smoothly distributed over the aircraft structural model, as the load distribution is not close to orthogonal for high frequency structural modes in that case, as they are characterized by significant spatial waviness. The structural modes selected for the projection of the aeroelastic equation correspond to the unloaded structure and therefore they are not influenced by the applied loads, leading to the selected modal basis to be not suitable to represent generic load distributions. It is common practice to consider the generalized and not the distributed aerodynamic forces for the aeroservoelastic system in the time domain in order to improve the least-squares solution by a reduced number of quantities to be fitted [8], leading to the application of the MDM for the cut loads recovery as a unique option, requiring the FSM method the knowledge of the aerodynamic loads distribution. As a consequence, modern techniques of the control theory are applied to aeroelastic systems where the loads recovery is done by means of the MDM [17], even though the FSM method has a superior convergence of the cut loads’ values with respect to the number of structural modes considered. For the FSM method, the challenge consists of properly modeling the unsteady distributed aerodynamic loads and this is achieved by the present method.

In general, all the RFA methods introduce an error in the description of the aerodynamic loads due to the least-squares fit. In order to avoid this error, different techniques based on rational interpolation may be used. More recently, the Eigensystem Realization Algorithm (ERA) algorithm [18] has been widely used for system identification and model reduction [19,20] in the time domain. As stated by Ma et al. [21], the ERA algorithm produces theoretically the same reduced-order models as balanced proper orthogonal decomposition (POD) [22] with no need of an adjoint system. It can also be regarded as the data-driven approximation to balanced truncation [23]. As a rational interpolation technique, it can be considered as a rational interpolant at infinite frequency [24]. There are two main difficulties associated with the ERA application. Firstly, its application is limited to proper systems and thus only the MDM method can be considered for the loads recovery. Secondly, if the system under consideration has a large number of inputs and outputs, then a big computational effort is required to compute the singular value decomposition (SVD) of the dense and large-size Hankel matrix [23].

For the consideration of linear time-invariant (LTI) multi-input multi-output (MIMO) systems with a non-proper transfer function matrix, a more general setting, namely that of the linear descriptor systems, is needed [25]. Additionally, the approach presented in this work is based on the Loewner framework, which is an extension of the rational interpolation at infinite frequency to a set of frequencies available at particular values [24].

In this work, a novel approach to generate an aeroservoelastic state-space model suited for loads computation and control design which overcomes the above limitations is presented. The proposed approach shall fulfill the following properties:

• There is no need for poles selection (lag states). On the one hand, the classical selection of real poles for the RFA techniques does not allow for describing phenomena which present resonance behaviour or peaks in the frequency-domain description at frequencies higher than zero. On the other hand, the set of real poles causes the RFA least-squares fit to be ill-conditioned, and care must be taken when increasing the number of poles. To overcome these limitations, within the present approach, neither a pole selection is required nor the transfer function is limited to a set of rational functions with real poles (see Section 3.1).

• It provides a small-size generalized state-space representation. The term generalized refers here to the fact that the theory of linear descriptor systems is needed for the present approach (see Section 2.1). This term should not be confused with the term generalized of the aeroelastic equation, where the physical equation is projected onto the set of generalized coordinates corresponding to the modes of the structure in vacuum. The proposed approach can be a regarded as providing a reduced order model (ROM) in the time domain, as it enables solving the aeroservoelastic system in a very efficient way. This is achieved by the condition of minimality of the rational interpolant within the Loewner framework theory [25].

• It is applicable to (input) delay systems, in particular when the excitation is due to gust disturbances. As described above, RFA techniques based on a least-squares fit of the frequency-domain data are not suited for a gust disturbance input. Two practical solutions in order to avoid this limitation of the RFA techniques, namely dividing the gust excitation in zones or applying a least-squares fit to the distributed aerodynamic nodal loads, may dramatically increase the size of the aeroelastic model in the time domain.

• It includes rigid-body modes, explicitly dealing with the singularity caused by the translational aircraft motion at zero frequency. This problem has been considered by Karpel et al. [26] in the frequency domain. In this work, the aerodynamic transfer function matrix is modified to include the derivatives of the translational motion previous to the time-domain realization.

• It is computationally efficient when generating the state-space. Unlike classical approaches where the precision is increased at the cost of an iterative approach [14], the current approach does not require any iterative process. In addition, the computational effort is drastically reduced compared to the ERA method by the consideration of tangential interpolation data within the Loewner pencil [25], avoiding the use of the dense and large-size Hankel matrix.

• It recovers the cut loads by means of the FSM. In order to achieve this, the unsteady aerodynamic loads distribution over the aircraft must be represented in the time domain. As described above, the FSM method is known to have a superior convergence for the cut loads prediction compared to the MDM method commonly used in applications for gust load alleviation (GLA) design [27].

In Section 2, the Loewner realization in connection with the theory of descriptor linear systems is presented and an application to the bidimensional unsteady incompressible flow over an airfoil is shown. A generalization to general three-dimensional aeroservoelastic systems is then presented in Section 3 with special focus on the representation of the distributed unsteady aerodynamic forces in the time domain. In Section 4, applications to the flutter reduced order assessment (FERMAT) model are shown and finally conclusions and future work are pointed out in Section 5.

2. Generalized Realization Problem

2.1. Tangential Interpolation

A linear time-invariant descriptor system ${\mathcal{S}}_n=\left(\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)$ is described in the time domain by a set of differential and algebraic equations (DAE), also denoted as a generalized state-space system:

(1)$\begin{array}{c}{\mathcal{S}}_n:\mathbf{E}\dot{\mathbf{x}}\left(t\right)=\mathbf{A}\mathbf{x}\left(t\right)+\mathbf{B}\mathbf{u}\left(t\right),\hfill \\ \mathbf{y}\left(t\right)=\mathbf{C}\mathbf{x}\left(t\right)+\mathbf{D}\mathbf{u}\left(t\right),\hfill \end{array}$

where $\mathbf{x}\left(t\right)\in {\mathbb{R}}^n$ is the state vector, $\mathbf{u}\left(t\right)\in {\mathbb{R}}^{n_u}$ the input vector, $\mathbf{E},\mathbf{A}\in {\mathbb{R}}^{n\times n}$, $\mathbf{B}\in {\mathbb{R}}^{n\times n_u}$, $\mathbf{C}\in {\mathbb{R}}^{n_y\times n}$ and $\mathbf{D}\in {\mathbb{R}}^{n_y\times n_u}$ are constant matrices with $\mathbf{E}$ possibly singular. In this work, regular systems are considered, that is, $\det \left(s\mathbf{E}-\mathbf{A}\right)\ne 0$, except for a finite number of eigenvalues (which may have an infinite value), denoted by the set $\sigma \left(\mathbf{E},\mathbf{A}\right)$. The resolvent set $\rho \left(\mathbf{E},\mathbf{A}\right)$ is given by $\rho \left(\mathbf{E},\mathbf{A}\right)=\mathbb{C}\backslash \sigma \left(\mathbf{E},\mathbf{A}\right)$.

The transfer function matrix $\mathbf{H}\left({\mathcal{S}}_n\right)\left(s\right)$ of the system ${\mathcal{S}}_n$ is:

$\begin{array}{cc}\hfill \mathbf{H}\left({\mathcal{S}}_n\right)\left(s\right)& =\mathbf{C}{\left(s\mathbf{E}-\mathbf{A}\right)}^{-1}\mathbf{B}+\mathbf{D}.\hfill \end{array}$

For the sake of clarity, the transfer function $\mathbf{H}\left({\mathcal{S}}_n\right)\left(s\right)$ for the system ${\mathcal{S}}_n$ is also denoted as $\mathbf{H}\left(s\right)$. The transfer function is said to be proper if $\underset{s\to \infty }{\lim }\mathbf{H}\left(s\right)$ is bounded (<∞ for each of its components) and strictly proper if $\underset{s\to \infty }{\lim }\mathbf{H}\left(s\right)=\mathbf{0}$. If the DAE system is split into the fast ${\mathbf{H}}_{\infty }\left(s\right)$ (including only the poles with infinite value) and slow subsystems ${\mathbf{H}}_p\left(s\right)$ (including only the poles with finite value) and restricted system equivalence relations are applied so that $\mathbf{H}\left(s\right)={\mathbf{H}}_p\left(s\right)+{\mathbf{H}}_{\infty }\left(s\right)$, the non-proper part of the transfer function matrix ${\mathbf{H}}_{\infty }\left(s\right)$ can be written as a Neumann series expansion [28]:

(2)$\begin{array}{cc}\hfill {\mathbf{H}}_{\infty }\left(s\right)& =\mathbf{C}{\left(s\mathbf{N}-\mathbf{I}\right)}^{-1}\mathbf{B}+\mathbf{D}=\mathbf{D}-\mathbf{C}\mathbf{B}-\sum _{j=1}^{\upsilon -1}{s}^j\mathbf{C}{\mathbf{N}}^j\mathbf{B},\hfill \end{array}$

where $\mathbf{N}$ is a nilpotent matrix with nilpotence index $\upsilon$. This provides a direct relation between the also called DAE index $\upsilon$ and the non-proper part of the transfer function matrix ${\mathbf{H}}_{\infty }\left(s\right)$, with $\nu$ one order less than that of the highest polynomial term.

The set of all possible minimal realizations of $\mathbf{H}\left({\mathcal{S}}_n\right)$ of size n is denoted by ${\mathbb{S}}_n$ such that ${\mathcal{S}}_n=\left(\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)\in {\mathbb{S}}_n$. The minimal realization of size n used here is based on the rational interpolation for tangential data [25] and the description of the tangential interpolation problem is based on that provided by Lefteriu et al. [29]. Tangential interpolation is a form of rational interpolation where the data is interpolated along particular directions. In particular, the data consist of the right interpolation data:

$\begin{array}{c}\hfill \left\{{\lambda }_i,{\mathbf{r}}_i,{\mathbf{w}}_i\mid {\lambda }_i\in \mathbb{C},{\mathbf{r}}_i\in {\mathbb{C}}^{n_y\times 1},{\mathbf{w}}_i\in {\mathbb{C}}^{n_u\times 1}\right\},\end{array}$

for $i=1,\dots ,n_r$:

$\begin{array}{c}\mathbf{\Lambda }=\mathrm{diag}\left[{\lambda }_1,\dots ,{\lambda }_{n_r}\right]\in {\mathbb{C}}^{n_r\times n_r},\hfill \\ \mathbf{R}=\left[{\mathbf{r}}_1,\dots ,{\mathbf{r}}_{n_r}\right]\in {\mathbb{C}}^{n_u\times n_r},\hfill \\ \mathbf{W}=\left[{\mathbf{w}}_1,\dots ,{\mathbf{w}}_{n_r}\right]\in {\mathbb{C}}^{n_y\times n_r},\hfill \end{array}$

and the left interpolation data:

$\begin{array}{c}\hfill \left\{{\mu }_j,{\mathbf{l}}_j,{\mathbf{v}}_j\mid {\mu }_j\in \mathbb{C},{\mathbf{l}}_j\in {\mathbb{C}}^{1\times n_u},{\mathbf{v}}_j\in {\mathbb{C}}^{1\times n_y}\right\},\end{array}$

for $j=1,\dots ,n_l$:

$\begin{array}{c}\mathbf{M}=\mathrm{diag}\left[{\mu }_1,\dots ,{\mu }_{n_l}\right]\in {\mathbb{C}}^{n_l\times n_l},\hfill \\ \mathbf{L}={\left[{\mathbf{l}}_1^T,\dots ,{\mathbf{l}}_{n_l}^T\right]}^T\in {\mathbb{C}}^{n_l\times n_y},\hfill \\ \mathbf{V}={\left[{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_{n_l}^T\right]}^T\in {\mathbb{C}}^{n_l\times n_u}.\hfill \end{array}$

The transfer function matrix is evaluated at the values ${\lambda }_i,{\mu }_j$ and the vectors ${\mathbf{r}}_i,{\mathbf{l}}_j$ are referred to as right and left tangential directions, while ${\mathbf{w}}_i,{\mathbf{v}}_j$ are the right and left tangential data. The rational interpolation problem for tangential data aims at finding a realization ${\mathcal{S}}_n=\left(\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)\in {\mathbb{S}}_n$ such that the associated transfer function $\mathbf{H}$ satisfies the right and left constraints,

(3)$\begin{array}{c}\mathbf{H}\left({\lambda }_i\right){\mathbf{r}}_i={\mathbf{w}}_i,i=1,\dots ,n_r,\hfill \\ {\mathbf{l}}_j\mathbf{H}\left({\mu }_j\right)={\mathbf{v}}_j,j=1,\dots ,n_l,\hfill \end{array}$

which is achieved by means of the Loewner and shifted Loewner matrices. Next, the structure of these two matrices is described.

A set Z of points in the complex plane $Z=\left\{z_1,\dots z_{n_r+n_l}\right\}$ and the corresponding values of the transfer function matrix is partitioned into the left and right data:

(4)$\begin{array}{cc}\hfill Z& =\left\{{\lambda }_1,\dots ,{\lambda }_{n_r}\right\}\cup \left\{{\mu }_1,\dots ,{\mu }_{n_l}\right\},\hfill \end{array}$

where the total number of sample points is $n_r+n_l$. The Loewner matrix is built as:

$\begin{array}{cc}\hfill \mathbb{L}& =\left[\begin{array}{ccc}\frac{{\mathbf{v}}_1{\mathbf{r}}_1-{\mathbf{l}}_1{\mathbf{w}}_1}{{\mu }_1-{\lambda }_1}& \cdots & \frac{{\mathbf{v}}_1{\mathbf{r}}_{n_r}-{\mathbf{l}}_1{\mathbf{w}}_{n_r}}{{\mu }_1-{\lambda }_{n_r}}\\ ⋮& \ddots & ⋮\\ \frac{{\mathbf{v}}_{n_l}{\mathbf{r}}_1-{\mathbf{l}}_{n_l}{\mathbf{w}}_1}{{\mu }_{n_l}-{\lambda }_1}& \cdots & \frac{{\mathbf{v}}_{n_l}{\mathbf{r}}_{n_r}-{\mathbf{l}}_p{\mathbf{w}}_{n_r}}{{\mu }_{n_l}-{\lambda }_{n_r}}\end{array}\right],\hfill \end{array}$

and satisfies the Sylvester equation:

$\begin{array}{cc}\hfill \mathbb{L}\mathbf{\Lambda }-\mathbf{M}\mathbb{L}& =\mathbf{L}\mathbf{W}-\mathbf{V}\mathbf{R}.\hfill \end{array}$

The Loewner matrix can be also expressed in terms of the tangential controllability and observability matrices [29]. If the directions ${\mathbf{r}}_i,{\mathbf{l}}_j$ are selected generically (random in practice), the rank of the Loewner matrix $\mathbb{L}$ is equal to the rank of the underlying matrix $\mathbf{E}$. The shifted Loewner matrix is the Loewner matrix corresponding to $s\mathbf{H}\left(s\right)$ and is built as:

$\begin{array}{cc}\hfill {\mathbb{L}}_{\sigma }& =\left[\begin{array}{ccc}\frac{{\mu }_1{\mathbf{v}}_1{\mathbf{r}}_1-{\lambda }_1{\mathbf{l}}_1{\mathbf{w}}_1}{{\mu }_1-{\lambda }_1}& \cdots & \frac{{\mu }_1{\mathbf{v}}_1{\mathbf{r}}_{n_r}-{\lambda }_{n_r}{\mathbf{l}}_1{\mathbf{w}}_{n_r}}{{\mu }_1-{\lambda }_{n_r}}\\ ⋮& \ddots & ⋮\\ \frac{{\mu }_{n_l}{\mathbf{v}}_{n_l}{\mathbf{r}}_1-{\lambda }_1{\mathbf{l}}_{n_l}{\mathbf{w}}_1}{{\mu }_{n_l}-{\lambda }_1}& \cdots & \frac{{\mu }_{n_l}{\mathbf{v}}_{n_l}{\mathbf{r}}_{n_r}-{\lambda }_{n_r}{\mathbf{l}}_{n_l}{\mathbf{w}}_{n_r}}{{\mu }_{n_l}-{\lambda }_{n_r}}\end{array}\right],\hfill \end{array}$

and satisfies the Sylvester equation:

$\begin{array}{cc}\hfill {\mathbb{L}}_{\sigma }\mathbf{\Lambda }-\mathbf{M}{\mathbb{L}}_{\sigma }& =\mathbf{L}\mathbf{W}\mathbf{\Lambda }-\mathbf{M}\mathbf{V}\mathbf{R}.\hfill \end{array}$

Assuming that $n_r=n_l$ and that $\det \left({\mathbb{L}}_{\sigma }-s\mathbb{L}\right)\ne 0$ a minimal realization ${\mathcal{S}}_n=\left(\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)\in {\mathbb{S}}_n$ for the linear descriptor system is given by [25]:

(5)$\begin{array}{cc}\hfill \mathbf{E}& =-\mathbb{L},\mathbf{A}=-{\mathbb{L}}_{\sigma },\mathbf{B}=\mathbf{V},\mathbf{C}=\mathbf{W},\mathbf{D}=\mathbf{0},\hfill \end{array}$

and the associated transfer function $\mathbf{H}\left(s\right)=\mathbf{W}{\left({\mathbb{L}}_{\sigma }-s\mathbb{L}\right)}^{-1}\mathbf{V}$ satisfies the right and left interpolation conditions of Equation (3). An additional SVD factorization of the Loewner matrix $\mathbb{L}$ allows for the further size reduction of the system, generating a reduced order model (ROM) of the descriptor system:

(6)$\begin{array}{cc}\hfill \mathbb{L}& =\left[\begin{array}{cc}{\mathbf{Y}}_1& {\mathbf{Y}}_2\end{array}\right]\left[\begin{array}{cc}{\mathit{\Sigma }}_1& \mathbf{0}\\ \mathbf{0}& {\mathit{\Sigma }}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{X}}_1^{*}\\ {\mathbf{X}}_2^{*}\end{array}\right],\hfill \end{array}$

where * denotes the conjugate transpose, ${\mathit{\Sigma }}_1\in {\mathbb{R}}^{r\times r}$, ${\mathit{\Sigma }}_2\in {\mathbb{R}}^{\left(n-r\right)\times \left(n-r\right)}$ and ${\mathbf{X}}_1,{\mathbf{X}}_2,{\mathbf{Y}}_1,{\mathbf{Y}}_2$ are of appropriate dimensions. The reduced system of size r defined by $S_r=\left({\mathbf{E}}_r,{\mathbf{A}}_r,{\mathbf{B}}_r,{\mathbf{C}}_r,{\mathbf{D}}_r\right)\in {\mathbb{S}}_r$ is obtained by a Petrov–Galerkin projection and represents the best approximation to the full Loewner matrix in the Frobenius or 2-norm:

(7)$\begin{array}{cc}\hfill {\mathbf{E}}_r& =-{\mathbf{Y}}_1^{*}\mathbb{L}{\mathbf{X}}_1,{\mathbf{A}}_r=-{\mathbf{X}}_1^{*}{\mathbb{L}}_{\sigma }{\mathbf{X}}_1,{\mathbf{B}}_r={\mathbf{Y}}_1^{*}\mathbf{V},{\mathbf{C}}_r=\mathbf{W}{\mathbf{X}}_1,{\mathbf{D}}_r=\mathbf{0}.\hfill \end{array}$

As a remark, note that the rational interpolation problem is different from the classical one of rational approximation. In the framework of rational interpolation presented here, the generalized system obtained by the Loewner realization exactly interpolates the transfer function matrix $\mathbf{H}$ given at the sample points in tangential directions, satisfying the conditions given in Equation (3). The rational approximation instead provides a least-squares fit to the data and thus a systematic error is introduced. As a result, some authors impose additional constraints (usually at zero or flutter frequencies) to improve the fit [13]. In order not to deteriorate the fit at other frequencies, the number of terms for the least-squares fit must then be increased.

Once the system $S_r=\left({\mathbf{E}}_r,{\mathbf{A}}_r,{\mathbf{B}}_r,{\mathbf{C}}_r,{\mathbf{D}}_r\right)\in {\mathbb{S}}_r$ has been determined, it may contain unstable poles. Assuming that the aeroservoelastic system under consideration is stable, the system is then approximated by a stable one. This is done by means of an approximation in the space of rational functions contained in the Hardy space $H_{\infty }$ as described in Section 2.2.

2.2. Optimal $R{H}_{\infty }$ Approximation

The following classes of descriptor systems are defined:

$\begin{array}{c}{\mathbb{S}}_0=\left\{{\mathbf{E}}_r,{\mathbf{A}}_r,{\mathbf{B}}_r,{\mathbf{C}}_r,{\mathbf{D}}_r\in {\mathbb{S}}_r\mid \rho \left({\mathbf{E}}_r,{\mathbf{A}}_r\right)\cap i\mathbb{R}=\varphi \right\},\hfill \\ {\mathbb{S}}_{+}=\left\{{\mathbf{E}}_{+},{\mathbf{A}}_{+},{\mathbf{B}}_{+},{\mathbf{C}}_{+},{\mathbf{D}}_{+}\in {\mathbb{S}}_r\mid \rho \left({\mathbf{E}}_{+},{\mathbf{A}}_{+}\right)\subset {\mathbb{C}}_{+}\right\},\hfill \\ {\mathbb{S}}_{-}=\left\{{\mathbf{E}}_{-},{\mathbf{A}}_{-},{\mathbf{B}}_{-},{\mathbf{C}}_{-},{\mathbf{D}}_{-}\in {\mathbb{S}}_r\mid \rho \left({\mathbf{E}}_{-},{\mathbf{A}}_{-}\right)\subset {\mathbb{C}}_{-}\right\},\hfill \end{array}$

where $\varphi$ is the empty set, $\rho \left(\mathbf{E},\mathbf{A}\right)$ the resolvent set as described in Section 2.1, ${\mathbb{S}}_r$ represents all realizations of size r as in Equation (7), ${\mathbb{S}}_{+}$ represents the class of stable descriptor systems and ${\mathbb{S}}_{-}$ represents the class of unstable descriptor systems, both of size r. The system obtained from tangential interpolation as specified in Equation (7) may contain artificial unstable poles, in the sense that the original system is stable, as the method described in Section 2.1 does not preserve stability. Assuming that the reduced system $S_r$ obtained by Equation (7) does not contain any poles in the imaginary axis ($S_r\subset {\mathbb{S}}_0$), the problem of finding the best approximation of the system $S_r$ by a stable system ${\widehat{S}}_{r+}\subset {\mathbb{S}}_{+}$ minimizes the distance between the transfer function matrices $\mathbf{H}$:

$\begin{array}{cc}\hfill {\widehat{S}}_{r+}& =\underset{S\subset {\mathbb{S}}_{+}}{\arg \min }{∥\mathbf{H}\left(S_r\right)-\mathbf{H}\left(S\right)∥}_p.\hfill \end{array}$

This continuous-time problem has been explicitly solved by Köhler [30], requiring additionally that the matrix ${\mathbf{E}}_{+}$ for the class of stable descriptor systems $S\subset {\mathbb{S}}_{+}$ be invertible, that is, with a proper transfer function matrix such that $\underset{s\to \infty }{\lim }\mathbf{H}\left(S\right)\left(s\right)$ is bounded (<∞ for each of its components). The subindex $p\in \left\{2,\infty \right\}$ represents the norm to be minimized in the corresponding Hardy space $H_p$. The space $R{H}_{\infty }$ includes all stable bounded rational transfer function matrices $\mathbf{H}$ in the Hardy space $H_{\infty }$ ($p=\infty$) of analytic functions in the (open) right half of the complex plane ${\mathbb{C}}_{+}$:

$\begin{array}{cc}\hfill H_{\infty }& =\left\{\mathbf{H}:\mathbb{C}\to {\mathbb{C}}^{n_y\times n_u}\mid {∥\mathbf{H}∥}_{\infty }<\infty \right\},{∥\mathbf{H}∥}_{\infty }=\underset{s\in {\mathbb{C}}_{+}}{\sup }{∥\mathbf{H}\left(s\right)∥}_{\infty }=\underset{\omega \in \mathbb{R}}{\sup }{\sigma }_{max}\left(\mathbf{H}\left(i\omega \right)\right),\hfill \end{array}$

with ${\sigma }_{max}$ the maximum singular value of the transfer function matrix $\mathbf{H}\left(i\omega \right)$.

The optimal solution in the $R{H}_{\infty }$ space is given by Köhler [30]:

(8)$\begin{array}{cc}\hfill {\widehat{S}}_{r+}& =S_{r+}\oplus \mathbf{P}{S}_{r-,{\sigma }_1}\mathbf{Q},\hfill \end{array}$

where the system $S_{r+}$ corresponds to the original system $S_r$ but keeping solely the stable poles (with negative real part). The system $S_{r-,{\sigma }_1}\subset {\mathbb{S}}_{-}$ is $S_{r-,{\sigma }_1}=\left({\mathbf{E}}_{r-}^T{\mathbf{R}}_1,-{\mathbf{A}}_{r-}^T{\mathbf{R}}_1-{\mathbf{C}}_{r-}^T{\mathbf{C}}_1,{\mathbf{E}}_{r-}^T\mathbf{Q}{\mathbf{B}}_{r-},{\mathbf{C}}_{r-}\mathbf{P}{\mathbf{E}}_{r-}^T,{\mathbf{D}}_{r-}\right)$ with:

$\begin{array}{c}{\mathbf{C}}_1={\mathbf{C}}_{r-}\mathbf{P}{\mathbf{E}}_{r-}^T,\hfill \\ {\mathbf{R}}_1=\mathbf{Q}{\mathbf{E}}_{r-}\mathbf{P}{\mathbf{E}}_r^T-{\sigma }_1^2\mathbf{I},\hfill \end{array}$

and $\mathbf{P}$ and $\mathbf{Q}$ are the unique infinite gramians that satisfy the following generalized Lyapunov equations:

$\begin{array}{c}{\mathbf{A}}_{r-}\mathbf{P}{\mathbf{E}}_{r-}^T+{\mathbf{E}}_{r-}^T\mathbf{P}{\mathbf{A}}_{r-}^T+{\mathbf{B}}_{r-}{\mathbf{B}}_{r-}^T=\mathbf{0},\hfill \\ {\mathbf{A}}_{r-}\mathbf{Q}{\mathbf{E}}_{r-}+{\mathbf{E}}_{r-}^T\mathbf{Q}{\mathbf{A}}_{r-}^T+{\mathbf{C}}_{r-}^T{\mathbf{C}}_{r-}^T=\mathbf{0}.\hfill \end{array}$

The optimal solution ${\widehat{S}}_{r+}$ in the $R{H}_{\infty }$ space verifies that ${∥\mathbf{H}\left(S_r\right)-\mathbf{H}({\widehat{S}}_{r+})∥}_{\infty }={\sigma }_1$, where ${\sigma }_1={∥S_{r-}∥}_{\infty }=\sqrt{\max \left(\mathbf{P}\mathbf{Q}\right)}$ is the greatest Hankel singular value of the system $S_{r-}\subset {\mathbb{S}}_{-}$ which results from the original system $S_r$ but keeping solely the unstable poles (with positive real part).

2.3. Application to Unsteady Incompressible Flow

As it is clear from Section 1, the main task for the application of the FSM method for the cut loads recovery in an aeroservoelastic framework remains to properly describe the distributed unsteady aerodynamic forces in the time domain. In order to show the convenience of the tangential interpolation method described in Section 2.1, an application for the incompressible unsteady flow over a bidimensional airfoil with heave and pitch degrees of freedom is presented, showing the connection between the DAE theory or descriptor systems and the unsteady aerodynamic model.

The classical unsteady model of Theodorsen [31] provides the lift and pitch moment coefficient (excluding the added-mass terms) to arbitrary input motions of the airfoil. The Wagner model in the time domain is equivalent to that of Theodorsen in the frequency domain, which uses the physical assumptions of inviscid flow, incompressibility and planar wake. As described by Brunton [32], several authors have constructed state-space models corresponding to the Theodorsen formulation. Note that excluding the models obtained by Peters [33], which requires eight states for an appropriate system representation, these state-space formulations do not include the added-mass terms. In this section, a time-domain model including also the added-mass terms is obtained.

The local lift $c_l$ (positive upwards) and local pitch moment at the quarter chord $c_m$ (positive nose up) can be described as a function of the heave $u_z$ (positive down) and pitch $u_{\theta }$ (positive nose up) motion in the frequency domain by means of an aerodynamic transfer function matrix:

$\begin{array}{cc}\hfill \left[\begin{array}{c}{c}_l\left(\omega \right)\\ c_m\left(\omega \right)\end{array}\right]& =\mathbf{H}\left(\omega \right)\left[\begin{array}{c}{u}_z\left(\omega \right)\\ u_{\theta }\left(\omega \right)\end{array}\right],\hfill \end{array}$

where the system output is $\mathbf{y}={\left[c_l{c}_m\right]}^T$, the system input is $\mathbf{u}={\left[u_z{u}_{\theta }\right]}^T$ and the transfer function matrix is [34]:

(9)$\begin{array}{cc}\hfill \mathbf{H}\left(\omega \right)& =\left[\begin{array}{cc}-{\omega }^2\frac{\pi L_{ref}}{U_{\infty }^2}+i\omega \frac{2\pi }{U_{\infty }}C\left(k\right)& {\omega }^2\frac{\pi a{L}_{ref}^2}{U_{\infty }^2}+i\omega \frac{\pi L_{ref}}{U_{\infty }}\left(1+\left(1-2a\right)C\left(k\right)\right)+2\pi C\left(k\right)\\ {\omega }^2\frac{\pi L_{ref}}{4U_{\infty }^2}& {\omega }^2\left(\frac{1}{8}-\frac{a}{2}\right)\frac{\pi L_{ref}^2}{2U_{\infty }^2}-i\omega \frac{\pi L_{ref}}{2U_{\infty }}\end{array}\right],\hfill \end{array}$

where k is the reduced frequency $k=\omega L_{ref}/U_{\infty }$ with $L_{ref}$ a reference length (the half-chord) and $U_{\infty }$ the true airspeed. The parameter a represents the pitch axis location with respect to the half-chord and ranges between $-1$ for the leading edge and 1 for the trailing edge. The Theodorsen function $C\left(k\right)$ is given by:

$\begin{array}{cc}\hfill C\left(k\right)& =\frac{H_1^{\left(2\right)}\left(k\right)}{H_1^{\left(2\right)}\left(k\right)+i{H}_0^{\left(2\right)}\left(k\right)},\hfill \end{array}$

with $H_0$ and $H_1$ the Hankel functions of first and second kind [31]. Interestingly, the pitch moment at the quarter chord position is independent of the Theodorsen function. Note that, in this work, all time-dependent dynamic variables such as $c_l$, $c_m$, $u_z$ and $u_{\theta }$ are incremental variables in the sense that they refer to deviations with respect to a steady-state reference.

It is clear from Equation (9) that the aerodynamic system is not proper, as the values $\underset{s\to \infty }{\lim }\mathbf{H}\left(s\right)$ are not bounded. According to Equation (2), the corresponding time-domain description must be in the form of a descriptor system of index $\nu -1=2$, that is, of index $\upsilon =3$. The generalized realization based on tangential interpolation of Section 2.1 can then be applied. However, it is well known that the numerical time-domain solution of DAE systems becomes very difficult for index values higher than 2 [35] and thus a procedure to reduce the DAE index for this particular aerodynamic system is proposed. With the transfer function matrix given by Equation (9), the DAE index can be reduced by simply reducing the polynomial order at $s\to \infty$. Even though the coefficients corresponding to the non-proper polynomial terms in Equation (9) are explicitly known, a general method based on the evaluation of the polynomial coefficients by means of the derivatives of the transfer function matrix $\mathbf{H}\left(\omega \right)$ with respect to the circular frequency $\omega$ is given. This results in a more general setting, as typically the aerodynamic system description is provided at a finite set of frequency values, without the explicit knowledge of the transfer function matrix. The resulting transfer matrix function ${\mathbf{H}}_2\left(\omega \right)$ which excludes the second order polynomial part can be obtained as:

$\begin{array}{cc}\hfill {\mathbf{H}}_2\left(\omega \right)& =\mathbf{H}\left(\omega \right)-\frac{{\omega }^2}{2}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }\frac{d^2\mathbf{H}\left(\omega \right)}{d{\omega }^2}\right\}.\hfill \end{array}$

For the present case, the derivatives of $\mathbf{H}\left(\omega \right)$ with respect to $\omega$ have been computed by (central) finite differences and the maximum reduced frequency value considered is $k_{max}=10$. Note that the high frequency values required for the proper estimation of the limit behaviour when $\omega \to \infty$ may involve frequencies beyond the range of interest for the application under consideration. Now, the tangential interpolation method can be applied to ${\mathbf{H}}_2\left(\omega \right)$ and a realization $\left({\mathbf{E}}_{2r},{\mathbf{A}}_{2r},{\mathbf{B}}_{2r},{\mathbf{C}}_{2r},{\mathbf{D}}_{2r}\right)$ as given by Equation (7) is obtained. For the tangential interpolation method, the set Z of Equation (4) is obtained as $Z=\left\{{\lambda }_1,\dots ,{\lambda }_{n_r}\right\}\cup \left\{{\mu }_1,\dots ,{\mu }_{n_l}\right\}$, where ${\lambda }_i=i{\omega }_i$ (the running index i and the complex number i should not be confused) and ${\mu }_j=i{\omega }_j$. Evaluating the transfer function matrix at this set of frequencies provides $\mathbf{H}\left({\lambda }_i\right)$ and $\mathbf{H}\left({\mu }_j\right)$. The tangential directions ${\mathbf{r}}_i,{\mathbf{l}}_j$ have been selected randomly, as the realization does not depend on them. In order to take into account the term corresponding to the coefficient of the ${\omega }^2$ term in $\mathbf{H}\left(\omega \right)$, the input vector is expanded including the second derivative terms as well, $\ddot{\mathbf{u}}={\left[{\ddot{u}}_z{\ddot{u}}_{\theta }\right]}^T$. Thus, the matrices ${\mathbf{B}}_2$ and ${\mathbf{D}}_2$ as obtained by Equation (7) are further modified to accommodate the second derivative term $\ddot{\mathbf{u}}$ and the resulting DAE system of reduced index 2 is:

(10)$\begin{array}{c}{\mathbf{E}}_{2r}{\dot{\mathbf{x}}}_2={\mathbf{A}}_{2r}{\mathbf{x}}_2+\left[\begin{array}{cc}{\mathbf{B}}_{2r}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \ddot{\mathbf{u}}\end{array}\right],\hfill \\ \left[\begin{array}{c}{c}_l\\ c_m\end{array}\right]={\mathbf{C}}_{2r}{\mathbf{x}}_2+\left[\begin{array}{cc}\mathbf{0}& -\frac{1}{2}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }\frac{d^2\mathbf{H}\left(\omega \right)}{d{\omega }^2}\right\}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \ddot{\mathbf{u}}\end{array}\right],\hfill \end{array}$

where the number of states or size of the matrix ${\mathbf{A}}_{2r}$ can be controlled by the number of singular values retained in the SVD decomposition given in Equation (6).

In a natural way, the same procedure can be applied to further reduce the system to a DAE of index 0 or set of ordinary differential equations (ODE). In that case, the transfer function is required to be strictly proper:

(11)$\begin{array}{cc}\hfill {\mathbf{H}}_0\left(\omega \right)& =\mathbf{H}\left(\omega \right)-\frac{1}{2}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }\frac{d^2\mathbf{H}\left(\omega \right)}{d{\omega }^2}\right\}-\mathrm{Im}\left\{\underset{\omega \to \infty }{\lim }\frac{d{\mathbf{H}}_2\left(\omega \right)}{d\omega }\right\}-\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }{\mathbf{H}}_1\left(\omega \right)\right\},\hfill \end{array}$

where:

$\begin{array}{cc}\hfill {\mathbf{H}}_1\left(\omega \right)& ={\mathbf{H}}_2\left(\omega \right)-i\omega \mathrm{Im}\left\{\underset{\omega \to \infty }{\lim }\frac{d{\mathbf{H}}_2\left(\omega \right)}{d\omega }\right\}.\hfill \end{array}$

Applying the tangential interpolation technique to the strictly proper transfer function ${\mathbf{H}}_0\left(\omega \right),$ the set of matrices $\left(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)$ representing an ODE system is obtained:

$\mathbf{A}={\mathbf{E}}_{0r}^{-1}{\mathbf{A}}_{0r},\mathbf{B}={\mathbf{E}}_{0r}^{-1}{\mathbf{B}}_{0r},\mathbf{C}={\mathbf{E}}_{0r}^{-1}{\mathbf{C}}_{0r},\mathbf{D}={\mathbf{E}}_{0r}^{-1}{\mathbf{D}}_{0r}=\mathbf{0}.$

In case any unstable poles exist in this system, a projection into the $R{H}_{\infty }$ space is done by applying Equation (8).

Additionally, the system input vector is expanded to include the first and second time derivatives $\dot{\mathbf{u}}$ and $\ddot{\mathbf{u}}$ and the matrices $\mathbf{B}$ and $\mathbf{D}$ are accordingly modified to account for the terms substracted in Equation (11):

(12)$\begin{array}{c}{\dot{\mathbf{x}}}_0=\mathbf{A}{\mathbf{x}}_0+\left[\begin{array}{ccc}{\mathbf{B}}_1& \mathbf{0}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \dot{\mathbf{u}}\\ \ddot{\mathbf{u}}\end{array}\right],\hfill \\ \left[\begin{array}{c}{c}_l\\ c_m\end{array}\right]=\mathbf{C}{\mathbf{x}}_0+\left[\begin{array}{ccc}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }{\mathbf{H}}_1\left(\omega \right)\right\}& \mathrm{Im}\left\{\underset{\omega \to \infty }{\lim }\frac{d{\mathbf{H}}_2\left(\omega \right)}{d\omega }\right\}& -\frac{1}{2}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }\frac{d^2\mathbf{H}\left(\omega \right)}{d{\omega }^2}\right\}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \dot{\mathbf{u}}\\ \ddot{\mathbf{u}}\end{array}\right],\hfill \end{array}$

where $\mathbf{B}=\left[{\mathbf{B}}_1\mathbf{0}\mathbf{0}\right]$ and $\mathbf{D}=[\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }{\mathbf{H}}_1\left(\omega \right)\right\}\mathrm{Im}\left\{\underset{\omega \to \infty }{\lim }\frac{d{\mathbf{H}}_2\left(\omega \right)}{d\omega }\right\}-\frac{1}{2}\mathrm{Re}\left\{\underset{\omega \to \infty }{\lim }\frac{d^2\mathbf{H}\left(\omega \right)}{d{\omega }^2}\right\}]$. Note that even though the system described by Equation (12) can be solved as an ODE, it is not an ODE in the classical sense as it includes the input time derivatives, which (opposite to the input time integrals) cannot be represented in a regular state-space system.

In order to show the suitability of this method, both the DAE of index 2 (DAE2) and ODE systems described by Equations (10) and (12) respectively have been solved numerically for the input:

$\begin{array}{c}\mathbf{u}\left(t\right)={\left[\begin{array}{cc}4u_{z0}{\left(t/t_0\right)}^2{e}^{\left(2-1/\left(1-t/t_0\right)\right)}& 4u_{\theta 0}{\left(t/t_0\right)}^2{e}^{\left(2-1/\left(1-t/t_0\right)\right)}\end{array}\right]}^T,0\le t\le t_0,\hfill \\ \mathbf{u}\left(t\right)={\left[\begin{array}{cc}0& 0\end{array}\right]}^T,t>t_0,\hfill \end{array}$

with $u_z$ in (m), $u_{\theta }$ in (rad), $u_{z0}=0.01$ (m) and $u_{\theta 0}=\pi /180$ (rad) as shown in Figure 1. The reference length is $L_{ref}=0.5$ (m) and the elastic axis is slightly behind the center of gravity as given by the parameter $a=0.1$. The respective solutions are compared against the reference solution obtained in the frequency domain and labeled as FD in Figure 2. For the numerical solution of the DAE of index 2, a solver [36] based on the algorithms described by Petzold [37] has been used and a total number of three states has been kept for the ODE representation. The time-domain results of the frequency domain (FD) solution correspond to the inverse Fourier transform of the output vector. Note the excellent agreement between the three solutions. The solution corresponding to the ODE requires less computational time and is the preferable option when computing the output response with a generalized state-space model in the time domain.

Note that, in the case of incompressible flow, the apparent mass is represented in the limit of $s\to \infty$ as a polynomial in s of order 2. For the case of compressible flow, the aerodynamic forces due to a sudden change in the structural motion are no longer impulsive but finite and time dependent [38]. Miles [39] indicated that the concept of added-mass is applicable in compressible flow only for the steady state case, but not for the unsteady compressible case. Instead of referring to the added-mass concept, the behaviour of the aerodynamic forces in the limit $s\to \infty$ can be considered. This limit has been analyzed by Vepa [40] for the more general doublet lattice method (DLM). He showed that the DLM method does not converge uniformly for all frequencies due to the lattice integration scheme. In particular, the DLM method does not converge to the results predicted by piston theory for increasing reduced frequencies and for certain mode shapes. For convergence, the distance between the sending and receiving point on each panel must vanish according to the piston theory. This would require, theoretically, an infinite number of boxes in the chordwise direction. He also showed that, for three-dimensional subsonic flow, there exists a Padé approximant sequence with $\underset{i,j}{\max }\left(deg\left(n_{ij}\right)-deg\left(d_{ij}\right)\right)=1$ (where $n_{ij}$ and $d_{ij}$ are respectively the components of the numerator and denominator of the aerodynamic transfer function matrix), which converges to the exact aerodynamic loads for an infinite number of terms. Thus, the behaviour of the subsonic aerodynamic loads may be described by a polynomial term in s in the limit $s\to \infty$. For the incompressible case, the highest order term in the denominator becomes zero and $\underset{i,j}{\max }\left(deg\left(n_{ij}\right)-deg\left(d_{ij}\right)\right)=2$, requiring a second order polynomial expansion in s to describe the behaviour in the limit $s\to \infty$, see Equation (9).

However, the dependence of the unsteady aerodynamic forces on the second time derivative resulting from the classical RFA fit [4] for compressible flow may be seen as a low frequency residualization of its high frequency behaviour [41]. In this work, a residualization of the aerodynamic forces is also done in the limit of $s\to \infty$ when considering compressible flow, see Section 3.1.

Application to the Theodorsen Function

In the case of bidimensional unsteady incompressible flow, a more general realization can be obtained which is then valid for all values of the parameters involved in Equation (9), namely $U_{\infty }$, $L_{ref}$ and a. This is not possible for the more general three-dimensional unsteady compressible flow and an extension of the method presented in Section 2.3 involving a residualization of the aerodynamic matrix transfer function at high frequency values will be applied in Section 3.1.

For a realization of the Theodorsen function, the following strictly proper aerodynamic transfer function is defined:

$\begin{array}{cc}\hfill H_c\left(k\right)& =C\left(k\right)-\frac{1}{2},\hfill \end{array}$

where $C\left(k\right)$ represents the Theodorsen function which depends on the reduced frequency k. Next, the Loewner and shifted Loewner matrices are built and a realization is carried out:

$\begin{array}{cc}\hfill {\mathbf{A}}_c& ={\mathbf{E}}_{cr}^{-1}{\mathbf{A}}_{cr},{\mathbf{B}}_c={\mathbf{E}}_{cr}^{-1}{\mathbf{B}}_{cr},{\mathbf{C}}_c={\mathbf{E}}_{cr}^{-1}{\mathbf{C}}_{cr},{\mathbf{D}}_c={\mathbf{E}}_{cr}^{-1}{\mathbf{D}}_{cr}=\mathbf{0},\hfill \end{array}$

where the system $\left({\mathbf{E}}_{cr},{\mathbf{A}}_{cr},{\mathbf{B}}_{cr},{\mathbf{C}}_{cr},{\mathbf{D}}_{cr}\right)$ is obtained after applying Equations (5) and (7) for a reduction on the number of states and Equation (8) for the suppression of possible unstable poles. The resulting generalized state-space system is:

(13)$\begin{array}{cc}\hfill {\dot{\mathbf{x}}}_c& ={\mathbf{A}}_c{\mathbf{x}}_c+{\mathbf{B}}_c\frac{2\pi }{U_{\infty }}\left[\begin{array}{cccccc}0& U_{\infty }& 1& L_{ref}\left(\frac{1}{2}-a\right)& 0& 0\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \dot{\mathbf{u}}\\ \ddot{\mathbf{u}}\end{array}\right]\hfill \\ \hfill \left[\begin{array}{c}{c}_l\\ c_m\end{array}\right]& =\left[\begin{array}{c}{\mathbf{C}}_c\\ \mathbf{0}\end{array}\right]{\mathbf{x}}_c+\left(\left[\begin{array}{cccccc}0& \pi & \frac{\pi }{U_{\infty }}& \frac{\pi L_{ref}}{U_{\infty }}\left(\frac{3}{2}-a\right)& \frac{\pi L_{ref}}{U_{\infty }^2}& -\frac{\pi a{L}_{ref}^2}{U_{\infty }^2}\\ 0& 0& 0& -\frac{\pi L_{ref}}{2U_{\infty }}& -\frac{\pi L_{ref}}{4U_{\infty }^2}& -\left(\frac{1}{8}-\frac{a}{2}\right)\frac{\pi L_{ref}^2}{2U_{\infty }^2}\end{array}\right]\right)\left[\begin{array}{c}\mathbf{u}\\ \dot{\mathbf{u}}\\ \ddot{\mathbf{u}}\end{array}\right].\hfill \end{array}$

Figure 3 shows different realizations of the Theodorsen function $C\left(k\right)$ in the complex plane for different values of the reduced frequency k. The Jones approximation requires two states. For the Loewner realization, the cases corresponding to a number of states equal to that of Jones ($n=2$) and eight are depicted. It becomes clear that the Loewner realization for an increasing number of poles is in agreement with the reference Theodorsen function computed directly in the frequency domain (FD). The number of poles in the Loewner realization is increased by simply modifying the truncation after the SVD decomposition of Equation (6). Note that, unlike for the classical interpolation techniques [42], the Loewner realization does not require the selection of any poles and no iterative optimization is involved. The matrices ${\mathbf{A}}_c$, ${\mathbf{B}}_c$, ${\mathbf{C}}_c$ and ${\mathbf{D}}_c$ obtained for the case where up to eight states are retained after the SVD decomposition of the Loewner matrix $\mathbb{L}$ are provided in Appendix A.

3. Aeroservoelastic System for Loads Analysis

In this section, the tangential interpolation method presented in Section 2.1 is applied to the more general case of an aircraft flying in subsonic compressible flow. Once the aerodynamic system has been obtained in generalized state-space form in Section 3.1, the structural model is also considered and a coupled aeroservoelastic system for dynamic loads prediction including both the MDM and the FSM methods is obtained in Section 3.2.

3.1. Generalized Realization of the Aerodynamic System

For the general case of three-dimensional unsteady compressible flow, the transfer function matrix describing the aerodynamic forces cannot be explicitly obtained as for the bidimensional incompressible case, given by Equation (9). A common technique for the computation of the subsonic compressible aerodynamic forces for dynamic loads applications is the doublet lattice method (DLM), which provides the matrix of aerodynamic influence coefficients (AIC) at a set of particular frequencies. Due to the assumption of linear and inviscid flow, the DLM method is not appropriate to describe the transonic flow. In this region, other methods such as the correction of the AIC matrices [43] or linearized computational fluid dynamics (CFD) solvers in the frequency domain can be used [44]. An application including CFD data considering a gust disturbance as input has been shown by Poussot-Vassal et al. [45].

The input of the aerodynamic system is given by the vertical atmospheric disturbance $w_g$ at a particular spatial location (the aircraft nose) and a set of generalized coordinates corresponding to the control surfaces ${\mathbf{u}}_c$ and the structural modes ${\mathbf{u}}_h$ of the aircraft structure (in vacuum) including the rigid-body and flexible modes, as the aeroservoelastic system is reduced by projection into the eigenvalues of the structure (modal truncation). Other directions of the atmospheric disturbance can be included by considering a linear combination of disturbances along three orthogonal axes. Within the usual assumptions for aeroservoelastic modeling, the translational degrees of freedom ${\mathbf{u}}_t$ of the aircraft structure do not cause any change in the aerodynamic forces and thus their time derivative is considered here as input to the aerodynamic system. The vector ${\mathbf{u}}_{rf}$ contains the generalized coordinates corresponding to rotational rigid-body motions and to the flexible modes. The aerodynamic forces due to the control surfaces are evaluated by means of a transpiration boundary condition by applying a local variation of the downwash velocity on the related aerodynamic panels but keeping the original panel geometry.

As output ${\mathbf{y}}_a$ of the aerodynamic system, the (incremental) local lift and local pitch moment coefficients for the aerodynamic strips defined over the lifting components ${\mathbf{y}}_a={\left[{\mathbf{c}}_l^T{\mathbf{c}}_m^T\right]}^T$ are chosen (see Section 4). With the aerodynamic distribution provided by these coefficients, the cut loads distribution can be recovered over the load reference axis (LRA) by means of the FSM. The aerodynamic transfer function matrix relates the input and output vectors in the frequency domain:

(14)$\begin{array}{cc}\hfill \left[\begin{array}{c}{\mathbf{c}}_l\left(\omega \right)\\ {\mathbf{c}}_m\left(\omega \right)\end{array}\right]& =\left[\begin{array}{cccc}{\mathbf{H}}_g\left(\omega \right)& {\mathbf{H}}_t\left(\omega \right)& {\mathbf{H}}_{rf}\left(\omega \right)& {\mathbf{H}}_c\left(\omega \right)\end{array}\right]\left[\begin{array}{c}{w}_g\left(\omega \right)\\ i\omega {\mathbf{u}}_t\left(\omega \right)\\ {\mathbf{u}}_{rf}\left(\omega \right)\\ {\mathbf{u}}_c\left(\omega \right)\end{array}\right],\hfill \end{array}$

where ${\mathbf{H}}_t\left(\omega \right))=\mathbf{H}\left(\omega \right){\mathbf{T}}_H/i\omega$ for $\omega >0$ with ${\mathbf{T}}_H$ a matrix selecting the columns from $\mathbf{H}\left(\omega \right)$ corresponding to the translational degrees of freedom ${\mathbf{u}}_t$. In this way, the singularity caused by the translational motion at frequency zero as described in Section 1 is dealt with systematically.

Within the DLM method, the integrals required to solve the AIC matrices have to be approximated [46]. When considering the aerodynamic transfer function matrix $\mathbf{H}\left(\omega \right)$ in the limit $\omega \to \infty$, the polynomial coefficients cannot be obtained from the derivatives of $\mathbf{H}\left(\omega \right)$ with respect to the circular frequency $\omega$ as done in the incompressible case (see Section 2.3). Figure 4 shows a typical element of the transfer function matrix $\mathbf{H}\left(\omega \right)$ computed numerically for increasing values of the reduced frequency k. Note that the high reduced frequency values shown are beyond the applications of practical interest but would allow for the numerical determination of the non-proper part. It is clear that, for a common DLM aerodynamic model discretization, the method is not able to handle very high frequencies. Indeed, a very high number of panels would be required to properly describe the behaviour at increasing reduced frequencies. Another approach would be to compute the values at very high reduced frequency by means of the piston theory [47]. In this work, a residualization of the transfer function matrix at high frequency values is done by representing it as the sum of a proper and a non-proper part instead.

First, a realization on the aerodynamic transfer function matrix as given in Equation (5) is done. Once the realization $\left(\widehat{\mathbf{E}},\widehat{\mathbf{A}},\widehat{\mathbf{B}},\widehat{\mathbf{C}},\widehat{\mathbf{D}}\right)$ has been obtained the eigenvalues of the system are obtained by solving $\det \left(s\widehat{\mathbf{E}}-\widehat{\mathbf{A}}\right)=0$. Note that, at this stage, there is no need for an SVD decomposition as this realization is not the final one representing the aerodynamic system. Similarly, no stabilization as described in Section 2.2 is enforced yet. The possible infinite poles due to the non-properness are excluded by setting a tolerance value $S_{\infty }$. The remaining finite poles are collected in the set $S_f$:

(15)$\begin{array}{cc}\hfill S_f& =\left\{s\in \sigma \left(\mathbf{E},\mathbf{A}\right)\mid \left|s\right|\le S_{\infty }\right\}.\hfill \end{array}$

The column of the aerodynamic transfer function matrix $\mathbf{H}\left(\omega \right)$ in Equation (14) corresponding to the gust input $w_g$ is proper and no split into proper and non-proper part is required. This is coherent with the fact that the Sears function for the gust input problem in the incompressible unsteady case tends to zero for $s\to \infty$ [48] (see Figure 5). Similarly, no residualization is required when considering the gust disturbance as unique system input [45]. In order to split the proper and non-proper parts for the rest of the columns of the aerodynamic transfer function matrix $\mathbf{H}\left(\omega \right)$ in a numerical stable way, a least-squares fit of a rational proper transfer function together with a polynomial part is done for each component of the transfer function matrix data (see Equation (16)):

(16)$\begin{array}{cc}\hfill H_{ij}\left(\omega \right)& =\sum _{k=1}^{N_{c,ij}}\frac{R_{ij}}{i\omega +p_{k,ij}}+P_{0,ij}+P_{1,ij}i\omega -P_{2,ij}{\omega }^2,\hfill \end{array}$

where $p_{k,ij},R_{ij}$ can either be real or come in conjugate pairs, $P_{0,ij},P_{1,ij},P_{2,ij}\in \mathbb{R}$ and the sum runs over the number $N_{c,ij}$ of finite poles $p_{k,ij}$ given by Equation (15) and which belong to the set $S_f$, $p_{k,ij}\in S_f$. The dependence of the finite poles with the indices $ij$ in $N_{c,ij}$ and $p_{k,ij}$ refers to the fact that each component of the transfer function matrix $\mathbf{H}\left(\omega \right)$ is dealt with separately. Note that, compared to the classical least-squares fit of the RFA technique which uses real coefficients [3,4,14], Equation (16) is able to represent peaks of the aerodynamic transfer functions in the frequency domain which appear in the transonic flow regime [49] by the presence of conjugate pair poles [6]. Additionally, the fit of Equation (16) does not suffer the ill-conditioning present when using an increasing number of real poles. In addition, the use of different poles for each element of the aerodynamic transfer function matrix $\mathbf{H}\left(\omega \right)$ as done here would lead to a prohibitive number of additional states if the classical approach of a least-squares fit were directly transformed into the time domain [50].

In the next step, the actual data after subtraction of the polynomial non-proper part instead of the strictly proper part of the fit given by the first term of Equation (16) is used. Thus, there is no error introduced in the proper part of the aerodynamic transfer function matrix and the representation of the non-proper part by the polynomial part can be seen as a residualization of the high frequency behaviour of the aerodynamic transfer function matrix $\mathbf{H}\left(\omega \right)$, as it is not available from the common aerodynamic solvers in the frequency domain (see Figure 4).

Next, the proper transfer function matrix ${\mathbf{H}}_a\left(\omega \right)$ is obtained, ${\mathbf{H}}_a\left(\omega \right)=\mathbf{H}\left(\omega \right)-{\mathbf{P}}_0-{\mathbf{P}}_{1}i\omega +{\mathbf{P}}_2{\omega }^2$. Similarly as done for the incompressible flow case, now the tangential interpolation method is applied to the proper transfer function matrix ${\mathbf{H}}_a\left(\omega \right)$ and the realization $\left({\mathbf{E}}_a=\mathbf{I},{\mathbf{A}}_a={\mathbf{E}}_r^{-1}{\mathbf{A}}_r,{\mathbf{B}}_a={\mathbf{E}}_r^{-1}{\mathbf{B}}_r,{\mathbf{C}}_a={\mathbf{C}}_r,{\mathbf{D}}_a={\mathbf{D}}_r\right)$ is obtained, where the SVD decomposition and the stabilization enforcement of Section 2.2 are applied afterwards. The matrices ${\mathbf{B}}_a$ and ${\mathbf{D}}_a$ together with the input vector have to be properly completed to accommodate the non-proper polynomial part (see Equation (17)). The matrices corresponding to the polynomial part have been split in column blocks corresponding to the translational motion (${\mathbf{P}}_{0t},{\mathbf{P}}_{1t},{\mathbf{P}}_{2t}$), rotational and flexible generalized coordinates (${\mathbf{P}}_{1rf},{\mathbf{P}}_{1rf},{\mathbf{P}}_{2rf}$) and control surface deflections (${\mathbf{P}}_{0c},{\mathbf{P}}_{1c},{\mathbf{P}}_{2c}$). The term corresponding to the third time derivative of the translational motion ${\stackrel{⃛}{\mathbf{u}}}_t$ is explicitly neglected by setting ${\mathbf{P}}_{2t}=0$. The submatrices in $\left[{\mathbf{B}}_{au1}{\mathbf{B}}_{ah2}{\mathbf{B}}_{ah1}{\mathbf{B}}_{au2}\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}\right]$ have been defined for convenience in a subsequent reordering of terms:

(17)$\begin{array}{c}\dot{{\mathbf{x}}_a}={\mathbf{A}}_a{\mathbf{x}}_a+\left[\begin{array}{cccccccccc}{\mathbf{B}}_{au1}& {\mathbf{B}}_{ah2}& {\mathbf{B}}_{ah1}& {\mathbf{B}}_{au2}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]{\mathbf{u}}_a\hfill \\ \left[\begin{array}{c}{\mathbf{c}}_l\\ {\mathbf{c}}_m\end{array}\right]={\mathbf{C}}_a{\mathbf{x}}_a+\left[\begin{array}{cccccccccc}\mathbf{0}& {\mathbf{P}}_{0t}& {\mathbf{P}}_{0rf}& {\mathbf{P}}_{0c}& {\mathbf{P}}_{1t}& {\mathbf{P}}_{1rf}& {\mathbf{P}}_{1c}& \mathbf{0}& -{\mathbf{P}}_{2rf}& -{\mathbf{P}}_{2c}\end{array}\right]{\mathbf{u}}_a,\hfill \end{array}$

where ${\mathbf{u}}_a={[w_g{\dot{\mathbf{u}}}_t^T{\mathbf{u}}_{rf}^T{\mathbf{u}}_c^T{\ddot{\mathbf{u}}}_t^T{\dot{\mathbf{u}}}_{rf}^T{\dot{\mathbf{u}}}_c^T{\stackrel{⃛}{\mathbf{u}}}_t^T{\ddot{\mathbf{u}}}_{rf}^T{\ddot{\mathbf{u}}}_c^T]}^T$. Dropping the terms corresponding to zero submatrices and reordering, the generalized state-space model for the unsteady aerodynamic loads is:

(18)$\begin{array}{c}{\dot{\mathbf{x}}}_a={\mathbf{A}}_a{\mathbf{x}}_a+\left[\begin{array}{cccc}{\mathbf{B}}_{au1}& {\mathbf{B}}_{au2}& \mathbf{0}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}{w}_g\\ {\mathbf{u}}_c\\ {\dot{\mathbf{u}}}_c\\ {\ddot{\mathbf{u}}}_c\end{array}\right]+\left[\begin{array}{cc}{\mathbf{B}}_{ah1}& {\mathbf{B}}_{ah2}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{rf}\\ \dot{{\mathbf{u}}_t}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{c}}_l\\ {\mathbf{c}}_m\end{array}\right]={\mathbf{C}}_a{\mathbf{x}}_a+\left[\begin{array}{cccc}\mathbf{0}& {\mathbf{P}}_{0c}& {\mathbf{P}}_{1c}& -{\mathbf{P}}_{2c}\end{array}\right]\left[\begin{array}{c}{w}_g\\ {\mathbf{u}}_c\\ {\dot{\mathbf{u}}}_c\\ {\ddot{\mathbf{u}}}_c\end{array}\right]+\left[\begin{array}{ccc}{\mathbf{P}}_{0rf}& {\mathbf{D}}_{ah2}& {\mathbf{D}}_{ah3}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{rf}\\ {\dot{\mathbf{u}}}_h\\ {\ddot{\mathbf{u}}}_h\end{array}\right],\hfill \end{array}$

where ${\mathbf{u}}_h={[{\mathbf{u}}_t^T{\mathbf{u}}_{rf}^T]}^T$ and:

$\begin{array}{c}{\mathbf{D}}_{ah2}=\left[\begin{array}{cc}{\mathbf{P}}_{0t}& {\mathbf{P}}_{1rf}\end{array}\right]\hfill \\ {\mathbf{D}}_{ah3}=\left[\begin{array}{cc}{\mathbf{P}}_{1t}& -{\mathbf{P}}_{2rf}\end{array}\right].\hfill \end{array}$

3.2. Aeroservoelastic System

In this section, the interaction between the aerodynamic and structural parts is considered in order to obtain a complete aeroservoelastic formulation. According to the virtual work principle, the forces transfer from the aerodynamic grid to the structural grid may be done with the corresponding transpose of the spline matrix which transfers the structural deformation to the aerodynamic grid [51]. However, and depending on the spline method for the transfer of forces from the aerodynamic grid to the structural grid, a set of nodal forces which is not smoothly distributed over the airframe but concentrated at particular positions can be obtained with such a transpose matrix. This is not a limitation when computing the generalized aerodynamic forces, but may lead to non realistic cut load distributions when applying the MDM method for the loads recovery. Assuming that the introduced error for not using the transpose of the spline matrix is negligible, the transfer of the aerodynamic forces to the structural grid is done by the summation matrix ${\mathbf{S}}_{gs}$, which transfers the aerodynamic coefficients acting over the lifting components to the nearest structural grid in the form of dimensional force and moment vectors. In order to obtain the total aerodynamic force acting over the aircraft structure, the resulting aerodynamic forces and moments are multiplied by the dynamic pressure q. The aeroservoelastic problem can be formulated in the time domain as:

$\begin{array}{cc}& {\mathbf{M}}_{hh}{\ddot{\mathbf{u}}}_h+{\mathbf{B}}_{hh}{\dot{\mathbf{u}}}_h+{\mathbf{K}}_{hh}{\mathbf{u}}_h\hfill \\ \hfill =& q{\mathit{\varphi }}_{gh}^T{\mathbf{S}}_{gs}\left({\mathbf{C}}_a{\mathbf{x}}_a+\left[\begin{array}{cccc}\mathbf{0}& {\mathbf{P}}_{0c}& {\mathbf{P}}_{1c}& -{\mathbf{P}}_{2c}\end{array}\right]\left[\begin{array}{c}{w}_g\\ {\mathbf{u}}_c\\ {\dot{\mathbf{u}}}_c\\ {\ddot{\mathbf{u}}}_c\end{array}\right]+\left[\begin{array}{ccc}{\mathbf{P}}_{0rf}& {\mathbf{D}}_{ah2}& {\mathbf{D}}_{ah3}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{rf}\\ {\dot{\mathbf{u}}}_h\\ {\ddot{\mathbf{u}}}_h\end{array}\right]\right),\hfill \end{array}$

where ${\mathbf{M}}_{hh}$, ${\mathbf{B}}_{hh}$ and ${\mathbf{K}}_{hh}$ represent the generalized mass, damping and stiffness matrices, respectively. The matrix ${\mathit{\varphi }}_{gh}$ contains a number $N_h$ (due to modal truncation) of the aircraft structural modes ${\mathbf{u}}_h$ (in vacuum) including the rigid-body and flexible modes. As in Section 3.1, the vector ${\mathbf{u}}_h={[{\mathbf{u}}_t^T{\mathbf{u}}_{rf}^T]}^T$ represents the total number of generalized coordinates (including rigid-body and structural flexible modes) and ${\mathbf{u}}_c$ the control surfaces. The subset ${\mathbf{u}}_t\subset {\mathbf{u}}_h$ represents the translational rigid-body degrees of freedom and the subset ${\mathbf{u}}_{rf}\subset {\mathbf{u}}_h$ the rotational rigid-body degrees of freedom together with the flexible modes. With no generalized stiffness corresponding to the translational degrees of freedom ${\mathbf{u}}_t,$ the relation ${\mathbf{K}}_{hh}{\mathbf{u}}_h={\mathbf{K}}_{hr}{\mathbf{u}}_{rf}={\mathbf{K}}_{hf}{\mathbf{T}}_f{\mathbf{u}}_{rf}$ can be used, where the subindex in the submatrix ${\mathbf{K}}_{hf}$ corresponds to the columns related to the subset of generalized coordinates representing the flexible modes ${\mathbf{u}}_f$, obtained with the matrix ${\mathbf{T}}_f$ from ${\mathbf{u}}_{rf}$.

In order to include the translational motion in the system output, the state vector $\mathbf{x}$ is augmented, $\widehat{\mathbf{x}}={\left[{\mathbf{u}}_t^T{\mathbf{u}}_{rf}^T{\dot{\mathbf{u}}}_h^T{\mathbf{x}}_a^T\right]}^T$, and the following equation of first order in the time derivatives is obtained:

(19)$\begin{array}{cc}\hfill \frac{d}{dt}\left(\left[\begin{array}{c}{\mathbf{u}}_t\\ {\mathbf{u}}_{rf}\\ {\dot{\mathbf{u}}}_h\\ {\mathbf{x}}_a\end{array}\right]\right)& =\left[\begin{array}{cccc}\mathbf{0}& \mathbf{0}& {\mathbf{T}}_t& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{T}}_{rf}& \mathbf{0}\\ \mathbf{0}& {\mathbf{A}}_{h1}& {\mathbf{A}}_{h2}& {\mathbf{A}}_{h3}\\ \mathbf{0}& {\mathbf{B}}_{ah1}& {\mathbf{B}}_{ah2}{\mathbf{T}}_t& {\mathbf{A}}_a\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_t\\ {\mathbf{u}}_{rf}\\ {\dot{\mathbf{u}}}_h\\ {\mathbf{x}}_a\end{array}\right]\hfill \\ & +\left[\begin{array}{cccc}\mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{0c}& {\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{1c}& -{\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{2c}\\ {\mathbf{B}}_{au1}& {\mathbf{B}}_{au2}& \mathbf{0}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}{w}_g\\ {\mathbf{u}}_c\\ {\dot{\mathbf{u}}}_c\\ {\ddot{\mathbf{u}}}_c\end{array}\right],\hfill \end{array}$

with:

$\begin{array}{c}{\mathbf{T}}_{hs}=q{\mathit{\varphi }}_{gh}^T{\mathbf{S}}_{gs},\hfill \\ \mathit{\alpha }={\mathbf{M}}_{hh}-{\mathbf{T}}_{hs}{\mathbf{D}}_{ah3},\hfill \\ {\mathbf{A}}_{h1}={\mathit{\alpha }}^{-1}\left(-{\mathbf{k}}_{hf}{\mathbf{T}}_f+{\mathbf{T}}_{hs}{\mathbf{P}}_{0rf}\right),\hfill \\ {\mathbf{A}}_{h2}={\mathit{\alpha }}^{-1}\left(-{\mathbf{B}}_{hh}+{\mathbf{T}}_{hs}{\mathbf{D}}_{ah2}\right),\hfill \\ {\mathbf{A}}_{h3}={\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{C}}_a,\hfill \end{array}$

where the matrix ${\mathbf{T}}_{rf}$ selects the non translational degrees of freedom ${\mathbf{u}}_{rf}$ from the generalized coordinates ${\mathbf{u}}_h$ (that is, excluding the subset ${\mathbf{u}}_t$ corresponding to the translational motion) and the matrix ${\mathbf{T}}_t$ selects the translational degrees of freedom from the vector of generalized coordinates ${\mathbf{u}}_h$.

Once the aeroservoelastic system has been written in a generalized state-space form of first order, the loads can be recovered by [52]:

• Mode displacement method (MDM), where the cut loads are given by:

  $\begin{array}{cc}\hfill {\mathbf{L}}_{MDM}& ={\mathbf{T}}_{cg}{\mathbf{K}}_{gg}{\mathit{\varphi }}_{gh}{\mathbf{u}}_h={\mathbf{T}}_{cg}{\mathbf{K}}_{gg}\left[\begin{array}{ccc}{\mathit{\varphi }}_{gt}& {\mathit{\varphi }}_{gr}& {\mathit{\varphi }}_{gf}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_t\\ {\mathbf{u}}_r\\ {\mathbf{u}}_f\end{array}\right]={\mathbf{T}}_{cg}{\mathbf{K}}_{gg}{\mathit{\varphi }}_{gf}{\mathbf{u}}_f,\hfill \end{array}$

  where the matrix ${\mathbf{K}}_{gg}$ corresponds to the physical stiffness matrix, ${\mathit{\varphi }}_{gf}$ includes the columns corresponding to the flexible modes and ${\mathbf{u}}_f$ their generalized coordinates. The matrix ${\mathbf{T}}_{cg}$ is a summation matrix which considers the nodal loads adding up to a particular location in order to obtain the resulting (incremental with respect to a steady-state reference) cut loads such as the shear forces together with the bending and torsional moments.

• Force summation method (FSM), where the cut loads are recovered by the equilibrium of forces:

  $\begin{array}{cc}\hfill {\mathbf{L}}_{FSM}& ={\mathbf{T}}_{cg}\left(-{\mathbf{M}}_{gg}{\mathit{\varphi }}_{gh}{\ddot{\mathbf{u}}}_h+q{\mathbf{S}}_{gs}{\mathbf{y}}_a\right),\hfill \end{array}$

  where the matrix ${\mathbf{M}}_{gg}$ represents the physical mass matrix. For the FSM method, not only the generalized aerodynamic forces but also the aerodynamic distribution ${\mathbf{y}}_a={\left[{\mathbf{c}}_l^T{\mathbf{c}}_m^T\right]}^T$ is required.

Both load recovery methods can be included in the output equation of the generalized state-space formulation of first order together with the aircraft response and the aerodynamic distribution:

(20)$\begin{array}{cc}\hfill \left[\begin{array}{c}{\mathbf{u}}_t\\ {\mathbf{u}}_{rf}\\ {\mathbf{y}}_a\\ {\mathbf{L}}_{MDM}\\ {\mathbf{L}}_{FSM}\end{array}\right]& =\left[\begin{array}{cccc}\mathbf{I}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{P}}_{0rf}+{\mathbf{D}}_{ah3}{\mathbf{A}}_{h1}& {\mathbf{D}}_{ah2}+{\mathbf{D}}_{ah3}{\mathbf{A}}_{h2}& {\mathbf{C}}_a+{\mathbf{D}}_{ah3}{\mathbf{A}}_{h3}\\ \mathbf{0}& {\mathbf{T}}_{cg}{\mathbf{K}}_{gg}{\mathit{\varphi }}_{gf}{\mathbf{T}}_f& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& q{\mathbf{T}}_{cg}{\mathbf{S}}_{gs}{\mathbf{P}}_{0rf}+\mathit{\beta }{\mathbf{A}}_{h1}& q{\mathbf{T}}_{cg}{\mathbf{S}}_{gs}{\mathbf{D}}_{ah2}+\mathit{\beta }{\mathbf{A}}_{h2}& q{\mathbf{T}}_{cg}{\mathbf{S}}_{gs}{\mathbf{C}}_a+\mathit{\beta }{\mathbf{A}}_{h3}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_t\\ {\mathbf{u}}_{rf}\\ {\dot{\mathbf{u}}}_h\\ {\mathbf{x}}_a\end{array}\right]\hfill \\ \hfill +& \left[\begin{array}{cccc}\mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{P}}_{0c}+{\mathbf{D}}_{ah3}{\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{0c}& {\mathbf{P}}_{1c}+{\mathbf{D}}_{ah3}{\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{1c}& -\left({\mathbf{P}}_{2c}+{\mathbf{D}}_{ah3}{\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}{\mathbf{P}}_{2c}\right)\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{\Gamma }{\mathbf{P}}_{0c}& \mathbf{\Gamma }{\mathbf{P}}_{1c}& -\mathbf{\Gamma }{\mathbf{P}}_{2c}\end{array}\right]\left[\begin{array}{c}{w}_g\\ {\mathbf{u}}_c\\ {\dot{\mathbf{u}}}_c\\ {\ddot{\mathbf{u}}}_c\end{array}\right],\hfill \end{array}$

where ${\mathbf{L}}_{MDM}$ and ${\mathbf{L}}_{FSM}$ represent the cut loads obtained by the MDM and FSM methods respectively and:

$\begin{array}{c}\mathit{\beta }=-{\mathbf{T}}_{cg}{\mathbf{M}}_{gg}{\mathit{\varphi }}_{gh}+q{\mathbf{T}}_{cg}{\mathbf{S}}_{gs}{\mathbf{D}}_{ah3},\hfill \\ \mathbf{\Gamma }=\mathit{\beta }{\mathit{\alpha }}^{-1}{\mathbf{T}}_{hs}+q{\mathbf{T}}_{cg}{\mathbf{S}}_{gs}.\hfill \end{array}$

Equations (19) and (20) represent the generalized state-space and output equations for the aeroservoelastic model description in the time domain.

In Figure 6 and Figure 7, the block diagrams corresponding to Equations (19) and (20) for both the MDM and FSM methods are shown respectively. In Figure 7, the matrix ${\mathbf{T}}_h$ selects the subset ${\mathbf{u}}_h$ from ${[{\mathbf{u}}_h^T{\mathbf{u}}_c^T]}^T$.

4. Application to the FERMAT Configuration

In this section, the FERMAT model created by Klimmek [53] has been used. The mass case C6 corresponding to a 100% of the fuel mass with a mass equal to the maximum takeoff weight $\mathrm{MTOW}$ = 260,000 (kg) and a center of gravity position in the x-direction (see Figure 8 for axes definition) with respect to the aircraft nose position $x_0$ of $x_{cg}-x_0=33.716$ (m) is considered. The spline from the structural model deformation to the aerodynamic grid is carried out by a thin-plate spline (TPS) method [54]. As described in Section 3.2, the TPS method generates a force transfer which may not smoothly distributed over the aircraft structure and thus the summation matrix ${\mathbf{S}}_{gs}$ from Equation (20) is used for the application of the FSM for the cut loads recovery. The cut loads presented in this section are incremental and do not include the loads acting at the 1g or level flight condition.

For the computation of the summation matrix ${\mathbf{S}}_{gs},$ the lifting surfaces are split up into strips, which are segments of constant width including all panels at a constant $\left(y,z\right)$ position, see Figure 8, where a typical strip aerodynamic element for the right wing component is shown and the panels conforming the control surface panels (ailerons) are highlighted. Neglecting the aerodynamic modeling of the fuselage component, there are a total of 1910 aerodynamic panels representing the complete configuration. The wing semi-span is 29.36 (m). A view of the right wing component is provided in Figure 9. Note that a strip can include panels from the wing component as well as from the aileron if they are topologically connected. The ailerons are taken into account in Section 4.2, where the effect of a gust loads alleviation (GLA) strategy is evaluated using the proposed aeroservoelastic model. In that case, the control surface modes corresponding to the ailerons are modeled with the transpiration approach as described in Section 3.1.

Regarding the dynamic structural model, two symmetric rigid-body modes corresponding to the heave and pitch (around the center of gravity) motion together with the 16 flexible modes ranging in frequencies from 2.03 (Hz) up to 7.30 (Hz) have been taken into account. The structural damping matrix has been set to zero, ${\mathbf{B}}_{hh}=\mathbf{0}$. The heave forms the vector ${\mathbf{u}}_t$ and the pitch together with the flexible modes form the vector ${\mathbf{u}}_{rf}$ in Equations (19) and (20). The aeroservoelastic system as described in Section 3.2 has been generated for the flightpoint corresponding to a Mach number of $M_{\infty }=0.85$ and an altitude $h=7640$ (m). The true airspeed is $U_{\infty }=263.147$ (m/s) and the reference length for the reduced frequency definition is $L_{ref}=7.005$ (m). The selected set of reduced frequencies for the Loewner realization has been chosen as $k=\left\{k_1\right\}\cup \left\{k_2\right\}\cup \left\{k_3\right\}$, where $k_1=\left\{0,\dots ,0.25\right\}$ with spacing $\mathrm{\Delta }{k}_1=0.005$, $k_2=\left\{0.29,0.335,0.575,0.415,0.46,0.5,0.585,0.665,0.75,0.8350,0.915,1\right\}$ and $k_3=\left\{1.165,1.335,1.5,1.665,1.835,2,2.165,2.335,2.5,2.665,2.835,3\right\}$. The interpolation of the precomputed AIC matrices to this set of reduced frequencies is done using a special linear interpolation technique [54]. The tolerance in Equation (15) has been set to $S_{\infty }=10\left(U_{\infty }/L_{ref}\right)$.

Figure 10 shows the element of the aerodynamic transfer function matrix corresponding to the (incremental) lift coefficient over the strip element at the wing station $y=13.48$ (m) shown in Figure 8 and Figure 9 in the complex plane due to a gust disturbance input for the different values of the reduced frequency k. In Figure 11, the pitch motion around the center of gravity is taken as input. For comparison purposes, the results obtained with six real poles applying the classical Roger’s approach (RFA) is also shown [4,5], where the poles $p_j$ have been computed as $p_j=-k_{max}/j$ for $j=1,\dots ,6$ with $k_{max}=3$. The data labeled as Loewner shows the present method and FD refers to the reference data in the frequency domain. It is clear from Figure 10 and Figure 11 that the present approach shows a better match compared to the Roger’s approach. Because of the finer resolution for the low reduced frequency range in the set of reduced frequencies selected, this better match is more evident in this region, which is of significant importance for aeroelastic applications.

4.1. Open Loop

In this section, the aeroservoelastic model is validated in open loop against the reference model available in the frequency domain. The model is excited by a 1-cosine gust $w_g\left(t\right)$ in (m/s) given by (at the aircraft nose) [1]:

$\begin{array}{c}{w}_g\left(t\right)=\frac{w_0}{2}\left(1-cos\left(\frac{\pi U_{\infty }}{H}t\right)\right),0\le t\le 2H/U_{\infty },\hfill \\ w_g\left(t\right)=0,t\geqq 2H/U_{\infty },\hfill \end{array}$

with $w_0$ the maximum gust amplitude (m/s), $U_{\infty }$ the translational speed (m/s) and H the gust gradient or half length (m). For this application, the gust amplitude has been set to ${\alpha }_{eq}=w_0/U_{\infty }$ equivalent to 1 (deg) and the gust gradient to $H=106.68$ (m) or $H=350$ (ft), the longest according to the regulations [1].

First, the rigid-body modes corresponding to the heave and pitch rotation about the center of gravity have been considered together with the first flexible mode of the structure (in vacuum) representing the symmetrical wing bending and with a natural frequency of $2.03$ (Hz), resulting in an aeroservoelastic model in the time domain with a total of 81 states. Even though the ailerons on the wing have been considered in the time-domain realization, they have been set to zero for the open loop validation. The corresponding incremental bending $MX$ and torsional $MY$ moments obtained by the FSM over the right wing component as predicted by both the aeroservoelastic model in the time domain (labeled as TD) and the reference model in the frequency domain (labeled as FD) are shown in Figure 12 and Figure 13 for the time instants at which the maximum and minimum vales are reached at the wing root location. It can be seen that the aeroservoelastic model with 81 states is able to reproduce the complete time history of the incremental cut loads distribution along the entire wing component.

Next, and in order to show the suitability of the present method for the cut loads recovery by means of the FSM, the effect of considering more flexible modes on the cut loads has been analyzed and a total number of 16 flexible modes ranging from frequencies $2.03$ up to $7.30$ (Hz) has been retained, obtaining a generalized state-space aeroservoelastic model with 231 states. Figure 14 shows the incremental bending moment $MX$ obtained in time domain (TD) with the resulting aeroservoelastic model with 231 states compared with the reference values computed in the frequency domain (FD) at the time instants where the maximum and minimum values are reached at the wing root. Again, both formulations are in agreement, which is due to the fact that the aerodynamic distribution can be represented in the time domain by the present method as detailed in Section 3.2, as shown in Figure 15 where the incremental lift coefficient distribution $c_l$ along the complete wing component is represented at the corresponding time instants.

4.2. Closed Loop

Once the aeroservoelastic model has been validated, the implementation of a simple GLA strategy is shown. An aeroservoelastic model including heave and pitch rigid-body modes, 16 flexible modes together with four ailerons (two inner and two outer) over the wing component has been chosen. As in Section 4.1, the same flightpoint has been selected and the total number of states used for the time-domain realization is 231. A simple proportional law between the gust speed at the aircraft nose and the outer ailerons, which are deployed symmetrically, has been implemented. An additional time delay between the gust speed detection and the aileron command of $t_d=0.06$ (s) has been taken into account. Thus, the vector ${\mathbf{u}}_c$ containing the control surface deflections is given in radian units by:

${\mathbf{u}}_c\left(t\right)={\left[\begin{array}{cccc}0& P{w}_g\left(t-t_d\right)& 0& P{w}_g\left(t-t_d\right)\end{array}\right]}^T,$

where the constant gain has been set to $P=-1/50$ (s/m) and the zero entries correspond to the inner ailerons. The maximum outer aileron deflection produced with this gain is 5.15 (deg). It is clear from Figure 16 that the maximum value of the incremental wing root bending moment $MX$ can be reduced in 19.01% at the expense of increasing the incremental torsional moment $MY$ in 14.01% (see Figure 17). Furthermore, the proposed aeroservoelastic model enables the evaluation of the effects of the applied GLA strategy on the cut loads over the complete aircraft structure.

5. Conclusions

In this work, an aeroservoelastic model for the prediction of the aeroelastic response together with the dynamic loads over the complete aircraft structure has been presented. It represents an interpolation instead of an approximation of the reference model in the frequency domain. The distributed aerodynamic forces created by the gust disturbance and induced by the aircraft motion can be obtained in the time domain and coupled with the structural model to obtain an aeroservoelastic state-space model of first order which includes the effect of control surfaces deflection. The force summation method can be used for the loads recovery, taking advantage of its superior convergence with the number of generalized coordinates in comparison with the mode displacement method.

The presented aeroservoelastic model can be used for further implementation of GLA strategies, including multivariable control design. As an application, a simple GLA logic with a proportional feedback closed loop architecture has been implemented to demonstrate the suitability of the present approach for aeroservoelastic design.

Related future work includes:

• Different GLA control law strategies. Within the presented method, the complete aerodynamic distribution together with the cut loads and combinations thereof can be chosen as objective functions.

• Consideration of different aerodynamic theories in the frequency domain which are appropriate for the transonic flow, such as the correction of the AIC matrices or linearized CFD solvers in the frequency domain. In particular, the piston theory as limit of the DLM method for high reduced frequency values may also be considered. In this case, the residualization could be eliminated by substracting the values predicted by the piston theory from the aerodynamic transfer function matrix.

• Inclusion of parametric generalized state-space aeroservoelastic models as an alternative to the classical gain scheduling approach.

• Extension to a nonlinear generalized state-space formulation for nonlinear aeroservoelastic systems. In this case, the Loewner framework in connection with a functional or Volterra series expansion theory can be followed.

Author Contributions

D.Q. developed and implemented the aeroservoelastic formulation and wrote the manuscript, P.V. and C.P.-V. implemented the Loewner framework algorithms and reviewed the manuscript.

Funding

This work has been funded within the frame of the Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology Demonstrator platform “AIRFRAME ITD” (contract N. CSJU-CS2-GAM-AIR-2015-15-01 Annex 1, Issue B04, 2 October 2015) being part of the Horizon 2020 research and Innovation framework programme of the European Commission.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Theodorsen Function Realization

Here, the matrices obtained after applying the Loewner realization given by Equation (5) to the Theodorsen function, followed by an SVD truncation retaining eight states, see Equation (7), and a projection into the $R{H}_{\infty }$ space in order to avoid the possible appearance of unstable poles, see Equation (8), are given:

$\begin{array}{cc}\hfill {\mathbf{A}}_c& =\left[\begin{array}{cccccccc}-70.2475& -15.2084& -4.3860& -1.5667& 0.6656& -0.2916& 0.1426& -0.0659\\ -127.6550& -60.8706& -23.0820& -8.5669& 3.5625& -1.5389& 0.7418& -0.3415\\ -168.3215& -105.6100& -59.3595& -26.9511& 11.5870& -4.9825& 2.3544& -1.0755\\ -225.3710& -146.1488& -100.3562& -60.3945& 30.6229& -13.8770& 6.5170& -2.9615\\ 315.4916& 200.8967& 143.6691& 101.8372& -64.7094& 34.1741& -16.5753& 7.5699\\ -451.2802& -282.2752& -200.3830& -148.7585& 109.8055& -72.2007& 40.6074& -19.7354\\ 681.3672& 420.8463& 294.1354& 217.1403& -167.1839& 127.0703& -87.6431& 50.0344\\ -967.5574& -594.0846& -411.6020& -301.1651& 232.6729& -187.3735& 150.9202& -108.2737\end{array}\right],\hfill \end{array}$

$\begin{array}{cc}\hfill {\mathbf{B}}_c& =\left[\begin{array}{c}-5.7098\\ -8.1427\\ -11.6790\\ -16.7748\\ 24.4703\\ -35.7673\\ 54.6494\\ -78.0051\end{array}\right],\hfill \end{array}$

$\begin{array}{cc}\hfill {\mathbf{C}}_c& =\left[\begin{array}{cccccccc}-5.7516& -0.9726& -0.3074& -0.1170& 0.0521& -0.0233& 0.0115& -0.0054\end{array}\right].\hfill \end{array}$

These matrices can then be used in Equation (13) to obtain a complete representation of the unsteady bidimensional incompressible flow over an airfoil with heave and pitch degrees of freedom. For an approximation with a number of states lower than eight, $n<8$, the first n rows and columns of the matrix ${\mathbf{A}}_c$ above should be taken. Similarly, the first n entries of the vectors ${\mathbf{B}}_c$ and ${\mathbf{C}}_c$ should be considered.

Figures

Figure 1. Heave uz and pitch uα pitch input.

Figure 2. Lift and pitch moment coefficients (cl,cm) as predicted by the system of ordinary differential equations (ODE) and the system of differential algebraic equations (DAE) of index 2. Reference results obtained in the frequency domain (FD) are also shown.

Figure 3. Theodorsen function in the complex plane and corresponding Jones approximation and Loewner realizations.

Figure 4. Typical element of the aerodynamic matrix transfer function (magnitude) against the reduced frequency.

Figure 5. Sears function Sk in the complex plane.

Figure 6. Mode displacement method (MDM).

Figure 7. Force summation method (FSM).

Figure 8. Typical strip element (blue) and control surfaces (red and orange) for the flutter reduced order assessment (FERMAT) configuration.

Figure 9. Typical strip element (blue) and control surfaces (red and orange) for the FERMAT configuration (view of the right wing component).

Figure 10. Element of the transfer function matrix from the gust input to the lift coefficient for the strip at the wing station y=13.48 (m).

Figure 11. Element of the transfer function matrix from the pitch input to the lift coefficient for the strip at the wing station y=13.48 (m).

Figure 12. Incremental bending moment MX over the right wing component at different time instants for M∞=0.85, h=7640 (m), H=106.68 (m) and αeq=1 (deg). Two rigid-body modes and one flexible mode considered. Results obtained by the FSM with the generalized state-space model in the time domain (TD) of 81 states compared to the reference model in the frequency domain (FD).

Figure 13. Incremental torsional moment MY over the right wing component at different time instants for M∞=0.85, h=7640 (m), H=106.68 (m) and αeq=1 (deg). Two rigid-body modes and one flexible mode considered. Results obtained by the FSM with the generalized state-space model in the time domain (TD) of 81 states compared to the reference model in the frequency domain (FD).

Figure 14. Incremental bending moment MX over the right wing component at different time instants for M∞=0.85, h=7640 (m), H=106.68 (m) and αeq=1 (deg). Two rigid-body modes and 16 flexible modes considered. Results obtained by the FSM with the generalized state-space model in the time domain (TD) of 231 states compared to the reference model in the frequency domain (FD).

Figure 15. Incremental lift coefficient cl over the complete wing component at different time instants for M∞=0.85, h=7640 (m), H=106.68 (m) and αeq=1 (deg). Two rigid-body modes and 16 flexible modes considered. Results obtained with the generalized state-space model in the time domain (TD) of 231 states compared to the reference model in the frequency domain (FD).

Figure 16. Time history of the incremental bending moment MX at the wing root with and without gust load alleviation (GLA) strategy.

Figure 17. Time history of the incremental torsional moment MY at the wing root with and without GLA strategy.