A frequency domain approach for reduced- order transonic aerodynamic modelling

Abstract

This paper describes a new efficient method for the construction of an approximately balanced aerodynamic Reduced Order Model (ROM) via the frequency domain using Computational Fluid Dynamics data. Time domain ROM construction requires CFD data, which is obtained from the DLR TAU RANS or Euler Linearised Frequency Domain (LFD) solver. The ROMs produced with this approach, using a small number of frequency simulations, are shown to exhibit a strong ability to reconstruct the system response for inviscid flow about the NLR7301 aerofoil and the FFAST wing; and viscous flow about the NASA Common Research Model. The latter demonstrates that the reduced order model approach can reconstruct the full order frequency response of a viscous aircraft configuration with excellent accuracy using a strip wise approach. The time domain models are built using the frequency domain, but also give promising results when applied to reconstruct non-periodic motions. Results are compared to time domain simulations, showing good agreement even with small ROM sizes, but with a substantial reduction in calculation time. The main advantage of the current model order reduction approach is that the method does not require the formation and storage of large matrices, such as in POD approaches.

Nomenclature

${\textbf{A}}$ , ${\textbf{B}},\;{\textbf{C}},{\textbf{D}}$

system matrices

$\widehat{\textbf{A}},\widehat{\textbf{B}},\widehat{\textbf{C}},\widehat{\textbf{D}}$

reduced system matrices

α

angle-of-attack

α ₀

amplitude of pitching motion

α _m

mean angle of pitching motion

c

aerofoil chord

C _p

pressure coefficient

C _L

lift coefficient

C _M

pitching moment coefficient

F _z

vertical force

$G$

transfer function

${\textbf{H}}$

Hankel matrix

$j$

$k $

reduced frequency

M _Y

pitching moment

$m$

number of inputs

$n$

system order

$p$

number of outputs

$s$

Laplace transform parameter

$t$

time

$T$

sample time interval

$\textbf{u}$

input vector

U_∞

freestream velocity

$\omega $

continuous domain frequency

$\widehat{\omega}$

discrete domain frequency

$\textbf{x}$

state vector

$\textbf{x}_{\textbf{d}}^{\ast}$

internal variable or state

$\textbf{y}$

output vector

$\textbf{y}^{\textbf{m}}$

modified output vector

$\widehat{y}$

discrete frequency output

$\tilde{y}$

continuous frequency output

$z$

Z transform parameter

Subscripts

d

discrete

1.0 Introduction

High-fidelity Computational Fluid Dynamics (CFD) simulations for unsteady aerodynamics remain too computationally expensive for routine use in aircraft design due to the high number of degrees of freedom and the thousands of parameter variations that require investigation. Therefore, historically simpler methods which are not able to predict all the features encountered in the flight regime of modern aircraft have been used e.g. doublet lattice method (DLM) [1, 2]. Thus, recently there has been a focus on Model Order Reduction (MOR) techniques to produce lower cost, but more accurate models of the unsteady CFD.

Physics-based MOR techniques have a long history for the modelling of physical systems, with early studies dating from the 1930s and 1940s [3–5]. In the field of CFD, construction of ROMs typically takes place in a discrete time or frequency domain since this is the form in which CFD codes are implemented. The efficiency of linear frequency domain solvers makes frequency domain ROMs creation particularly efficient [6]. Snapshot Proper Orthogonal Decomposition (POD) was first introduced by Sirovich in 1987 [7] and has become widely used in fluid dynamics [8–10]. The method projects large CFD solutions onto a smaller subspace such that the projection error is minimised and efficiently finds the projection kernel. Implementation of the approach requires a suitable selection of frequency solutions or time domain simulations [11] to create the snapshots. If a balanced truncation is used, then the method is called BPOD [8] and it can be shown that this approach has links to another popular approach based on the Eigensystem Realisation Algorithm (ERA) [12–15]. ERA builds on the work of Kung et al. [16] and has been applied to a range of fluid problems including flow control [12], gust modelling [17] and flutter analysis [18]. ERA uses forward simulation in the time domain, and BPOD requires forward and backward (adjoint) simulation in the time domain. This means that BPOD produces the normal and adjoint modes, whereas ERA only produces the normal modes.

Non-intrusive artificial intelligence or machine learning methods are also candidates to produce reduced order models of unsteady CFD, but there are many issues that remain unresolved with purely data driven models [19], including the amount of training data needed for accuracy and the associated high computational costs. Instead, physics informed methods can make use of physical laws to constrain the admissible solution space so that for example solutions that do not conserve mass can be discarded [20]. Hybrid MOR methods that combine physics-based ROMs with machine learning have also been proposed and this is an emerging area of study [21, 22], which will incorporate developments in both fields.

The work described in this paper is a physics-based MOR method that modifies the system identification method of McKelvey et al. [23], which can be adapted for MOR. One algorithm based on equi-spaced frequencies and a further algorithm with variable frequency spacing are presented [23]. Vakilzadeh et al. [24] note that the latter method of system identification can be numerically sensitive to the frequencies used and discuss how the excitation frequencies should be selected to obtain a well-conditioned estimation problem. The main advantage of the current model order reduction approach is that the method is balanced and does not require the formation and storage of large matrices, such as in POD approaches. However, unlike balanced POD [8] the method does not allow the reconstruction of adjoint modes.

Here the system matrix ${\textbf{D}}$ is preserved on MOR as it has a physical meaning for the unsteady fluid flow, which would be lost in the basic algorithm of McKelvey et al. [23]. The method can be used to calculate both periodic and non-periodic flows and results are presented for inviscid flow about the NLR7301 aerofoil and the FFAST wing, and viscous flow about the NASA Common Research Model aircraft. The necessary CFD data is obtained from the DLR TAU code [25–27].

2.0 Systems

For many unsteady motions in aerospace engineering, the dynamic system behaviour is nearly linear about a non-linear mean flow. Therefore, the RANS or Euler equations in spatially discrete form in a CFD code can be approximated by a MIMO linear time-continuous time invariant state space system of $n$ -th order, with $m$ inputs and $p$ outputs, where the future behaviour of the system depends on its past, its present condition and its inputs:

(1) \begin{align}\dot{\textbf{x}}(t) &= {\textbf{A}}\textbf{x}(t) + {\textbf{B}}\textbf{u}(t)\nonumber \\[5pt] \textbf{y}(t) & = {\textbf{C}}\textbf{x}(t) + {\textbf{D}}\textbf{u}(t)\end{align}

