Modelling the unsteady lift of a pitching NACA 0018 aerofoil using state-space neural networks

Abstract

The development of simple, low-order and accurate unsteady aerodynamic models represents a crucial challenge for the design optimisation and control of fluid dynamical systems. In this work, wind tunnel experiments of a pitching NACA 0018 aerofoil conducted at a Reynolds number $Re = 2.8 \times 10^5$ and at different free-stream turbulence intensities are used to identify data-driven nonlinear state-space models relating the time-varying angle of attack of the aerofoil to the lift coefficient. The proposed state-space neural network (SS-NN) modelling technique explores an innovative methodology, which brings the flexibility of artificial neural networks into a classical state-space representation and offers new insights into the construction of reduced-order unsteady aerodynamic models. The work demonstrates that this technique provides accurate predictions of the nonlinear unsteady aerodynamic loads of a pitching aerofoil for a wide variety of angle-of-attack ranges and frequencies of oscillation. Results are compared with a modified version of the Goman–Khrabrov dynamic stall model. It is shown that the SS-NN methodology outperforms the classical semi-empirical dynamic stall models in terms of accuracy, while retaining a fast evaluation time. Additionally, the proposed models are robust to noisy measurements and do not require any pre-processing of the data, thus involving only a limited user interaction. Overall, these features make the SS-NN technique an excellent candidate for the construction of accurate data-driven models from experimental fluid dynamics data, and pave the way for their adoption in applications entailing design optimisation and real-time control of systems involving lift.

1. Introduction

In recent years, the importance of unsteady aerodynamics in a wide range of fluid dynamic applications has led to a substantial increase in research efforts. Typical examples of unsteady flows include wind turbines (Leishman 2002; Sebastian & Lackner 2013; Le Fouest & Mulleners 2022), helicopter rotors (Ganesh et al. 2005; Mulleners, Kindler & Raffel 2012; Tan & Wang 2013), micro-aerial vehicles (Ho et al. 2003; Golubev & Visbal 2012) and bio-inspired aerodynamics (Birch & Dickinson 2001; Dong, Bode-Oke & Li 2018). The need for accurate predictions of the unsteady aerodynamic loads has promoted the development of high-fidelity techniques, such as computational fluid dynamics (CFD) and wind tunnel experiments. These approaches, however, are time intensive and costly, which makes them unsuited whenever fast mathematical models are needed. Indeed, engineering applications involving design optimisation or real-time control require low-order, fast and accurate models.

As far as pitching aerofoils are concerned, the first analytical unsteady aerodynamics models trace back to Wagner (1925), who investigated the lift response due to a step change in angle of attack, and Theodorsen (1935), who developed an unsteady aerodynamics model for a pitching–plunging flat plate. However, these theories, based on potential flows, assume small perturbations and cannot be used for aerofoils at high angles of attack and pitch rates. To overcome these limitations, semi-empirical dynamic stall models have been developed to account for the unsteady nonlinear lift response associated with nonlinear aerodynamic phenomena such as dynamic stall, flow separation and vortex shedding. Examples of dynamic stall models include the Leishman–Beddoes (Leishman & Beddoes 1989), the Goman–Khrabrov (Goman & Khrabrov 1994) and the ONERA (McAlister, Lambert & Petot 1984) models. Although these techniques are of utmost importance within unsteady aerodynamic modelling, they present major shortcomings. In particular, classical dynamic stall models typically have limited accuracy and involve a cumbersome calibration of the model parameters. The latter are generally tuned on static and dynamic experiments, which are therefore necessary for the model estimation.

In recent times, the abundance of numerical and experimental datasets has stimulated the use of data-driven methods. Many different data-driven modelling strategies have been investigated for the identification of unsteady aerodynamics systems, such as linear/linearised state-space models (Brunton, Rowley & Williams 2013; Dawson et al. 2015), polynomial nonlinear state-space models (Decuyper et al. 2018; Siddiqui et al. 2022), block-oriented models (Kou, Zhang & Yin 2016) and the sparse identification of nonlinear dynamics method (Brunton, Proctor & Kutz 2016; Loiseau, Noack & Brunton 2018; Sun et al. 2021). A promising data-driven technique is the state-space neural network (SS-NN) methodology (Schoukens 2021), which exploits the flexibility of artificial neural networks within a classical state-space representation. The state-space formulation proves particularly useful since it facilitates the integration of the model within larger simulation and control frameworks. Previous work (Damiola et al. 2023a) has demonstrated the huge potential of SS-NN models on a test case where numerical CFD data were used for the estimation of the unsteady lift coefficient of a pitching aerofoil.

The present contribution represents a substantial extension of that work, and aims at assessing the capability of state-space neural networks for the construction of reduced-order models from noisy and highly nonlinear wind tunnel experimental data. The use case being considered involves the modelling of the unsteady lift of a pitching NACA 0018 aerofoil for two distinct free-stream turbulence intensities, 0.3 % and 8.2 %. The analysis of different inflow turbulence levels aims at reproducing different scenarios often encountered in real-life applications, where aerofoils experience a variety of free-stream turbulence intensities that can dramatically change their stall characteristics. The SS-NN models are trained using broadband sine-sweep signals, and they are validated on sine-sweep motions and simple harmonic motions performed at different oscillating frequencies and for different angle-of-attack ranges, including pre-stall and post-stall flow conditions. Results are benchmarked with the output of a modified version of the Goman–Khrabrov dynamic stall model (Williams et al. 2017). The Goman–Khrabrov (GK) model is a first-order model that is gaining more and more attention within the fluid dynamics community, especially for closed-loop flow control applications, such as lift hysteresis reduction during pitching manoeuvres (Williams et al. 2015) and gust load alleviation (Williams & King 2018; Sedky, Jones & Lagor 2020). Recent efforts towards the enhancement of the GK model include the work of An et al. (2021), who estimated the instantaneous lift coefficient of a pitching aerofoil by assimilating the predictions of a GK model with limited surface pressure measurements, and the studies of Williams et al. (2019) and De Troyer et al. (2022), who extended the formulation of the GK model to account for active flow control by slot blowing and plasma actuation, respectively. Finally, Ayancik & Mulleners (2022) proposed to revisit and generalise the GK model by replacing its empirical parameters with physics-based time constants.

