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

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. (Reference Bruschetta28) 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 [Reference Al-Saggaf and Franklin29], 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 [Reference McKelvey, Akçay and Ljung23]. 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. [Reference McKelvey, Akçay and Ljung23]. It should be noted that a key difference between the application here, and that targeted in Ref. (Reference McKelvey, Akçay and Ljung23) 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. (Reference McKelvey, Akçay and Ljung23) 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. (Reference McKelvey, Akçay and Ljung23) 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. [Reference Vakilzadeh, Yaghoubi, McKelvey, Abrahamsson and Ljung24]. 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. [Reference McKelvey, Akçay and Ljung23], 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 [Reference McKelvey, Akçay and Ljung23] 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 [Reference Ma, Ahuja and Rowley12, Reference Juang and Pappa13, Reference Wales, Gaitonde and Jones17]. 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. (Reference McKelvey, Akçay and Ljung23) 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 (Reference McKelvey, Akçay and Ljung23):

(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. (Reference McKelvey, Akçay and Ljung23) 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. [Reference McKelvey, Akçay and Ljung23], 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 [Reference Rowley30].

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 [Reference Thormann and Widhalm27]. 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 [Reference Thormann and Widhalm27]. 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 [Reference Thormann and Widhalm27]. 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.