where $\textbf{y}(t)\in \mathbb{R}^{p}$ is the vector of outputs of interest from the system, $\textbf{u}(t)\in \mathbb{R}^{m}$ is the vector of inputs and $\textbf{x}(t)\in \mathbb{R}^{n}$ is the state matrix. The time invariant system matrices have the following sizes: $\textbf{A}\in \mathbb{R}^{n\times n}$ , $\textbf{B}\in \mathbb{R}^{n\times m}$ , $\textbf{C}\in \mathbb{R}^{p\times n}$ and $\textbf{D}\in \mathbb{R}^{p\times m}$ . It should be noted that a reduced order model reduces the large dimension $n$ , but $p$ and $m$ remain unchanged. Thus for the current CFD application matrices ${\textbf{A}}$ , ${\textbf{B}}$ and C will have a smaller dimension in the ROM, but the matrix D will stay the same small size $p \times m$ and the reduced order model will be forced to keep this matrix unchanged as it has a physical meaning, see Section 3.

2.1. Continuous frequency domain solution

Then using Laplace transforms it can be shown that the transfer function of Equation (1) is

(2) \begin{align} G(s) = {\textbf{C}}{({s{\textbf{I}} - {\textbf{A}}})^{-1}}{\textbf{B}} + {\textbf{D}} \end{align}

If a frequency domain solution of (1) is sought with $\textbf{u}(t) = {e^{j\omega t}}$ then the frequency domain output solution $\tilde{y}$ where $\textbf{y}(t)=\tilde{y}e^{j\omega t} $ is given by

(3) \begin{equation} \tilde{y}=\textbf{C}(j\omega\textbf{I}-\textbf{A})^{-1}\,\textbf{B}+\textbf{D}\end{equation}

and will equal the transfer function provided that $s = j\omega $ .

2.2. Discrete time system

One approach to obtain a discrete time system approximation is to integrate (1) from $t = kT$ to $ = \!\left( {k + 1} \right)T$ :

(4) \begin{align}\textbf{x}({k + 1}) = \textbf{x}(k) + \mathop \int \nolimits_{kT}^{({k + 1})T} ({{\textbf{A}}\textbf{x}(t) + {\textbf{B}}\textbf{u}(t)})dt \end{align}

where $\textbf{x}(k)$ indicates the exact continuous-time solution at time level $kT$ . Then approximating the integral using the trapezium rule and introducing the variable $\textbf{x}_{\textbf{d}}(k)$ which approximates the sampled state $\textbf{x}(k)$ yields the following discrete system

(5) \begin{align}{\textbf{x}_{\textbf{d}}}\!\left( {k + 1} \right) &= {\textbf{x}_{\textbf{d}}}(k) + \frac{T}{2}\!\left( {{\textbf{A}}{\textbf{x}_{\textbf{d}}}\!\left( {k + 1} \right) + {\textbf{B}}\textbf{u}\!\left( {k + 1} \right)} \right) + \frac{T}{2}\!\left( {{\textbf{A}}{\textbf{x}_{\textbf{d}}}(k) + {\textbf{B}}\textbf{u}(k)} \right)\nonumber\\ {\textbf{y}_\textbf{d}}(k) &= {\textbf{C}}{\textbf{x}_{\textbf{d}}}(k) + {\textbf{D}}\textbf{u}(k)\end{align}

Note here that even if $\textbf{u}(t)$ is piecewise linear so that the integration of $\textbf{u}$ is exact, $\textbf{x}_{\textbf{d}}(k) \ne \textbf{x}(k)$ because the integration of $\textbf{x}(t)\;$ in Equation (5) is also achieved using the trapezium rule. Following Ref. (28) a midpoint series is first introduced and then a new internal variable

(6) \begin{equation}\textbf{x}_{\textbf{d}}^*(k) = \frac{1}{{\sqrt T }}\!\left( {{\textbf{I}} - \frac{T}{2}{\textbf{A}}} \right){\textbf{x}_{\textbf{d}}}(k) - \frac{{\sqrt T }}{2}\;{\textbf{B}}\;\textbf{u}(k)\end{equation}

Then, simplifying leads to the following system of equations

(7) \begin{align} \textbf{x}_{\textbf{d}}^{\ast}\!\left( {k + 1} \right) &= {{\textbf{A}}_d}\textbf{x}_{\textbf{d}}^{\ast}(k) + {{\textbf{B}}_d}\textbf{u}(k)\nonumber \\{\textbf{y}_{\textbf{d}}}(k) & = {{\textbf{C}}_d}\textbf{x}_d^{\ast}(k) + {{\textbf{D}}_d}\textbf{u}(k)\end{align}

where $\textbf{y}_{\textbf{d}}(k)$ is a discrete approximation to $\textbf{y}(t)$ at $t = kT$ and the discrete system matrices are

(8) \begin{align} {{{\textbf{A}}_d}} & = {\!\left( {\frac{2}{T}{\textbf{I}} + {\textbf{A}}} \right){\rm{\;}}{{\!\left( {\frac{2}{T}{\textbf{I}} - {\textbf{A}}} \right)}^{-1}}}\nonumber \\[5pt]{{{\textbf{B}}_d}} & = {\frac{2}{{\sqrt T }}{{\left( {\frac{2}{T}{\textbf{I}} - {\textbf{A}}} \right)}^{-1}}{\textbf{B}}}\nonumber \\[5pt]{{{\textbf{C}}_d}} &= {\frac{2}{{\sqrt T }}{\textbf{C}}{{\left( {\frac{2}{T}{\textbf{I}} - {\textbf{A}}} \right)}^{-1}}}\nonumber \\[5pt]{{{\textbf{D}}_d}} &= {{\textbf{C}}{{\left( {\frac{2}{T}{\textbf{I}} - {\textbf{A}}} \right)}^{ - 1}}{\textbf{B}} + {\textbf{D}}}\end{align}

This form for the discrete system matrices is also given by Al-Saggaf and Franklin [29], who show that the corresponding continuous-time system is controllable and observable if and only if the discrete-time system is also controllable and observable. Further a balanced truncation in one domain will map to a balanced truncation in the other domain [23]. This transformation is called the bilinear, Cayley or Tustin transformation and the inverse is given by:

(9) \begin{align}{{\textbf{A}}} &= {\frac{2}{T}\!\left( {{\textbf{I}} + {{\textbf{A}}_d}} \right){^{-1}}\!\left( {{{\textbf{A}}_d} - {\textbf{I}}} \right)}\nonumber \\[5pt]{{\textbf{B}}} & = {\frac{2}{{\sqrt T }}\!\left( {{\textbf{I}} + {{\textbf{A}}_d}} \right){^{-1}}{{\textbf{B}}_d}}\nonumber \\[5pt]{{\textbf{C}}} &= {\frac{2}{{\sqrt T }}{{\textbf{C}}_d}\!\left( {{\textbf{I}} + {{\textbf{A}}_d}} \right){^{-1}}}\nonumber\\[5pt]{\textbf{D}} & = {{{\textbf{D}}_d} - {{\textbf{C}}_d}\!\left( {{\textbf{I}} + {{\textbf{A}}_d}} \right){^{ - 1}}{{\textbf{B}}_d}}\end{align}

2.3. Discrete frequency domain solution

For the discrete system (7) the transfer function can be found directly or from Equation (5) since the input/output relationship is unchanged by the introduction of a new internal variable. Hence the transfer function of the system (7) is given by