The remainder of the manuscript is structured as follows. Section 2 details the methodology proposed in this work, which includes a description of the experimental set-up and a thorough explanation of the implementation of the SS-NN model and the GK model. Results are compared and discussed in § 3, and the major outcomes of the work are summarised in § 4.

2. Methodology

The methodology proposed in this work relies on a black-box modelling approach, in which unsteady aerodynamic data are obtained through dedicated wind tunnel experiments. Results are then compared with the predictions of a modified version of the GK dynamic stall model.

2.1. Wind tunnel experiments

Experiments are performed in the open-circuit wind tunnel of the Vrije Universiteit Brussel. This facility is characterised by a rectangular inlet test section with a width of 2 m and a height of 1 m. The test chamber has a length of 12 m and features a divergent ceiling (opening angle equal to 0.29 degrees) which accounts for the development of the boundary layer. Figure 1 depicts the experimental set-up, which features a rectangular wing with NACA 0018 aerofoil shape, chord length of 300 mm and aspect ratio of 1.8. The wing comprises a central 3D printed part and two equal carbon-fibre sections on the sides. A circular rod connects the wing to the sidewall of the wind tunnel, and an actuation mechanism is used to impose a rotation about the pitch axis, which is located at the mid-chord. To minimise three-dimensional effects, two equal endplates of circular shape are positioned at the tip and root of the wing. This set-up is capable of testing a wide range of pitching motions, including simple harmonic signals, sine sweeps and random-phase multisines. The solid blockage ratio is between 4.5 % and 7 % for every imposed motion. The wing features 47 pressure taps, as illustrated in figure 1(b), which are situated at mid-span. Data are acquired using a ZOC33/64Px Scanivalve^® pressure measurement system which samples at 200 Hz. As there is no tap in the rearmost part of the aerofoil, the pressure at the trailing edge is determined using linear extrapolation. The aerodynamic forces exerted on the wing are determined by numerically integrating the readings of the pressure taps.

Figure 1. (a) CAD model of the experimental set-up, (b) graphical representation of the 47 pressure taps located at the mid-span of the wing and (c) square-mesh grid used to increase the free-stream turbulence level (dimensions are in centimetres). Figure taken from Damiola et al. (2023b).

By adjusting the distance between the grid illustrated in figure 1(c) and the wing model, a range of longitudinal free-stream turbulence intensities between 0.3 % and 8.2 % can be achieved. For more details regarding the experimental set-up, and for a comprehensive analysis of the free-stream turbulence characteristics and flow uniformity, the reader is referred to Damiola et al. (2023b), in which dual-component hot-wire measurements were used to evaluate the inflow conditions.

2.2. The aerodynamic characteristics of a NACA 0018 aerofoil at $Re = 2.8 \times 10^{5}$

Before proceeding with the modelling task, it is important to have a good understanding of the system under consideration. For this purpose, the aerodynamic characteristics of the aerofoil are evaluated under quasi-static conditions for the two different free-stream turbulence levels, $I_u=0.3\,\%$ and $I_u=8.2\,\%$. Measurements are acquired by pitching the aerofoil around its mid-chord in a sinusoidal manner at a reduced frequency $\kappa = ({\rm \pi} f c)/U_{\infty } = 0.0006$, which corresponds to an equivalent reduced pitch rate $r_{eq}=\alpha _1 ({\rm \pi} /180) \kappa = 2.8\times 10^{-4}$. This value is sufficiently small to ensure an accurate evaluation of the static aerodynamic loads. Results are obtained by computing the phase average over three subsequent periods. Figure 2 illustrates the quasi-static lift coefficient obtained for the considered NACA 0018 aerofoil at a Reynolds number $Re = 2.8 \times 10^5$. Note that the values shown in the figure represent averages over a sliding window containing 130 samples. The different free-stream turbulence intensity completely alters the static stall characteristics of the aerofoil. While for the clean tunnel ($I_u=0.3\,\%$) an abrupt stall occurs at $21^\circ$, at high turbulence intensity ($I_u=8.2\,\%$) stall is delayed to $23^\circ$ and exhibits a much more gentle behaviour. The sudden occurrence of lift stall at a low level of free-stream turbulence is attributed to the bursting of a leading-edge laminar separation bubble (LSB) on the suction side of the aerofoil, as supported by the experimental results of Gerakopulos, Boutilier & Yarusevych (2010), Mueller-Vahl et al. (2015) and Damiola et al. (2023b). The larger free-stream disturbances also produce an increase in the maximum lift coefficient by approximately 13 %. Another crucial difference is the presence of a sizeable stall hysteresis at low turbulence level in the angle-of-attack range $13<\alpha <22$, which is related to the breakdown of the LSB and indicates the existence of two stable configurations. This flow feature is absent in the highly turbulent case, in which the large free-stream disturbances trigger early laminar–turbulent transition and prevent the formation of the LSB.

Figure 2. Quasi-static lift coefficient of a NACA 0018 profile at $Re = 2.8 \times 10^5$ for two different free-stream turbulence levels. Solid line indicates upstroke motion, dashed line indicates downstroke and the coloured band represents standard deviation. Note the positive deviation from the linear curve at $\alpha = 8^\circ$ for the low $I_u$ case, betraying the presence of a LSB. This is absent in the high $I_u$ case.

2.3. State-space neural network model

Among the many existing data-driven modelling tools available, the state-space representation is selected because it possesses universal approximation properties which make this approach suitable for a large range of engineering applications (Schoukens & Ljung 2019). The classical formulation of a nonlinear state-space model in discrete time is given by

(2.1)\begin{equation} \left.\begin{gathered} x(k+1) = f\left( x(k), u(k) \right) \\ y(k) = g\left( x(k), u(k) \right), \end{gathered}\right\} \end{equation}

where $x(k)$ indicates the states of the system at time $k$, $u(k)$ is the input of the system and $y(k)$ is the output of the system. The relation presented above can be reformulated by explicitly separating a linear part from a nonlinear part. In the context of this work, the latter is described using a single-hidden-layer neural network. The resulting state-space representation is denominated SS-NN, and can be written as follows (Schoukens 2021):

(2.2a)\begin{gather} x(k+1) = \left[A\quad B\right] \begin{bmatrix} x(k) \\ u(k) \end{bmatrix} + W_x \tanh \left( \left[W_{fx} \quad W_{fu}\right]\begin{bmatrix} x(k) \\ u(k) \end{bmatrix} + b_f \right) + b_x \end{gather}

(2.2b)\begin{gather}y(k) = \left[C \quad D\right] \begin{bmatrix} x(k) \\ u(k) \end{bmatrix} + W_y \tanh \left( \left[ W_{gx}\quad W_{gu} \right] \begin{bmatrix} x(k) \\ u(k) \end{bmatrix} + b_g \right) + b_y , \end{gather}

where $A$, $B$, $C$ and $D$ are the linear state-space matrices, $\tanh ({\cdot })$ is the hyperbolic tangent activation function and $W_x$, $W_y$, $W_{fx}$, $W_{fu}$, $W_{gx}$, $W_{gu}$, $b_x$, $b_y$, $b_{f}$ and $b_g$ denote the weights and biases of the neural network. The current study considers a single input single output system, in which the input $u(k)$ consists in the time-dependent angle of attack of the aerofoil, and the output $y(k)$ is the lift coefficient. A time step $\Delta t = 0.005 \ \textrm {s}$ is selected, corresponding to the sampling time of the pressure measurement system. It is important to remark that the time step selected during the training phase becomes an intrinsic property of the model, and thus must be kept unchanged when evaluating the model. Considering a model structure with $n$ states and with two single-hidden-layer neural networks of equal size ($n_x = n_y$ neurons) to model the nonlinear contributions in the state and output equations, the free parameters of the model assume the dimensions $A \in \mathbb {R}^{n \times n}$, $B \in \mathbb {R}^{n \times 1}$, $C \in \mathbb {R}^{1 \times n}$, $D \in \mathbb {R}$, $W_x \in \mathbb {R}^{n \times n_x}$, $W_y \in \mathbb {R}^{1 \times n_y}$, $W_{fx} \in \mathbb {R}^{n_x \times n}$, $W_{fu} \in \mathbb {R}^{n_x \times 1}$, $W_{gx} \in \mathbb {R}^{n_y \times n}$, $W_{gu} \in \mathbb {R}^{n_y \times 1}$, $b_{x} \in \mathbb {R}^{n}$, $b_{y} \in \mathbb {R}$, $b_{f} \in \mathbb {R}^{n_x}$, $b_{g} \in \mathbb {R}^{n_y}$. Figure 3 provides a graphical illustration of the state-space equations reported in (2.2a) and (2.2b). The free parameters of the model, which consist in the linear state-space matrices and in the weights and biases of the neural network, are obtained by minimising the mean-squared-error cost function using the Levenberg–Marquardt algorithm provided in the neural network toolbox of Matlab (Levenberg 1944; Beale, Hagan & Demuth 2018). Following the methodology proposed by Schoukens (2021), the model parameters are initialised with a linear approximation of the nonlinear system. This approach is chosen in the attempt to prevent the nonlinear optimisation algorithm from reaching a local minimum.

Figure 3. Illustration of the SS-NN model structure described in (2.2a) and (2.2b): layers 1 and 2 represent the state equation, layers 3 and 4 represent the output equation. Note that a nonlinear activation function (‘tansig’) is used in the first and third layers, whereas a linear activation function (‘purelin’) is used in the second and fourth layers.

2.4. State-space neural network model training

For an accurate model estimation, it is essential to select a suitable training dataset which effectively covers the parameter space of interest. Within the system identification community, it has been established that broadband signals represent excellent candidates for data-driven modelling applications (Pintelon & Schoukens 2012; Decuyper et al. 2018). The present work makes use of sine-sweep signals, which are excitations rich in information and designed to effectively cover the desired frequency band. From a practical perspective, a further advantage of using sine sweeps is that they closely resemble simple harmonic motions. The latter are frequently encountered in numerous engineering applications, including wind turbines and helicopter aerodynamics. The sine-sweep signals considered in this study excite the frequency range $f = [0- 2]$ Hz, which corresponds to a reduced frequency $\kappa = [0- 0.126]$. These values represent an appropriate and sufficiently wide range for the analysis of unsteady aerodynamics. The mathematical description of the considered sine-sweep motions is given as