(10) \begin{equation}{G_d}\!\left( z \right) = {{\textbf{C}}_d}{\left( {z{\textbf{I}} - {{\textbf{A}}_d}} \right)^{ - 1}}{{\textbf{B}}_d} + {{\textbf{D}}_d}\end{equation}

or alternatively

(11) \begin{equation}{G_d}\!\left( z \right) = {\textbf{C}}{\left[ {\frac{2}{T}\frac{{z - 1}}{{z + 1}}{\textbf{I}} - {\textbf{A}}} \right]^{ - 1}}{\textbf{B}} + {\textbf{D}}\end{equation}

Comparing Equations (11) and (2) it can be observed that the discrete and continuous transfer functions are identical if

(12) \begin{align}s = \frac{2}{T}\!\left( {\frac{{z - 1}}{{z + 1}}} \right){\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}or{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}z = \frac{{1 + sT/2}}{{1 - sT/2}}\end{align}

If it is assumed that the discrete input and outputs have the form

(13) \begin{align}\textbf{x}_{\textbf{d}}^*(k)&=\hat{x}e^{j\hat{\omega}k}\nonumber \\\textbf{u}(k)&=e^{j\hat{\omega}k}\\\textbf{y}_{\textbf{d}}(k)&=\hat{y}e^{j\hat{\omega}k} \nonumber\end{align}

where $\widehat{\omega}$ is the discrete frequency, then the discrete system (7) becomes

(14) \begin{align}\hat{x}e^{j\hat{\omega}(k+1)}&=\textbf{A}_d \hat{x}e^{j\hat{\omega}k}+\textbf{B}_d e^{j\hat{\omega}k}\\[5pt]\hat{y}e^{j\hat{\omega}k}&=\textbf{C}_d \hat{x}e^{j\hat{\omega}k}+\textbf{D}_d e^{j\hat{\omega}k} \end{align}

and it can be shown that

(15) \begin{equation}\hat{y}=[\textbf{C}_d(e^{j\hat{\omega}}\textbf{I}-\textbf{A}_d)^{-1}\textbf{B}_d+\textbf{D}_d]\end{equation}

Comparing to Equation (10), if $z=e^{j\hat{\omega}}$ then the discrete frequency domain output function $\widehat{y}$ is equal to the transfer function ${{\mathrm{G}}_d}(z).$ Consequently, it is possible to link the discrete and continuous frequency solutions using Equation (12) so $\hat{y}=\tilde{y}$ if

(16) \begin{equation}\omega=\frac{2}{T}tan\left({\frac{\hat{\omega}}{2}}\right)\end{equation}

The reduced order modelling approach described here starts from a set of discrete frequencies and Equation (16) allows the required transfer functions to be calculated from Linear Frequency Domain (LFD) solver within DLR code TAU, which works in the continuous frequency domain. The aim of ROMs is to build them from existing, widely used codes with minimal changes and this link between solutions in various domains and transfer functions enables that approach.

3.0 Reduced order models

A discrete time-domain ROM is first constructed using frequency domain CFD data. After the discrete-time ROM has been created, a corresponding continuous-time ROM is obtained using Equation (9), allowing the simulation of any motion in the time domain, including non-periodic motions. The reduced order modelling technique developed in this work is based on a method developed for system identification by McKelvey et al. [23]. It should be noted that a key difference between the application here, and that targeted in Ref. (23) is that the system here is fully known and hence the application is one of Model Order Reduction (MOR) rather than identification. This means that for the current application the system matrix ${\textbf{D}}$ is a small known matrix and has a physical meaning in that it is a linear term resulting from the integration of the mean pressure over the perturbed mesh due to the aerofoil motion. In the methods described in Ref. (23) the matrix ${\textbf{D}}$ would be identified and the link to the physical meaning of this term in the final reduced order model would be lost. There is no requirement for knowledge of the high order system matrices, which are often not explicitly available in CFD codes or when using this type of method with experimental data.

Here a method based on equi-spaced frequencies is used as it is shown in Ref. (23) to be strongly consistent if the noise has zero mean and a covariance function which is uniformly bounded. In the present study, noise is not an issue as the full order system is known. The long-term aim of the work is an automated method that does not require an expert user. The idea is thus to start from an equi-spaced ROM, where the user has only to decide on a frequency range and number of frequencies to use in ROM creation. This approach is used since preliminary studies with a variable spacing approach has shown that for the current target application the quality of the ROM produced is highly dependent on the selected frequencies and hence on expert knowledge. This issue of sensitivity to the frequencies is also reported by Vakilzadeh et al. [24]. In future, a monitoring parameter will be identified such that additional frequencies can be added in where the frequency domain solution is changing rapidly in an automated manner.

In what follows a uniformly spaced set of discrete frequencies between 0 and $\pi$ is specified and the method is identical to that of McKelvey et al. [23], until the method deviates to define a new approach to defining $\widehat{\textbf{B}}_d$ and $\widehat{\textbf{D}}_d$ , the reduced order system matrices. The frequencies are:

(17) \begin{equation}\hat{\omega}_K=\frac{K\pi}{N},\qquad K\in[0,N]\end{equation}

If the corresponding frequency domain solution $\hat{y}_K$ could be found, then the discrete transfer function $\hat{y}_K=G_d(e^{j\omega K})=G_d(K)$ would be known. The transfer functions obtained for $\widehat{\omega}_K \in [0,\pi]$ can be used to extend the model input data to the interval $\widehat{\omega}_K \in [\pi,2\pi]$ since

(18) \begin{equation}{G_d}\!\left( {K + N} \right) = {\rm{\;}}G_d^{{\rm{*\;}}}\!\left( {N - K} \right),\qquad K{\rm{\;}} \in {\rm{\;}}\left[ {1{\rm{\;}},{\rm{\;}}N - 1} \right]\end{equation}

where $\widehat{\omega}_{K+N}=\frac{(K+N)\pi}{N},K \in [1, N-1]$ . Thus, the whole unit circle is utilised within the algorithm [23] and the frequency response is sampled at $2N$ points. The solutions $\hat{y}_K$ are not calculated directly, instead these are found from continuous frequency domain solutions $\tilde{y}_K$ calculated using the DLR-TAU code using the correct corresponding continuous frequency specified by Equation (16). The sampled frequency response data can be used with a 2N-point Inverse Discrete Fourier Transform (IDFT) to obtain entries for a Hankel matrix. The coefficients of the IDFT are defined by

(19) \begin{align}h_i & = \frac{1}{2N}\ \sum^{2N-1}_{K=0}{{{\rm G}}_d(K) e^{j2\pi iK{\rm /2}N}}=\frac{1}{2N} \sum^{2N-1}_{K=0}{{{\rm G}}_d(K) e^{j{\widehat{\omega }}_K{i/2}N}} \nonumber \\[5pt] & \qquad\qquad\qquad\quad i\in {[0\ ,\ 2N -1]} \end{align}

It should be noted that as $N \to \infty $