(2.3) \begin{equation} \alpha(t) = \left\{\begin{array}{@{}ll} \displaystyle\alpha_0 + \alpha_1 \sin \left[(a_1 t + b_1)t \right] & \text{for \ $0 \leq t < T_0$}\\ \displaystyle\alpha_0 + \alpha_1 \sin \left[(a_2 (t-T_0) + b_2)(t-T_0) \right] & \text{for \ $T_0 \leq t \leq 2 \,T_0$}, \end{array} \right. \end{equation}

where $\alpha _0$ represents the mean angle of attack, $\alpha _1$ is the amplitude of oscillation, $T_0$ is the half-sweep period, $a_1= {\rm \pi}(f_{max}-f_{min})/T_0$, $a_2= {\rm \pi}(f_{min}-f_{max})/T_0$, $b_1= 2 {\rm \pi}f_{min}$ and $b_2= 2 {\rm \pi}f_{max}$. Note that this signal is denominated by linear sine sweep, since the frequency is varied over time in a linear fashion. The complete training dataset is generated through the concatenation of eight sine-sweep motions of equal length covering distinct angle-of-attack ranges, as described in table 1. We choose to concatenate multiple signals in order to obtain one single time series to train the model across the desired parameter space. The relative root-mean-square error, $e_{rms}$, is introduced as a measure of the accuracy of the model, according to the following definition:

(2.4)\begin{equation} e_{rms} = \sum_{i=1}^m e_{{rms}_{i}} \, \frac{N_i}{N} \quad \text{where}\ e_{{rms}_{i}} = \sqrt{\frac{\sum_{k=1}^{N_i} \left( y(k) -\hat{y}(k|\theta) \right)^2} {\sum_{k=1}^{N_i} \left( y(k) - \mathbb{E} (\kern0.7pt y_{i})\right)^2 }} . \end{equation}

In (2.4), $y(k)$ represents the true output of the system and $\hat {y}(k|\theta )$ the output of the SS-NN model, with $\theta$ being the vector containing all model parameters (i.e. the weights and biases of the neural networks, and the coefficients of the linear state-space matrices). Moreover, $N$ denotes the total number of samples, $m$ indicates the number of different parts which compose the data, $N_i$ is the number of samples associated with the $i\textrm {th}$ data part and $\mathbb {E} (\kern0.7pt y_{i})$ is the expected value related to the $i\textrm {th}$ data part.

Table 1. Parameters defining the eight swept sines used for model training.

Two independent SS-NN models are estimated for the two different levels of free-stream turbulence intensity, but the underlying structure of the model is maintained identical. For both datasets, a one-state model ($n=1$) with $n_x = n_y = 20$ neurons defining the neural network represents a good compromise between accuracy and complexity of the model. The choice of a one-state model is also motivated by the desire to provide a fair comparison with the GK model, which is a first-order model. On the other hand, the number of neurons was determined after evaluating multiple models with increasing size of the neural network. Figure 4 illustrates the modelling error on the training data as a function of the number of neurons, for the two considered turbulence intensities. It can be noted that even a small neural network composed of 10 neurons allows for the identification of a very accurate model, whereas a model with zero neurons, which is effectively a linear model, provides a poor description of the highly nonlinear dynamical system under consideration. Moreover, the root-mean-square error remains almost constant when the number of neurons is increased from 10 to 50. It is important to note that the models trained with data at high turbulence intensity converge to a larger error compared with the models obtained at low turbulence level, because of the higher noise level in the measurements. The opposite trend is observed for a number of neurons equal to zero (linear models), where the model obtained at high turbulence intensity performs better than the one obtained at low turbulence level. This result indicates that the nonlinearity of the system is greater at low turbulence intensity, confirming what could be inferred from the quasi-static lift curves, figure 2.

Figure 4. Modelling error on the training data as a function of the number of neurons. Note that the models with zero neurons represent linear models.

Based on this analysis, a one-state model with 20 neurons is chosen for both turbulence levels, and these settings are retained for the remainder of the paper. According to (2.4), the selected models exhibit a root-mean-square error $e_{rms} = 0.125$ for low turbulence intensity, figure 5(a), and $e_{rms} = 0.247$ for high turbulence intensity, figure 5(b). A comparison of the two different scenarios demonstrates that, at elevated free-stream turbulence intensity, the overall error is larger by approximately a factor two, since the incoming flow is characterised by strong velocity fluctuations which result in an extremely high noise floor. Additionally, the performance of the selected models on the training data is graphically illustrated in figure 5. For clarity, the modelling error is visualised as the difference between the true experimental output and the SS-NN model output at all time instants. The error is generally smaller at low angles of attack, whereas it becomes larger at the high angle-of-attack ranges where the unsteadiness and nonlinearity of the flow are greater.

Figure 5. Model training for the two different free-stream turbulence intensities. The top figure illustrates the input signal used for training, i.e. a concatenation of 8 sine sweeps covering the angle-of-attack range $[-5^\circ, 28^\circ ]$. The central figure shows the corresponding output signal used for training, together with the prediction of the SS-NN model. The bottom figure depicts the modelling error defined as $C_{l_{exp}}\small {(k)} - C_{l_{model}}\small {(k)}$; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

2.5. Goman–Khrabrov model

In the present work, the SS-NN technique is compared with a classical dynamic stall model, i.e. the modified version of the GK model developed by Williams et al. (2017). In 1994, Goman & Khrabrov (1994) proposed an empirical low-order state-space model to represent the lift force of a wing in unsteady motion at high angles of attack. The model employs a first-order state equation for the internal dynamic variable $x(t)$, which indicates the degree of flow attachment over the wing. More specifically, the value $x = 1$ corresponds to fully attached flow, $x = 0$ corresponds to fully separated flow and $0< x<1$ indicates partially attached flow states. The evolution of the state variable $x(t)$ is defined by the following differential equation:

(2.5)\begin{equation} \tau_1 \frac{{{\rm d}\kern0.7pt x}(t)}{{\rm d}t} + x(t) = x_0(\alpha(t) - \tau_2\dot{\alpha}(t)). \end{equation}

In (2.5), the function $x_0(\alpha )$ indicates the degree of flow separation when the wing undergoes a quasi-static pitching motion at sufficiently small reduced pitch rate. Moreover, the time constants $\tau _1$ and $\tau _2$ are introduced to scale the dynamic effects of the pitching motion. They represent distinct physical processes:

1. (i) $\tau _1$ reflects the transient aerodynamics of separated flow: any disturbance at fixed $\alpha$ will cause the flow to adjust to the steady state over time in a relaxation process that can be approximated by a first-order differential equation;

2. (ii) $\tau _2$ concerns circulation and boundary-layer convection lags that affect flow separation and reattachment; those lags are proportional to $\dot {\alpha }$, so that these effects can be modelled through an argument shift of the quasi-static separation point, $x_0(\alpha - \tau _2\dot {\alpha }).$

In their original paper, Goman and Khrabrov suggested to use linear cavitation theory to define the lift coefficient in terms of the angle of attack and the point of separation

(2.6)\begin{equation} C_l(\alpha, x) = \frac{\rm \pi}{2} \sin\alpha (1+\sqrt{x})^2. \end{equation}

Later, Williams et al. (2017) generalised the output equation to extend the validity of the model also in the presence of hysteresis in the static lift curve. In the modified formulation, the lift coefficient is given by

(2.7)\begin{equation} C_l(\alpha, x) = 2{\rm \pi}\left[F(\alpha)x(t)+G(\alpha)(1-x(t))\right]. \end{equation}

The functions $F(\alpha )$ and $G(\alpha )$ in (2.7) are defined as

(2.8)\begin{equation} \left.\begin{gathered} F(\alpha) = m_{pre}\alpha,\\ G(\alpha) = m_{post}(\alpha-\alpha_{off}), \end{gathered}\right\} \end{equation}

where $m_{pre}$ is derived from a linear fit of the quasi-static lift curve in the attached flow region (pre-stall), whereas $m_{post}$ and $\alpha _{off}$ are computed through a linear fit in the fully separated flow region (post-stall). Note that, for the considered symmetrical NACA 0018 aerofoil, no offset is required for $F(\alpha )$. Figure 6 graphically illustrates the curve fitting of the lift coefficient data in the pre-stall and post-stall regions. For low turbulence intensity, figure 6(a), linear approximations are computed in the angle-of-attack ranges $\alpha \in [0, 11]$ and $\alpha \in [23, 26]$, whereas for high turbulence intensity, figure 6(b), the linear fitting is performed in the ranges $\alpha \in [0, 11]$ and $\alpha \in [25, 26]$. The obtained parameters are summarised in table 2.

Figure 6. Quasi-static lift coefficient as a function of the angle of attack: solid line indicates pitch-up motion, dashed line indicates pitch-down motion. Linear approximations are shown for pre-stall (blue line) and post-stall (red line); (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

Table 2. Parameters obtained through linear fit of the lift coefficient data for pre-stall and post-stall.

Once the values for $m_{pre}$, $m_{post}$ and $\alpha _{off}$ have been established, the function $x_0(\alpha )$ can be evaluated by inserting the measured quasi-static lift coefficient data into the output equation, according to

(2.9)\begin{equation} x_0(\alpha) = \frac{C_l/2{\rm \pi}-G(\alpha)}{F(\alpha)-G(\alpha)} = \frac{C_l/2{\rm \pi}-m_{post}(\alpha-\alpha_{off})}{m_{pre}\alpha-m_{post}(\alpha-\alpha_{off})}. \end{equation}

The obtained function $x_0(\alpha )$, which represents the steady-state dependence of the separation point on the angle of attack, is illustrated in figure 7. We did not enforce strict bounds on $x_0$, allowing minor deviations from the $[0,1]$ interval. The curve associated with the low free-stream turbulence intensity, figure 7(a), features a large hysteresis loop, which is linked to the bistable behaviour of the lift coefficient for $13^\circ < \alpha < 23^\circ$. By defining two distinct functions for the upstroke and the downstroke motion, the GK model is able to account for the static hysteresis, as suggested by Williams et al. (2017). The switching condition described in table 3 is used to select the appropriate branch of the curve in the evaluation of $x_0$. Conversely, for high free-stream turbulence, as shown in figure 7(b), the absence of hysteresis in the quasi-static lift curve causes the function $x_0$ to be similar during pitch-up and pitch-down, removing the need to distinguish between upstroke and downstroke. It should be noted that $x_0$ is very sensitive to small changes in the parameters defining the linear fits ($m_{pre}, m_{post}, \alpha _{off}$), making it necessary to carefully select the appropriate ranges for the linear approximations.

Figure 7. Quasi-steady point of separation $x_0$ as a function of the angle of attack: solid line indicates pitch-up motion, dashed line indicates pitch-down motion; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

Table 3. Switching condition in the presence of static hysteresis for low free-stream turbulence intensity ($I_u = 0.3\,\%$). Note that $\alpha _{reattach}^{st}=13.1^\circ$ and $\alpha _{stall}^{st}=21.3^\circ$.

In principle, $\tau _1$ and $\tau _2$ can be obtained from flow measurements directly. In practice, however, this is not straightforward, e.g. when only load measurements are available. Therefore, it is customary to determine $\tau _1$ and $\tau _2$ empirically, by numerically minimising the error between the model output and one dynamic experiment (a training dataset). This approach was successfully followed, amongst others, by Williams et al. (2017) and An et al. (2021). On the other hand, Ayancik & Mulleners (2022) recently proposed a different approach in which physics-based information is exploited to evaluate the time constants. This alternative methodology was tested on different aerofoil geometries and flow conditions, yielding positive results. This approach looks particularly valuable whenever only quasi-static measurements are available, and no dynamic experiment can be used for error minimisation.

In the present work, it was decided not to use physics-based time constants, for several reasons. Firstly, the physics-based formulation was derived either for simple harmonic motions or for constant pitch rate motions, whereas in the present study sine-sweep motions are considered. These signals are intrinsically characterised by a continuously varying oscillating frequency, which would require $\tau _2$ to vary in a continuous manner during the sweep. Moreover, this study considers different free-stream turbulence levels which cause dramatic changes in the flow physics over the aerofoil, and thus make it necessary to define different time constants for each turbulence level. For the reasons listed above, it was chosen to evaluate the time constants $\tau _1$ and $\tau _2$ based on a numerical procedure which minimises the root-mean-square error between the model output and a training dataset. This approach enhances the accuracy of the GK model by utilising additional information to fine tune the time constants. The selected training dataset consists of two sine-sweep motions which are also included in the training data of the SS-NN model. More specifically, we choose the sine-sweep signals labelled no. $2$ and no. $6$, according to table 1. Figure 8 illustrates the root-mean-square error on the selected training dataset as a function of the time constants of the GK model, for the two distinct free-stream turbulence intensities. It can be noted that the two cases require different time constants for the minimisation of the cost function. The optimal ($\tau _1$, $\tau _2$) combinations used in the remainder of the paper are listed in table 4.

Figure 8. Contour plots of the root-mean-square error evaluated on the selected training dataset, as a function of the time constants of the GK model. The white cross defines the optimal ($\tau _1, \tau _2$) combination; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

Table 4. Time constants which minimise the root-mean-square error of the GK model on the selected training dataset ($t_c$ represents the convective time $c/U_{\infty }$).

To conclude, the GK model is a simple empirical model capable of reproducing dynamic stall effects on a pitching aerofoil, using physically interpretable parameters. Yet this simplicity limits the range of nonlinear phenomena it can reproduce.

2.6. Behaviour of the SS-NN model in quasi-static conditions

A desirable feature for a dynamical model is its ability to revert to linearised perturbations about quasi-steady conditions for small amplitudes and low frequencies. This property is inherent to the GK model, as the nonlinear function $x_0(\alpha )$ described in (2.9) is constructed using the quasi-static lift coefficient. Indeed, the GK model is a static nonlinear model built upon the quasi-static curve, and it is then augmented with a time delay to account for the system dynamics. In the case of the proposed black-box SS-NN model, which is a dynamic nonlinear model, this property is not strictly enforced through the model structure and the model need not necessarily reduce to a static nonlinear system for low pitch rates. This feature, however, can be implicitly included in the model by providing a sufficiently rich training dataset. It is thus important to analyse the phase space covered by the training data, so that one can identify the bounds within which the model is expected to be accurate. One way to visualise the phase space is to graphically represent the internal state variable of the model ($x$) as a function of the model input, i.e. the angle of attack ($\alpha$), as shown in figure 9. The grey dotted lines in the figure describe the phase-space coverage of the purely dynamic data used for model training, whereas the orange line indicates the prediction of the model when evaluated in quasi-static conditions at reduced frequency $\kappa =0.0006$. It is observed that, at low angles of attack, the phase space is narrow, indicating that the dynamic effects are less pronounced in the attached flow regime, whereas at higher angles of attack the phase space expands. The enlargement of the phase space at high incidence is more apparent at low free-stream turbulence level, figure 9(a). In this scenario, the phase space reveals a central region which is never visited by the training data, reflecting the significant lift hysteresis associated with the considered aerofoil shape and flow conditions. It is noteworthy that the internal state variable predicted by the SS-NN model when evaluated in quasi-static conditions (orange line) consistently remains close to the dynamic training data. Furthermore, it is apparent that a qualitative resemblance exists between the predicted state variable for quasi-static conditions and the $x_0$ function of the GK model presented in figure 7. The similarity is observed for both turbulence levels, and represents an interesting result, as it was not enforced in any way through the model structure. When representing the coverage of the phase space in terms of the output–input relationship (lift coefficient vs angle of attack), it can be noted that the predicted lift coefficient in quasi-static conditions is in good agreement with the true experimental data, as illustrated in figure 10. This result is notable and underlines the model's capability to reproduce the static behaviour, despite being predominantly trained on high-frequency pitching motions.

Figure 9. Representation of the internal state variable of the model as a function of the angle of attack. The grey lines indicate the phase-space coverage of the dynamic dataset used to train the model, whereas the orange line indicates the prediction of the model when evaluated in quasi-static conditions; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

Figure 10. Representation of the lift coefficient as a function of the angle of attack. The prediction of the SS-NN model in quasi-static conditions (orange line) is compared against the true experimental quasi-static result (black line). The grey dotted lines in the background illustrate the dynamic dataset used to train the model; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

3. Results

It has been shown in the previous sections that the aerodynamic characteristics of the selected NACA 0018 aerofoil at $Re = 2.8 \times 10^5$ are completely different at different levels of free-stream turbulence intensity. Therefore, engineering models should be able to correctly reproduce the physics related to the distinct physical processes. In the present paper, two independent models are estimated for the two different free-stream turbulence intensities, both for the SS-NN model and for the GK model. The obtained models are then validated at the considered turbulence levels using different kinds of motion, i.e. sine-sweep signals and simple harmonic motions.

3.1. State-space neural network model validation: sine-sweep motions

The chosen SS-NN model is first tested on sine-sweep motions. The major difference with respect to the sine-sweep signals included in the training dataset is the sweep rate, which is five times larger (from $2.4$ to $12\ \textrm {Hz}\ \min ^{-1}$). The pitch amplitude of the motions is set to $6^\circ$ and, to test the performance of the model in a wide range of angles of attack, two different cases with a distinct mean angle are considered. The selected mean angles, $\alpha _0 = 12^\circ$ and $\alpha _0 = 18^\circ$, are chosen since these two experiments constitute a pre-stall and a post-stall case, respectively. Therefore, they represent excellent validation cases which feature different flow behaviours.

The experiment with mean angle of attack of $12^\circ$, figure 11, reaches a maximum pitch angle well below the static stall angle, and thus is characterised by a moderately separated flow and exhibits an almost linear lift response. No large-scale vortices are shed and the flow is generally well behaved. At low turbulence intensity, figure 11(a), the SS-NN model captures very accurately the unsteady lift generated by the aerofoil. On the other hand, the GK model tends to overestimate the maximum lift coefficient and underestimate the minimum lift coefficient, and also exhibits some deviations from the experimental data at low pitching frequency. At elevated turbulence level, figure 11(b), the time series of the instantaneous lift coefficient presents severe load fluctuations. In this highly turbulent case and with such noisy measurements, the performance of the two models is much more similar, and both approaches can predict the mean value of the time-varying lift coefficient to a good approximation.

Figure 11. Model validation on sine-sweep motion ($\alpha _0 = 12^\circ$). The error is defined as $C_{l_{model}}\small {(k)} - C_{l_{exp}}\small {(k)}$; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

A more challenging scenario is encountered when the mean angle of attack of the sine-sweep motion is increased to $\alpha _0 = 18^\circ$, figure 12. At low turbulence intensity, figure 12(a), experimental data show the abrupt occurrence of stall. The highly nonlinear response of the system increases the complexity of the modelling task and highlights the weaknesses of the GK model. The largest deviations from the experimental values appear at low pitching frequencies, where the GK model predicts a too early stall occurrence. The SS-NN model, on the other hand, follows the experimental data much more accurately and reproduces the stall location to a very good degree of approximation. It is interesting to note that the experimental data in this angle-of-attack range contain important cycle-to-cycle variations. Two clear examples can be seen around $t=7.5$ s and $t=14.1$ s: in these cases, the lift coefficient assumes surprisingly high values, which are approximately 60 % larger compared with the values obtained at the same phase angle in the neighbouring cycles. This odd behaviour only occurs for this angle-of-attack range, since in this case the maximum pitch angle is close to the stall angle. In this case, even a small perturbation in the free-stream conditions is sufficient to modify the stall behaviour and cause the flow to reattach earlier, with a less severe lift reduction. This behaviour is clearly not captured by the models, which are not able to predict such a random phenomenon.

Figure 12. Model validation on sine-sweep motion ($\alpha _0 = 18^\circ$). The error is defined as $C_{l_{model}}\small {(k)} - C_{l_{exp}}\small {(k)}$; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

When considering the experiment performed at high turbulence level, figure 12(b), the turbulent boundary layer on the aerofoil surface prevents the formation and bursting of the leading-edge LSB, producing a more gentle stall behaviour. In this case, the modelling errors of the SS-NN model and the GK model are comparable. Both approaches manage to represent the general evolution of the lift coefficient, but leave out small details which are hard to capture due to the extremely large noise level.

A summary of the modelling error for each case is illustrated in table 5. Overall, it can be noted that, at low free-stream turbulence level, the SS-NN model outperforms the GK model at both angle-of-attack ranges, with a relative root-mean-square error which is approximately half compared with that of the GK approach. Conversely, when the free-stream disturbances are high, the performance of the two techniques becomes more comparable, with the SS-NN model showing an error only a few percentage points lower than that of the GK model.

Table 5. Comparison of the modelling error evaluated on the sine-sweep motions used for validation.

3.2. State-space neural network model validation: simple harmonic motions

In this section, the selected data-driven models are validated on simple harmonic motions which are intrinsically different from the signals that compose the training dataset. These motions hold significant relevance in fluid dynamics due to their frequent utilisation in real-world applications, such as helicopters or wind turbines. Several motions are imposed on the aerofoil in order to test the models under a wide variety of flow conditions, including pre-stall and post-stall. Three distinct frequencies of oscillation are analysed, corresponding to reduced frequencies $\kappa =0.025$, $\kappa =0.063$ and $\kappa =0.100$. It is important to note that the experimental results have been obtained by phase averaging a large number of periods (31 periods for $\kappa =0.025$, 39 periods for $\kappa =0.063$ and 47 periods for $\kappa =0.100$); the figures reported in the remainder of the paper show the mean value of the experimental lift coefficient, together with the associated standard deviation band depicted as a shaded region.

3.2.1. Case 1: $\alpha (t) = 12^\circ + 4^\circ \sin (2 {\rm \pi}f t)$

The first test case under consideration consists of a sinusoidal motion characterised by a small pitch amplitude $\alpha _1 = 4^\circ$. The SS-NN model has been trained on sine sweeps with amplitudes equal to $6^\circ$ and $10^\circ$, and thus this case effectively represents an extrapolation of the results to a smaller amplitude of oscillation. In this example, the maximum angle of attack of the aerofoil is significantly lower than the stall angle in static conditions. Consequently, the flow over the aerofoil exhibits moderate flow separation in the trailing-edge region, and its lift response is characterised by a low degree of nonlinearity.

Considering the experiments conducted at low free-stream turbulence intensity, figure 13, the SS-NN model is able to reproduce very accurately the enlargement of the hysteresis loop and the increase in maximum lift coefficient related to the increased frequency of oscillation. On the other hand, the GK model generally predicts a larger hysteresis loop and overestimates the maximum lift coefficient at high pitching frequencies.

Figure 13. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 0.3\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0017)$, (b) $\kappa = 0.063\ (r_{eq} = 0.0044)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0070)$.