(20) \begin{align}\begin{array}{*{20}{c}}{\mathop {\lim }\limits_{N \to \infty } {h_i} = {\rm{\;}}g\!\left( i \right){\rm{\;\;}}}\\{{\rm{\;\;}}}\end{array}\end{align}

where

(21) \begin{align} g(k) = \left\{ \begin{array}{c@{\quad}c} \textbf{D}_{d} & k = 0 \\[5pt] \textbf{C}_{d}\textbf{A}_{d}^{k-1}\textbf{B}_{d} & k > 0 \end{array} \right. \end{align}

i.e. the terms tend to the time-domain unit impulse responses of the discrete linear system. The ${h_{i\;}}$ terms are used to construct a block, real Hankel matrix, which has $l$ block rows and $\;s$ block columns

(22) \begin{align}\textbf{H}&=\!\left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}h_1 & h_2 & \cdots & h_s \\h_2 & h_3 & \cdots & h_{s+1} \\\vdots & \vdots & \ddots & \vdots \\h_l & h_{l+1} & \cdots & h_{l+s-1} \end{array}\right) {\in {\mathbb R}}^{lp{\rm \times }sm}\\ & {\rm{\;with\;}}l \gt r,{\rm{\;}}s \ge r{\rm{\;and\;}}l + s \le 2N\nonumber\end{align}

where $r$ is the rank of the discrete-time ROM to be derived. Equation (20) means that for $N \to \infty $ this tends to the Hankel matrix used in other reduced order model methods [12, 13, 17]. Following a Singular Value Decomposition (SVD) of ${\textbf{H}}$ , model reduction is achieved by keeping the $r$ largest singular values. This is done by partitioning ${\textbf{U}}$ , ${\textbf{V}}$ and ${\boldsymbol\Sigma}$ so that

(23) \begin{align}\textbf{H} = \left[ \textbf{U}_{r}{\textbf{U}}_{o} \right] \left[ \begin{array}{c@{\quad}c}{\boldsymbol\Sigma }_{r} & 0\\0 & {\boldsymbol\Sigma}_{o}\end{array} \right] \left[ \begin{array}{c}\textbf{V}_{r}^{T}\\[5pt]\textbf{V}_{o}^{T}\end{array} \right]\end{align}

where the diagonal block matrix ${\boldsymbol\Sigma}_{r}\;$ contains the $r$ largest singular values, $r$ being the chosen size of the reduced order model. A reduced order approximation to the discrete-time linear state-space system then has the form

(24) \begin{align}{\widehat{{\textbf{x}}}}^{{\,\boldsymbol\ast}}_{{\textbf{d}}}\!\left(k+1\right) & ={\widehat{{\textbf A}}}_d{\widehat{{\textbf{x}}}}^{{\,\boldsymbol\ast}}_{{\textbf{d}}}\!\left(k\right)+{\widehat{{\textbf B}}}_d{\textbf{u}}\!\left(k\right) \nonumber \\[5pt]{\widehat{{\textbf{y}}}}_{{\textbf{d}}}\!\left(k\right) & ={\widehat{{\textbf{C}}}}_d{\widehat{{\textbf{x}}}}^{{\,\boldsymbol\ast}}_{{\textbf{d}}}\!\left(k\right)+{\widehat{{\textbf D}}}_d{\textbf{u}}(k) \end{align}

where ${{\widehat{{\textbf{x}}}}^{{\textbf *}}_{{\textbf{d}}}\in {\mathbb R}}^r$ , ${{\widehat{{\textbf A}}}_d\in {\mathbb R}}^{r{\rm \times }r}$ , ${\widehat{{\textbf B}}}_d{{\rm \ }\in {\mathbb R}}^{r{\rm \times }m}, {\widehat{{\textbf C}}}_d{{\rm \ }\in {\mathbb R}}^{p{\rm \times r}}, {\widehat{{\textbf D}}}_d{{\rm \ }\in {\mathbb R}}^{m{\rm \times }p}$ , and ${\widehat{{\textbf{y}}}}_{{\textbf{d}}}{\in {\mathbb R}}^p$ is the model approximation of ${\textbf{y}_{\textbf{d}}}\;.$

The new system identification and reduction scheme for $\widehat{\textbf{A}}_d, \widehat{\textbf{B}}_d, \widehat{\textbf{C}}_d$ and $\widehat{\textbf{D}}_d$ is modified from Ref. (23) to ensure that $\widehat{\textbf{D}}_d$ the system matrix maintains the physical meaning it has in the original unreduced system. The method is applied here to a known CFD system, where the matrices ${\textbf{A}}$ , ${\textbf{B}}$ and ${\textbf{C}}$ are large and not explicitly constructed within the code. The matrix ${\textbf{D}}$ is however a small known matrix of size $p \times m$ and the algorithm is modified so that the continuous time reduced order model matrix $\widehat{\textbf{D}}=\textbf{D}$ First define a modified output variable

(25) \begin{align}\begin{array}{*{20}{c}}{{\textbf{y}^m}(t) = \textbf{y}(t) - {\textbf{D}}\textbf{u}(t){\rm{\;\;}}}\\{{\rm{\;\;}}}\end{array}{\rm{\;}}\end{align}

So that effectively equation (1) becomes

(26) \begin{align}\begin{array}{c}\dot{\textbf{x}}(t) = {\textbf{A}}\textbf{x}(t) + {\textbf{B}}\textbf{u}(t)\\[5pt]{\textbf{y}^m}(t) = {\textbf{C}}\textbf{x}(t)\end{array}\end{align}

i.e. for this modified system ${\textbf{D}} = {\textbf{0}_{p \times m}}$ . The discretisation process is the same as described above, but ${\textbf{D}} = {\textbf{0}_{p \times m}}$ in the equations. The discrete and reduced matrices $\widehat{\textbf{A}}_d$ and $\widehat{\textbf{C}}_d$ are then as given by (23):

(27) \begin{align}\begin{array}{c}{\widehat{{\textbf A}}}_d=\dfrac{{\left({{\textbf J}}_1\ {{\textbf U}}_r\right)}^{\dagger }\left({{\textbf J}}_1\ {{\textbf U}}_r\right)}{{\left({{\textbf J}}_1\ {{\textbf U}}_r\right)}^{\dagger }\left({{\textbf J}}_2\ {{\textbf U}}_r\right)}\ \ and\ {\widehat{{\textbf C}}}_d=\left({{\textbf J}}_3\ {{\textbf U}}_r\right) \\\ \ \end{array}\end{align}

where

(28) \begin{align}{{\textbf{J}}_1} &= \left[ {{{\textbf{I}}_{\!\left( {l - 1} \right)p}}{\rm{\;}}{\textbf{0}_{\!\left( {l - 1} \right)p \times p}}} \right]\nonumber \\[4pt] {{\textbf{J}}_2} & = \left[ {{\textbf{0}_{\!\left( {l - 1} \right)p \times p}}{\rm{\;}}{{\textbf{I}}_{\!\left( {l - 1} \right)p}}} \right]\nonumber \\[4pt] {{\textbf{J}}_3} &= \left[ {{{\textbf{I}}_p}{\rm{\;}}{\textbf{0}_{p \times \!\left( {l - 1} \right)p}}} \right]\end{align}