As for the elevated free-stream turbulence level, figure 14, both models are capable of estimating the evolution of the lift coefficient in spite of the high level of flow unsteadiness. Predictions always lie within the experimental standard deviation band, with the SS-NN model giving better results at $\kappa =0.025$ and $\kappa =0.063$, and the GK model providing a slightly better match at the highest reduced frequency $\kappa =0.100$. It is important to underline the robustness and the ability of the SS-NN model to cope with very noisy training data, which were not pre-processed in any way. This aspect is particularly valuable as it suggests a limited user interaction required for the construction of the model.

Figure 14. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 8.2\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0017)$, (b) $\kappa = 0.063\ (r_{eq} = 0.0044)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0070)$.

The relative root-mean-square error is used as a metric to quantitatively compare the accuracy of the SS-NN model and the GK model. Table 6 demonstrates that the SS-NN model dramatically reduces the error at low turbulence levels, whereas results at elevated free-stream turbulence show a more comparable level of accuracy.

Table 6. Modelling error on the simple harmonic motion $\alpha (t) = 12^\circ + 4^\circ \sin (2 {\rm \pi}ft)$.

3.2.2. Case 2: $\alpha (t) = 16^\circ + 8^\circ \sin (2 {\rm \pi}f t)$

The present section considers simple harmonic motions in which the maximum pitch angle, $\alpha _{max}=24^\circ$, is slightly larger than the static stall angle. This case exhibits a qualitatively different flow behaviour between low and high free-stream turbulence levels, as a direct consequence of their distinct static stall characteristics. While at $I_u = 0.3\,\%$ the flow presents an abrupt stall and a strong hysteretic behaviour, at $I_u = 8.2\,\%$ the lift coefficient features a much smaller hysteresis loop. The hysteresis loop is due to the pitching motion of the aerofoil, which tends to delay flow separation in the upstroke motion and promotes late flow reattachment in the downstroke.

As for the low turbulence intensity case, figure 15, the SS-NN model predictions provide excellent agreement with the experimental data. The GK model, conversely, qualitatively captures the shape of the lift hysteresis, but misses critical flow features such as an accurate prediction of the stall angle, the reattachment point and the maximum lift coefficient. As highlighted in table 7, the SS-NN model drastically reduces the root-mean-square error by at least 60 % compared with the classical GK model, reaching an accuracy improvement of 85 % at the smallest reduced frequency $\kappa =0.025$.

Figure 15. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 0.3\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0035)$, (b) $\kappa = 0.063\ (r_{eq} = 0.0088)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0140)$.

Table 7. Modelling error on the simple harmonic motion $\alpha (t) = 16^\circ + 8^\circ \sin (2 {\rm \pi}f t)$.

As for the high turbulence intensity case, figure 16, predictions from the SS-NN model and the GK model exhibit a consistent behaviour in general agreement with the experimental data. The only relevant deviation from the wind tunnel measurements is observed at the highest pitching frequency, where both models estimate a slightly smaller hysteresis loop and underpredict the maximum lift coefficient by $3\,\%$.

Figure 16. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 8.2\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0035)$, (b) $\kappa = 0.063 (r_{eq} = 0.0088)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0140)$.

3.2.3. Case 3: $\alpha (t) = 20^\circ + 8^\circ \sin (2 {\rm \pi}f t)$