In these equations ${{\textbf{I}}_i}\;$ is a $\!\left( {i \times i} \right)$ identity matrix, ${\textbf{0}_{i \times j}}\;$ is a $\!\left( {i \times j} \right)$ matrix with zero entries and ${{\textbf{X}}^\dagger } = {\!\left( {{{\textbf{X}}^T}{\textbf{X}}} \right)^{ - 1}}{{\textbf{X}}^T}$ is the Moore-Penrose pseudoinverse of matrix ${\textbf{X}}$ . The process for defining $\widehat{\textbf{B}}_d$ and $\widehat{\textbf{D}}_d$ needs to ensure that when using the inverse transformation (9) with the reduced order matrices, the continuous time reduced order modified system matrix $\widehat{\textbf{D}}$ also has $\widehat{\textbf{D}}=\textbf{0}_{p\times m}$ so that

(29) \begin{align}\widehat{{\textbf{y}}}\left(t\right)={\widehat{{\textbf{y}}}}^m\left(t\right)+{\textbf D}{\textbf{u}}(t) \end{align}

Then from Equation (9) require that

(30) \begin{align}\begin{array}{c}{\widehat{{\textbf D}}=\widehat{{\textbf D}}}_d-{\widehat{{\textbf C}}}_d{\left({\textbf I}+{\widehat{{\textbf A}}}_{{\textbf{d}}}\right)\ }^{-1}{\widehat{{\textbf B}}}_d={{\textbf 0}}_{p\times m}\ \\\ \ \end{array}\end{align}

which means that the following relationship must hold

(31) \begin{align}\begin{array}{c}{\widehat{{\textbf D}}}_d={\widehat{{\textbf C}}}_d{\left({\textbf I}+{\widehat{{\textbf A}}}_{{\textbf{d}}}\right)\ }^{-1}{\widehat{{\textbf B}}}_d\ \ \\\ \ \end{array}\end{align}

So the process for calculating $\widehat{\textbf{B}}_d$ is modified from that described in Ref. (23) as follows. For $\;N\; \ge \;r$ , then the matrix ${\boldsymbol{\Omega }}$ is here defined as