Lastly, model validation is conducted on sinusoidal motions reaching angles of attack well above the static stall angle and up to $28^\circ$, where the lift response is highly nonlinear. Therefore, this case is particularly challenging and demanding from a modelling perspective. Similar to the previous case, free-stream turbulence plays a crucial role. The large free-stream disturbances change the lift response of the aerofoil, producing a less abrupt stall and promoting early flow reattachment during the downstroke motion. In this highly nonlinear context, table 8 demonstrates that, for all tested turbulence intensities and pitching frequencies, the root-mean-square error obtained with the SS-NN technique is at least 62 % lower compared with that of the GK model, with performance enhancements up to 82 %. At the lowest free-stream turbulence intensity, figure 17, the SS-NN model provides a significantly better accuracy, especially in the prediction of the stall angle. The GK model badly reproduces this crucial feature and estimates its occurrence at an incidence which is approximately two degrees smaller compared with the actual stall angle. Moreover, the SS-NN model greatly improves the prediction of the shape and size of the hysteresis loop, providing a very good agreement with the experimental results at every angular position. As for the highest free-stream turbulence level, figure 18, the SS-NN approach significantly improves the prediction of the lift coefficient, especially during the downstroke motion, where the GK model overestimates the lift generated by the aerofoil.

Table 8. Modelling error on the simple harmonic motion $\alpha (t) = 20^\circ + 8^\circ \sin (2 {\rm \pi}f t)$.

Figure 17. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 0.3\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0035)$, (b) $\kappa = 0.063\ (r_{eq} = 0.0088)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0140)$.

Figure 18. Model validation on simple harmonic motions characterised by different reduced frequencies ($I_u = 8.2\,\%$). Panels show (a) $\kappa = 0.025 \ (r_{eq} = 0.0035)$, (b) $\kappa = 0.063\ (r_{eq} = 0.0088)$ and (c) $\kappa = 0.100 \ (r_{eq} = 0.0140)$.

To better clarify the flow physics at play, figure 19 provides a spatio-temporal graph depicting the pressure coefficient on the suction side of the aerofoil over four consecutive representative periods obtained at a pitching frequency of 1 Hz ($\kappa = 0.063$). At low free-stream turbulence levels, the passage of a dynamic stall vortex at every pitching cycle can be identified as an oblique line characterised by a significantly higher suction. Conversely, there is no clear evidence of a strong coherent vortical structure at high free-stream turbulence. According to Sharma & Visbal (2019), who conducted a numerical investigation of the flow around four symmetric NACA aerofoils with different thickness-to-chord ratios at $Re = 2 \times 10^5$, the onset of dynamic stall for a thick NACA 0018 profile at low free-stream turbulence might be attributed to the interaction between trailing-edge separation and the LSB, rather than solely being related to the bursting of the LSB, as observed in thinner symmetric aerofoils. Figure 19 also illustrates the lift time history associated with the sinusoidal pitching motion, further highlighting that, in the case of elevated free-stream disturbances, the lift coefficient is generally higher and characterised by large and rapid fluctuations. Despite the high noise level, the SS-NN model consistently provides accurate estimates of the aerodynamic loads.

Figure 19. Spatio-temporal representation of the pressure coefficient on the upper side of the aerofoil over four consecutive representative periods for $\alpha (t)=20^\circ +8^\circ \sin (2{\rm \pi} f t)$ at $f=1$ Hz ($\kappa = 0.063$). The corresponding experimental lift time history (black line) together with the predictions of the SS-NN model (orange line) and the GK model (blue line) are also reported as a function of the non-dimensional time $t/\tau$, with $\tau$ being the period of oscillation; (a) $I_u = 0.3\,\%$ and (b) $I_u = 8.2\,\%$.

Overall, such accurate predictions in this highly nonlinear scenario show the huge potential of the SS-NN technique, which is able to effectively reproduce the hard-to-model nonlinear dynamics of the system at high angles of attack. Although the identification of nonlinear systems is known to be a very difficult task, unsteady aerodynamics models must be able to account for nonlinear effects. Good models must be able to capture flow features such as the stall angle, the reattachment angle and the shape of the hysteresis loop in a precise and reliable manner, since these predictions have a crucial impact in performance estimations and in the assessment of fatigue loads.

4. Conclusions

This work has explored the identification of data-driven nonlinear state-space models in fluid dynamics through the SS-NN modelling technique. This novel methodology, which combines the flexibility of neural networks with a classical state-space representation, provides a new insight into the construction of reduced-order unsteady aerodynamic models. The technique has been applied to experimental wind tunnel data of a pitching NACA 0018 aerofoil at a chord-based Reynolds number of $2.8 \times 10^5$, considering two distinct levels of free-stream turbulence intensity, 0.3 % and 8.2 %. The obtained SS-NN models have a general model structure and relate the time-varying angle of attack of the aerofoil to the lift coefficient. They represent fast and accurate mathematical models which are well suited for applications involving aerodynamic optimisation and control of systems based on lift.

The SS-NN models were trained using a concatenation of eight sine-sweep motions performed at different angle-of-attack ranges. Sine sweeps were found to be suitable for this purpose given their ability to efficiently cover the desired frequency range. Model validation was conducted on distinct sine-sweep motions and on simple harmonic motions performed in a wide range of angles of attack and reduced frequencies. Results were compared with the output of a modified version of the GK dynamic stall model, whose parameters were set using the experimental static lift curve together with two sine-sweep experiments.

Overall, both approaches succeeded in representing the aerodynamic loads on the aerofoil at the different angle-of-attack ranges, pitching frequencies and free-stream turbulence intensities. However, the SS-NN models demonstrated a greater accuracy, especially at low free-stream turbulence level, where static experiments exhibited a large stall hysteresis. Moreover, performance enhancements were particularly evident at high angles of attack, where the SS-NN models significantly outperformed the GK model, providing better predictions in terms of stall angle, reattachment angle and overall shape of the dynamic hysteresis loop. Although the GK model might be satisfactory for applications which do not require high accuracy, its inherent simplicity makes it not suitable whenever load variations must be captured in detail. In these scenarios, the SS-NN technique represents a fast and powerful tool. It is, however, worth mentioning that the presented SS-NN models are constructed using a black-box approach, and thus they trade the (albeit limited) physical interpretability of the GK model for better agreement with the data. Findings also revealed that SS-NN models were robust to noisy measurements and did not require any pre-processing of the data, thus involving a limited user interaction. The entire model estimation only required one set of dynamic experiments and a single training session. Overall, these properties make the SS-NN technique an excellent candidate for the construction of parsimonious and accurate data-driven models, which can lead to substantial performance gains in fluid dynamics systems involving design optimisation and real-time control.