(32) \begin{align}\begin{array}{c}{\boldsymbol\Omega }=\left[ \begin{array}{c}{\widehat{{\textbf C}}}_d\left[{\left(z_0{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}+{\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}\right] \\[4pt]{\widehat{{\textbf C}}}_d\left[{\left(z_1{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}+{\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}\right] \\[4pt]\vdots \\{\widehat{{\textbf C}}}_d\left[{\left(z_N{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}+{\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}\right] \end{array}\right]\in \ {{\mathbb R}}^{\left(N+1\right)p\times r},{\rm \ }\ \\ \\{\ z}_K=e^{j{\omega }_Kt}{\rm \ }\forall {\rm \ K\ }\in {\rm [0\ ,\ N]} \end{array}\end{align}

Then it can be observed that

(33) \begin{align} \begin{array}{c}{\boldsymbol \Omega }{\widehat{{\textbf B}}}_d=\left[ \begin{array}{c}{\widehat{{\textbf C}}}_d{\left(z_0{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{{\widehat{{\textbf C}}}_d\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d \\[4pt] {\widehat{{\textbf C}}}_d{\left(z_1{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{{\widehat{{\textbf C}}}_d\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d \\[4pt] \vdots \\[4pt] {\widehat{{\textbf C}}}_d{\left(z_N{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{\widehat{{\textbf C}}}_d{\left({\textbf I}{\textbf +}{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d \end{array}\right]\ \ \\\ \ \end{array}\end{align}

and using Equation (30)

(34) \begin{align} \begin{array}{c}{\boldsymbol \Omega }{\widehat{{\textbf B}}}_d=\left[ \begin{array}{c}{\widehat{{\textbf C}}}_d{\left(z_0{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{\widehat{{\textbf D}}}_d \\[4pt] {\widehat{{\textbf C}}}_d{\left(z_1{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{\widehat{{\textbf D}}}_d \\[4pt] \vdots \\[4pt] {\widehat{{\textbf C}}}_d{\left(z_N{\textbf I}-{\widehat{{\textbf A}}}_d\right)}^{-1}{\widehat{{\textbf B}}}_d+{\widehat{{\textbf D}}}_d \end{array}\right]=\widehat{{\mathcal G}}\ \ \\\ \ \end{array}\end{align}

where

(35) \begin{align} \begin{array}{c}\widehat{{\mathcal G}}={\left[G_d(0),\ G_d(1),\ G_d(2),\cdots G_d(N) \right]}^T\ \\\ \ \end{array}\end{align}

Now introducing $\breve{\boldsymbol{\Omega }}$ , which consists of the real and imaginary parts of ${\boldsymbol{\Omega }}$ .

(36) \begin{align} \begin{array}{*{20}{c}}{\mathop {\breve{\boldsymbol{\Omega }}}\limits = {\rm{\;}}\left[ {\begin{array}{*{20}{c}}{Re\!\left( {\boldsymbol{\Omega }} \right){\rm{\;}}}\\{Im\!\left( {\boldsymbol{\Omega }} \right)}\end{array}} \right]{\rm{\;\;}}}\\{{\rm{\;\;}}}\end{array}\end{align}

means $\widehat{\textbf{B}}_d \in \mathbb{R}^{r\times m}$ can are calculated by

(37) \begin{align} \begin{array}{c}{\ \widehat{{\textbf B}}}_d=({\check{{\boldsymbol \Omega^T }}}\check{{\boldsymbol \Omega }}\ )^{-1}{\check{{\boldsymbol \Omega^T }}}\ \left[ \begin{array}{c}Re(\widehat{{\mathcal G}})\ \\Im(\widehat{{\mathcal G}}) \end{array}\right]\ \\\ \ \end{array}\end{align}

and $\widehat{\textbf{D}}_d \in \mathbb{R}^{p\times m}$ can be found from Equation (30) if the discrete scheme matrices are required. From these it is possible to find continuous reduced order system matrices $\widehat{\textbf{A}}, \widehat{\textbf{B}}, \widehat{\textbf{C}}$ , and $\widehat{\textbf{D}}$ using Equation (9).

As shown in McKelvey et al. [23], in the limit of an infinite set of samples the reduction is balanced since the observability and controllability Gramians converge to the diagonal matrix of system singular values (24). From this it can be seen that the model reduction using the finite series is approximately balanced in a similar way to balanced snapshot POD [30].

4.0 Results

All the LFD results required for construction of reduced order models and time-domain solutions for comparison are obtained using the DLR TAU suite of codes [25–27].

4.1. Sampling parameter T

The parameter T is the sampling time interval governing parameter of the bilinear transformation (16). The choice of T is important as it has a major importance in the relationship between a continuous frequency ( $\omega $ ) and a discrete frequency ( $\hat{\omega}$ ). As it is easier to think in terms of continuous reduced frequency, $\omega $ is transformed into a reduced frequency $k$ , where

(38) \begin{align} k{\rm{\;}} = {\rm{\;}}\frac{{{\rm{\omega }}\cdot{\rm{c}}}}{{{U_\infty }}}\end{align}

T has to be chosen such that the continuous reduced frequencies are in the range of interest of the model input. Figure 1 illustrates the impact of increasing T by a factor of 10 starting from T = 0.006. This shows that for T = 0.006 the corresponding frequencies are mainly in the range [1,100], for T = 0.06 in the range [0.1, 10] and for T = 0.6 in the range [0.01,1].

Figure 1. Impact of choice of T on bilinear transformation.

4.2. Inviscid results for NLR7301

The system reduction theory mentioned earlier is first applied to two dimensional flows about an aerofoil in a pitching motion at different frequencies. The motion is sinusoidal about 40% chord with mean incidence ${\alpha _m} = {0.5^0}$ and amplitude ${\alpha _{0{\rm{\;}}}} = {0.15^0}$ . The Mach number ${M_\infty } = 0.68$ . The mesh has 211 points on the aerofoil, 34,784 cells and the mesh near the aerofoil is shown in Fig. 2. The plots in Figs 3 and 4 show the frequency response reconstructed by two ROMs of size 9 and 23, compared to the values obtained with the LFD solver. Figure 3 results were constructed using 32 frequencies, and Fig. 4 results using 64 frequencies. The frequency response constructed from the LFD solutions has a distinct feature for $k = 5$ in magnitude and phase. For $N = 32$ , only the larger ROM closely matches this behaviour, but for $N = 64$ it is captured by both ROMs.

Figure 2. Details of inviscid mesh for NLR7301.

Figure 3. Comparison between the frequency responses of two different ROM, N = 32.

Figure 4. Comparison between the frequency responses of two different ROM, N = 64.

In practice for periodic solutions built with 32 frequencies, even ROM sizes as small as 3 produce good agreement for lift and pitching moment, see for example Fig. 5.

Figure 5. Lift and moment reconstruction, k = 0.5, N = 32, r = [3,15] vs LFD.

The reduced order model built is now used to reconstruct a non-periodic 1-cosine pitching motion. In frequency domain models such as this one, or frequency-based POD, the input data used to build the ROM creation impacts on the frequency responses that the ROM can accurately reproduce. Responses to non-periodic input will only be well predicted if the frequency decomposition of the signal is well represented by the ROM.

The influence of reduced order model size is shown in Fig. 6, where ROM solutions for pitching moment are compared to time-domain simulations. All the ROMs are produced using 128 frequencies. Only the pitching moment is shown because it is more sensitive to ROM size and for this case the lift coefficient shows good agreement for all cases. The pitching moment given by the ROM seems to be accurate for almost all the ROM sizes. Only $r =5$ seems to be too low to capture the entire behaviour.

Figure 6. Pitching moment, Euler vs ROM, k = [0.05.0.5,5], various r, N = 128.

The influence of the number of frequencies used to construct the ROM is then investigated, where $N = \left[ {16;\,32;\,64;\,128} \right]$ . The lift and the pitching moment are reconstructed and compared with Euler time domain simulations in Fig. 7. The reduced frequency chosen is $k = 5$ , since it is the most challenging case, with more dynamics to be captured by the reduced order model. $N = 16$ does not give a good prediction in any case, but from $N = 32$ , the reduced order model gives a satisfactory value. For both lift and moment, a size of 9 ensures a good reconstruction. Moreover, the peak value is well represented. An analysis of the relative error in the maximum value predicted by the ROM and the time at which this maximum occurs, shows errors of about 0.1% for ROMs of size greater than 10, independent of the number of frequencies used in construction (N > 16).

Figure 7. 1-cosine pitching, k = 5, Euler vs ROM based on LFD.

4.3. Inviscid results for the FFAST wing

The FFAST wing model is based upon a typical civil aircraft wing and was developed to provide test cases in the EC FP7 project Future Fast Aeroelastic Simulation Technologies (FFAST) [31]. The surface mesh used for the simulations is shown in Fig. 8 and the freestream conditions are: M = 0.85, reference temperature = 216.65K, reference density = 0.365 kg/m³ and α _m =2degrees. The upper surface C _p and the Mach number for the mean static flow solution is shown in Fig. 9, where the strong shock, at approximately 70% chord, is evident.

Figure 8. Euler mesh, FFAST wing.

Figure 9. Static C _p and Mach number on upper surface of FFAST wing, α_m = 2degrees.

Experience of the authors has shown that rather than creating a global ROM for the whole wing directly, more reliable ROMs that better handle shock motion are obtained by taking a strip approach to its construction. A strip-wise approach allows the behaviour of lift and moments for all individual strips to be represented in the final combined ROM. For a global ROM of the same size, certain strips may have more dominant behaviour and the smaller changes that impact force coefficients on some sections of the wing may not be captured within the ROM. Further it is possible to extract strip wise forces suitable for aeroelastic simulation, which may include beam stick models, from the strip wise model. Run times for global or strip wise ROMs of the same size are identical. Hence results for a strip-based ROM are thus presented here. The aerodynamics of the FFAST wing is therefore described in terms of a set of spanwise strips, see Fig. 10.

Figure 10. Definition of the strip numbers.

The total forces and moments on each strip are calculated using LFD for different frequencies. Surface cells contained entirely within a strip are associated with that strip. Any cell in more than one strip is split into the correct number of sub-cells with each associated with the relevant strip; each part of the split cells maintains the properties of the original cell. Forces and moments are then calculated for each strip by summing the contributions of pressure from all the cells associated with that strip.

The spanwise variation of the aerodynamic coefficients for the test case is shown in Fig. 11. Here the arrow represents the increasing values of reduced frequency, $k$ . ${\mathrm{Re}}\!\left( {{F_Z}} \right)$ shows the typical distribution of the vertical force on a wing, and both values of ${\mathrm{Re}}\!\left( {{F_Z}} \right)$ and ${\mathrm{Im}}({F_Z}$ ) decrease with the reduced frequency. Due to the shock strength near y = 19m, ${\mathrm{Re}}\!\left( {{M_Y}} \right)$ reaches a minimum and ${\mathrm{Im}}\!\left( {{M_Y}} \right)$ reaches a maximum. For selected strips (1, 5, 10, 15, 25) or (2%, 18%, 38%, 58%, 100% of semi span) the magnitude and the phase of ${F_Z}$ is plotted as a function of the reduced frequency so the dependence of these on the strip location can be shown, see Fig. 12. Different trends are observed as far as the magnitude is concerned, especially for Strip 1 that shows one more inflection point than the others.

Figure 11. ${\mathrm{Re}}\!\left( {{F_Z}} \right)$ , ${\mathrm{Im}}\!\left( {{F_Z}} \right)$ , ${\mathrm{Re}}\!\left( {{M_Y}} \right)$ , ${\mathrm{Im}}\!\left( {{M_Y}} \right)$ vs wingspan.

Figure 12. $abs\!\left( {{F_Z}} \right)$ , $ang\!\left( {{F_Z}} \right)$ vs reduced frequency, different strips.

A reference set of data is built with 128 input frequencies. The entire data set is used for comparisons, but only subsets are used to create ROMs. This means the ROMs can be tested outside of the data used in its construction. Models are presented that have been built with 16, 32 and 64 frequencies. After consideration of error levels [32], a reference value of T = 0.02 was selected for the sample spacing.

To be able to have on overview of the trend of the errors depending on the ROM size, the Mean Absolute Percentage Error (MAPE) is calculated. For each ROM size, it actually represents the averaged sum of the relative errors at the 128 reduced frequencies of the reference LFD data set.

(39) \begin{equation}MAPE = {\rm{\;}}\frac{1}{{128}}{\rm{\;}}\mathop \sum \limits_{k = {k_1}}^{k = {k_{128}}} \frac{{\left| {Refererence{{\rm{\;}}_{LFD}} - Predictio{n_{ROM}}} \right|}}{{Refererenc{e_{{\rm{\;}}LFD}}}}{\rm{*}}100\end{equation}

A comparison of the average MAPE against the ROM size, for different numbers of inputs is given in Fig. 13.

Figure 13. Strip average MAPE vs ROM size, ${F_Z}$ magnitude, T = 0.02.

The MAPE decreases as ROM size increases, but then plateaus. Increasing $N$ decreases the MAPE, but at a higher cost of simulations for ROM generation. Hence, $N = 32$ represents a good trade-off since the error reaches 0.1% in average in this case. Further if the MAPE is considered for each individual strip then similar trends are observed. This sensitivity study has shown that choosing $N = 32$ and ${T} = 0.02$ will enable efficient and accurate ROMs to be built.

The ability of the reduced order model to reconstruct the vertical force and pitching moment on the wing, for one particular frequency (k = 0.2) is investigated Figs 14 and 15. In each case, the reference magnitude and phase of ${F_z}$ and ${M_Y}$ are represented, and can be compared to the prediction given by ROMs of different sizes.

Figure 14. F_Z magnitude and phase, with relative error vs ROM size, T = 0.02, N = 32.

Figure 15. M_Y magnitude and phase, with relative error vs ROM size, T = 0.02, N = 32.

Figure 16. Relative error in total F_Z magnitude vs ROM size, T = 0.02, N = 32.

The ROMs of size 3 are not accurate, particularly when reconstructing the phase of the aerodynamic coefficients. The error in this case is around 1%, with a peak close to 5% close to the wing root. ROMs of size 9 give a good prediction, with an error between 0.01% and 1% depending on the cases. However, a much better prediction is achieved for $r = 15$ , with a relative error in each case lower than 0.001%. Increasing the ROM size to 30 does not improve the accuracy of the model.

The total force on the wing is given by the sum of the forces on each slice. For 0.01 < k < 1, the error of prediction of the total force on the wing is around 1% for a ROM of size 3, 0.1% for r = 9, and 0.0001% for $r \ge 15 $ (Fig. 16). As far as high frequencies are concerned for $r \ge 15$ , the relative error does not exceed 0.1%. This increase in error is due to the low number of inputs in this frequency range. However, a maximum error lower than 0.1% guarantees a very good accuracy. Increasing the size of the ROM does not increase the accuracy for $r \ge 15$ . It confirms what has been observed when checking the accuracy of the ROM strip by strip.

4.4. Viscous NASA CRM aircraft results

The test case geometry is the Nasa Common Research Model (CRM), see Fig. 17 and the motion considered is a rigid pitching motion. The mesh used is publicly available and was provided for the fourth AIAA drag prediction workshop. The mesh consists of around 3.7 million points with 100,014 surface nodes. The freestream conditions for the transonic test cases are: Mach Number = 0.86, reference temperature = 228.724K, reference density = 0.4588 kg/m³ and α_m = 2 degrees. In all the results shown, the reference point for the calculation of moments is at middle of the fuselage in the $x$ and $z$ direction, with $y = 0$ .

Figure 17. Viscous mesh for NASA CRM.

Convergence of the residual and lift for the steady state calculation are shown in Fig. 18 and the surface values of y ⁺ and the pressure coefficient are presented in Fig. 19. TAU [21] enables a hybrid-Re treatment of turbulent viscous walls. It enables high values of y ⁺ to be used (up to y ⁺ = 50) for transonic flight conditions. The pressure coefficient shows a standard evolution on the surface for such conditions.

Figure 18. Steady calculation, residuals. Ma = 0.86, α = 2 degrees.

Figure 19. Steady calculation, surface values of y+ and pressure coefficient.

To demonstrate the accuracy of the LFD solver for this geometry, comparisons are made with harmonics constructed using time domain simulations where the incidence varies periodically. TAU calculates harmonics from the time-domain simulation after each period of the sinusoidal motion, but it is necessary to wait a few periods so that the values of the harmonics are converged. The time domain simulations had an amplitude of ${0.1^o}$ and the following reduced frequencies have been considered: 0.1, 0.5, 1, 1.5, 2 and 5. Figures 20 and 21 show results for the magnitude and phase of the aerodynamic coefficients. The LFD results were obtained for $N = 64$ different reduced frequencies, but due to the high computational costs only a few values were calculated for the time domain harmonics. There is good agreement between the two methods for both lift and pitching moment.

Figure 20. LFD vs time domain simulation, lift magnitude and phase.

Figure 21. LFD vs time domain simulation, pitching moment magnitude and phase.

4.4.1 Decomposing the surface mesh into components for ROM generation

To build a low-cost model of the time variation of the aerodynamic coefficients on the CRM geometry, the surface is decomposed into component parts, as shown in Fig. 22. Three different reduced order models are built: one for the fuselage, slicing it in different strips along the $x\;$ axis; a second one to model the forces and moments on the horizontal tail plane and finally one for the wing, slicing it into strips in the y direction. The number of frequencies used in the ROM generation process is specified as $N = 64$ and the number of frequencies used for checking is given by $N = 256.$ The sampling parameter used is $T = 0.02$ .

Figure 22. Strip approach.

Figure 23. Real part of the vertical force and pitching moment along the wing, k = 0.2.

Examples of the outputs obtained from the LFD solver are shown in Figs 23 and 24. Figure 23 shows span wise variation of the real part of the vertical force and pitching moment for strips along the wing at reduced frequency k = 0.2. The distribution of ${\mathrm{Re}}\!\left( {{F_Z}} \right)$ , is typical for a wing, with ${\mathrm{\;Re}}\!\left( {{F_Z}} \right)$ decreasing with the wingspan and the real part of the moment reaches a minimum around $y=24\mathrm{m}$ , where the shock is quite strong. For the fuselage in Fig. 24, ${\mathrm{\;Re}}\!\left( {{F_Z}} \right)$ shows a peak where the wing is (between $x=20$ and $x\;=40\mathrm{m}$ ), whereas the moment behaviour is more complex.

Figure 24. Real part of the vertical force and pitching moment along the fuselage.

4.4.2 Reduced order model results: forces and moments on the fuselage

For each fuselage strip, reduced order models of size $r = 3,\;9$ and $13$ are created using $N = 64\;$ frequencies. The model predictions of the response of the aircraft to a sinusoidal pitch oscillation as a function of frequency are compared to reference results obtained from the LFD solver. The prediction of these ROMs can be compared to the reference solutions, which are the results from the LFD solver at $256\;$ frequencies. The ROMs give results for each strip in the form of predictions of the harmonics of the vertical force. For example, for a strip in the middle of the fuselage results are shown in Fig. 25. This shows that the magnitude of the vertical force on Strip 7 is accurately reproduced by all ROMs and the phase of vertical force is accurate except for the ROM with $r = 3$ . This ROM underestimates the magnitude of the phase angle for low frequencies. The phase is more challenging to model, but from $r = 9$ , the ROM prediction is good. Having obtained predictions for each strip, it is possible to model the behaviour of the whole fuselage in a pitching motion (Fig. 26).

Figure 25. Reconstruction of ${F_Z}$ magnitude and phase, Strip 7 – Fuselage.

Figure 26. Reconstruction of ${F_Z}$ magnitude and phase, Fuselage.

The fuselage is a complex part of the aircraft to model since there are large variations in pressure coefficient distributions between the different strips. However, for both magnitude and phase, small ROMs ( $r = 9$ ) give an accurate prediction of the aerodynamic coefficient.

4.4.3 Reduced order model results: forces and moments on the wing

For the wing a similar approach is taken, with ROMs for each strip created first (see for example the results for Strip 3 in Fig. 27) and then the results from these are combined to give predictions of total force and moment. The predictions of the magnitude and phase of the total aerodynamic coefficients on the wing are shown in Figs 28 and 29.

Figure 27. Reconstruction of ${F_Z}$ magnitude and phase, Strip 3 – Wing.

Figure 28. Reconstruction of ${F_Z}$ magnitude and phase.

Figure 29. Reconstruction of ${M_Y}$ magnitude and phase.

In this case, the model with $ r = 3$ is inaccurate for magnitude and phase of both aerodynamic coefficients. However, the ROMs with $r = 9$ and $r = 13$ both give accurate solutions compared to the reference. Even for the phase of ${M_Y}$ which exhibits complex behaviour, the $r = 13$ results match the entire phase graph and the $r = 9$ ROM only has a small discrepancy in the final peak.

4.4.4 Reduced order model results: total forces and moments on the HTP

The horizontal tail plane is not sliced in different strips, so the reduced order models are built for the whole part and results are shown in Fig. 30. Again, good predictions are obtained except for the smallest ROM with $r = 3$ , which does not capture all the dominant dynamics and hence cannot recreate the high order model results. Similar observations are made for the pitching moment.

Figure 30. Reconstruction of ${F_Z}$ magnitude and phase.

5.0 Discussion of efficiency

The efficiency of reduced order models depends not only on the cost to run the model compared to the high-fidelity model it replaces, but also on the costs of model construction. The reduced order models have negligible cost to run: for a typical unsteady simulation of a full aircraft the ROM is run on a laptop in seconds as compared to needing to use many hours of HPC time. The main efficiency question thus relates to the costs of construction of the ROM compared to running an unsteady simulation. The LFD solver in TAU has been the subject of considerable development to improve its efficiency so that using a single grid ILU, GMRES approach gives significant computational advantages compared to running a periodic time domain solution to produce the periodic solution [27]. For an aerofoil case with 11,800 grid nodes 1 LFD simulation takes the same CPU time as 1.86 time steps of the unsteady time stepping scheme in TAU [27]. Therefore, the costs of construction of a ROM using 128 frequencies would be equivalent to a single time domain simulation with 237 time steps. Thus, one unsteady time domain solution has a similar order of cost to the ROM construction, but the ROM can be used for many simulations. A similar analysis for a swept wing with $1.2 \times {10^6}$ grid nodes shows that 1 LFD simulation takes the same CPU as 4.3 time steps [27]. Thus, the costs of construction of a ROM using 128 frequencies would be equivalent to a single time domain simulation with 550 time steps. Again, one unsteady time domain solution has a similar order of cost to the ROM construction. It should be noted that these are illustrative examples, and the specific costs will depend on the geometry, mesh size and flow conditions since these will all impact convergence and computational time requirements.

To allow accurate modelling of non-dimensional frequencies up to 10 using a time domain ERA it would be expected to use 80–240 timesteps per pulse response. As 2 pulse responses are required for each motion (displacement and velocity) and 4 pulse responses may be required if the code is not linearised, then an equivalent ERA based model would require between 160 and 960 timesteps to build the ROM. This makes the cost of construction of both the developed method and the time domain ERA very similar. A detailed analysis using both methods would be required to decide which were the more efficient. However, it should be noted that whilst TAU is optimised for LFD simulations it has not been optimised to produce the linearised time domain solutions needed for ERA ROM construction. One can see that an optimised preconditioner for linearised time domain solutions would significantly decrease the CPU required to produce the pulse responses.

6.0 Conclusion

The model order reduction method presented here modifies one of the system identification approaches of McKelvey et al. [23] so that the known system matrix ${\textbf{D}}$ is preserved, so as to retain the physics of the term since the original system is known. The ROMs produced with this approach using a small number of frequency simulations show a strong ability to reconstruct the system response for test cases ranging from aerofoils to full aircraft. The main advantage of the current model order reduction approach is that the method does not require the formation and storage of large matrices, such as in POD approaches. In this respect the relationship of the current method to other frequency-based POD methods, is similar to that Eigensystem Realisation Algorithm methods have with time domain snapshot approximate balanced POD schemes. However, unlike balanced POD [8] the method does not allow the reconstruction of adjoint modes. Future research will consider a frequency domain equivalent of balanced POD.

Using a strip wise approach to ROM construction for an aircraft or wing yields a reduced order model that enables loads prediction at negligible computational cost, since strip forces are known. A strip wise approach also gives better models than a global ROM of the same size, since some areas of the wing or fuselage may dominate and the details of the force responses in other areas of the wing may not be captured. A sensitivity analysis has been carried out on the optimal choices for the number of frequencies needed and the ROM size providing the best trade-off between accuracy and efficiency established. This has shown that for the inviscid wing ROMs built with only 32 input frequencies enable the reconstruction of the frequency response with a maximum error of 0.1%, but for the viscous full aircraft 64 input frequencies was needed.