Hostname: page-component-6898cfcb55-cbk24 Total loading time: 0 Render date: 2026-11-12T02:58:50.565Z Has data issue: false hasContentIssue false
  • English
  • Français

Construction of reduced-order models for fluid flows using deep feedforward neural networks

Published online by Cambridge University Press:  14 June 2019

Hugo F. S. Lui and
William R. WolfOpen the ORCID record for William R. Wolf[Opens in a new window]
Hugo F. S. Lui
Affiliation:
School of Mechanical Engineering, Universidade Estadual de Campinas, Campinas, SP13083-860, Brazil
William R. Wolf*
Affiliation:
School of Mechanical Engineering, Universidade Estadual de Campinas, Campinas, SP13083-860, Brazil
*
†Email address for correspondence: wolf@fem.unicamp.br
Rights & Permissions [Opens in a new window]

Abstract

We present a numerical methodology for construction of reduced-order models (ROMs) of fluid flows through the combination of flow modal decomposition and regression analysis. Spectral proper orthogonal decomposition is applied to reduce the dimensionality of the model and, at the same time, filter the proper orthogonal decomposition temporal modes. The regression step is performed by a deep feedforward neural network (DNN), and the current framework is implemented in a context similar to the sparse identification of nonlinear dynamics algorithm. A discussion on the optimization of the DNN hyperparameters is provided for obtaining the best ROMs and an assessment of these models is presented for a canonical nonlinear oscillator and the compressible flow past a cylinder. Then the method is tested on the reconstruction of a turbulent flow computed by a large eddy simulation of a plunging airfoil under dynamic stall. The reduced-order model is able to capture the dynamics of the leading edge stall vortex and the subsequent trailing edge vortex. For the cases analysed, the numerical framework allows the prediction of the flow field beyond the training window using larger time increments than those employed by the full-order model. We also demonstrate the robustness of the current ROMs constructed via DNNs through a comparison with sparse regression. The DNN approach is able to learn transient features of the flow and presents more accurate and stable long-term predictions compared to sparse regression.

JFM classification

Mathematical Foundations: Computational methodsNonlinear Dynamical Systems: Low-dimensional models

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 872 , 10 August 2019, pp. 963 - 994
Copyright
© 2019 Cambridge University Press 

1 Introduction

Current supercomputers allow the application of high fidelity numerical simulations of turbulent flows over large-scale industrial configurations. The results from these simulations certainly improve the understanding of complex physical phenomena such as mixing enhancement, drag reduction, heat transfer and noise generation, to name a few. Such simulations may lead to discretizations with billions of degrees of freedom in order to resolve the energetically relevant spatial and temporal scales. In these cases, the fluid flow data is generally obtained for long periods using small time steps to compute meaningful converged statistics of the flow with sufficient accuracy.

The analysis of unsteady flows by time-resolved simulations and experiments require the acquisition and treatment of large datasets. In recent years, data-driven algorithms have been developed and applied to perform statistical post-processing of such datasets of unsteady fluid flows allowing the investigation of complex physical mechanisms in turbulent flows and improving their analyses. For example, one can cite techniques of flow modal decomposition such as proper orthogonal decomposition (POD) and its variations (Lumley Reference Lumley1967; Sieber, Paschereit & Oberleithner Reference Sieber, Paschereit and Oberleithner2016; Ribeiro & Wolf Reference Ribeiro and Wolf2017; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018), dynamic mode decomposition (DMD) and variations (Schmid Reference Schmid2010; Tu et al. Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014; Clainche & Vega Reference Clainche and Vega2017), Lagrangian coherent structures (LCS) (Green, Rowley & Haller Reference Green, Rowley and Haller2007; Haller Reference Haller2015; Rockwood, Taira & Green Reference Rockwood, Taira and Green2016) and resolvent analysis (Luhar, Sharma & McKeon Reference Luhar, Sharma and McKeon2014; Sharma et al. Reference Sharma, Moarref, McKeon, Park, Graham and Willis2016), among others. Recently, Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017) published a review of such techniques in the context of fluid flows.

The previous techniques of modal analysis can also be employed for the construction of reduced-order models (ROMs) which are appealing since they can be used in the preliminary stages of design due to their inherent reduction in the computational costs compared to large-scale simulations. Such techniques should also be useful in the context of optimization analyses and studies of flow control (Brunton & Noack Reference Brunton and Noack2015). However, in order to be employed for these applications, ROMs should be able to reproduce the main physical aspects of the full scale models. Several ROM techniques have been discussed in the literature such as Galerkin projection (Rowley, Colonius & Murray Reference Rowley, Colonius and Murray2004), least-squares Petrov–Galerkin projection (Carlberg, Bou-Mosleh & Farhat Reference Carlberg, Bou-Mosleh and Farhat2011), eigensystem realization analysis (Juang & Pappa Reference Juang and Pappa1985) and sparse regression of nonlinear dynamics (Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016). The previous techniques can be used to reduce the original sets of partial differential equations to sets of ordinary differential equations. Sparse regression has also been applied for discovering sets of partial differential equations through spatio-temporal data collection (Rudy et al. Reference Rudy, Brunton, Proctor and Kutz2017). These techniques have mostly been applied for canonical problems and their application to more complex turbulent flows is still a challenging task.

In some cases, ROMs constructed using some of the above techniques may exhibit unstable behaviour (Carlberg et al. Reference Carlberg, Bou-Mosleh and Farhat2011). For example, ROMs developed based on Galerkin projection employ POD to rewrite the Navier–Stokes equations as sets of dynamical systems for the evolution of the POD temporal modes. Due to the basis truncation typical of such methods, numerical instabilities can be created from the imbalance in the turbulent kinetic energy budget in the ROM. This issue may be addressed using turbulence models (Cazemier, Verstappen & Veldman Reference Cazemier, Verstappen and Veldman1998; Östh et al. Reference Östh, Noack, Krajnovi, Barros and Borée2014; Protas, Noack & Östh Reference Protas, Noack and Östh2015), however, this approach can destroy consistency between the original partial differential equations and the ODE system of the ROM. Other alternatives have been also proposed to deal with this problem via minimal rotation of the projected subspace (Balajewicz, Tezaur & Dowell Reference Balajewicz, Tezaur and Dowell2016). This approach is able to account for the contribution of the truncated modes while keeping the consistency between the full- and reduced-order models. Recently, San & Maulik (Reference San and Maulik2018) employed machine learning via neural networks to compute optimal coefficients for an eddy viscosity model which is able to stabilize a POD–Galerkin ROM.

Machine learning is a field that has applications in several areas from data classification to pattern recognition and nonlinear function approximation. Several groups have applied machine learning for investigations concerning fluid flows. Ling & Templeton (Reference Ling and Templeton2015) evaluated different machine learning techniques to classify regions of high uncertainty in RANS calculations. Ling, Kurzawski & Templeton (Reference Ling, Kurzawski and Templeton2016) presented a novel deep neural network (DNN) architecture to improve RANS turbulence models. In a similar context, Wang, Wu & Xiao (Reference Wang, Wu and Xiao2017) developed a machine learning approach based on random forests for predicting discrepancies in Reynolds stresses obtained by RANS calculations. Ströfer et al. (Reference Ströfer, Wu, Xiao and Paterson2019) employed convolutional neural networks, a machine learning technique well suited for image pattern recognition, to identify features in fluid flows.

Machine learning is also a natural candidate for the development of ROMs of large-scale dynamical systems, typical of numerical simulations of unsteady flows. A discussion of the application of deep learning in the context just described was presented by Kutz (Reference Kutz2017). Recently, several authors have proposed algorithms for the prediction of high-dimensional complex dynamical systems using neural networks and their variants (Pan & Duraisamy Reference Pan and Duraisamy2018; Rudy, Kutz & Brunton Reference Rudy, Kutz and Brunton2018; San & Maulik Reference San and Maulik2018; Vlachas et al. Reference Vlachas, Byeon, Wan, Sapsis and Koumoutsakos2018; Wan et al. Reference Wan, Vlachas, Koumoutsakos and Sapsis2018). Rudy et al. (Reference Rudy, Kutz and Brunton2018) presented a methodology that is able to learn the dynamics of a particular system and estimate the noise from measurements using feedforward neural networks (FNN). Pan & Duraisamy (Reference Pan and Duraisamy2018) performed a comparison between FNNs and sparse regression for modelling of dynamical systems and demonstrated the benefits of the former in terms of adaptability. They computed the Frobenius norm of the Jacobian of the neural network as the regularization term of the cost function to improve the accuracy and robustness of a framework. Wan et al. (Reference Wan, Vlachas, Koumoutsakos and Sapsis2018) developed a ROM methodology capable of modelling extreme events occurring in dynamical systems. These authors employed a long short-term memory (LSTM) approach to regularize a recurrent neural network which obtains the complementary dynamics of a nonlinear Galerkin projection of the dynamical system. Vlachas et al. (Reference Vlachas, Byeon, Wan, Sapsis and Koumoutsakos2018) combined a LSTM neural network with a mean stochastic model to propose a data-driven algorithm which has desirable short- and long-term prediction capabilities.

In this work we present a numerical methodology which combines flow modal decomposition via POD and regression analysis using machine learning and FNNs. The framework is implemented in a context similar to that of the sparse identification of nonlinear dynamics (SINDy) algorithm (Brunton et al. Reference Brunton, Proctor and Kutz2016). In order to improve the regression step in the current approach, a spectral proper orthogonal decomposition (SPOD) (Sieber et al. Reference Sieber, Paschereit and Oberleithner2016) is employed to extract a low-dimensional representation of the full-order model (FOM). The SPOD technique is able to filter the temporal modes while preserving the information of the FOM by a redistribution of the energy of the system to higher POD modes. The current method is applied for the construction of ROMs of unsteady compressible flows. First, we test the method for a simple two degree of freedom nonlinear dynamical system. Then, we evaluate the capability of the method to reproduce the compressible flow past a cylinder including its noise generation. In this case, we show that the current DNN approach is also able to reproduce the transient stages of the flow. Finally, the method is tested for a turbulent flow involving dynamic stall of a plunging airfoil. We present a discussion of the optimization of hyperparameters for obtaining the best models for the deep neural networks. For both the cylinder and plunging airfoil cases a comparison between DNN and sparse regression techniques is presented in terms of their predictive capability. The approach presented in this work allows us to predict the flow field beyond the training window and with larger time increments than those used by the FOM, demonstrating the robustness of the current ROMs constructed via DNNs. The current numerical tool for the construction of ROMs via DNNs can be downloaded from http://cces.unicamp.br/software/.

Figure 1.

Schematic of the current method used for construction of ROMs via DNNs and flow modal decomposition.

2 Theoretical and numerical formulation

2.1 Flow modal decomposition

The equations governing an unsteady three-dimensional compressible flow contain partial derivatives with respect to both the spatial coordinates $\boldsymbol{x}=\left[\begin{array}{@{}ccc@{}}x & y & z\end{array}\right]^{\text{T}}$ and time $t$ . Using the method of lines one can first approximate the spatial derivatives producing a system of nonlinear ordinary differential equations (ODEs). In the most general notation, for each mesh point, these ODEs would be expressed in the form

(2.1) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{q}}{\text{d}t}=\boldsymbol{F}(\boldsymbol{q},t),\end{eqnarray}$$

where $\boldsymbol{F}$ is a nonlinear operator, $\boldsymbol{q}=\left[\begin{array}{@{}ccccc@{}}\unicode[STIX]{x1D70C} & \unicode[STIX]{x1D70C}u & \unicode[STIX]{x1D70C}v & \unicode[STIX]{x1D70C}w & p\end{array}\right]^{\text{T}}$ is the vector of flow variables, $\unicode[STIX]{x1D70C}$ is the density, $u$ , $v$ and $w$ are respectively the $x$ -, $y$ - and $z$ -components of the velocity vector and $p$ is the pressure.

Here, we consider that a dynamical system given by (2.1) is solved at each mesh point. To determine the nonlinear operator $\boldsymbol{F}$ from data, we follow the ideas of the SINDy algorithm developed by Brunton et al. (Reference Brunton, Proctor and Kutz2016) with some modifications. A schematic of the current method can be seen in figure 1. First, we collect snapshots of the variables which will be our training data. The data set is then arranged into a matrix $\unicode[STIX]{x1D64C}$ as

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D64C}=\left[\begin{array}{@{}cccc@{}}\boldsymbol{q}(\boldsymbol{x}_{\mathbf{1}},t_{1}) & \boldsymbol{q}(\boldsymbol{x}_{\mathbf{2}},t_{1}) & \cdots \, & \boldsymbol{q}(\boldsymbol{x}_{\boldsymbol{N}_{\boldsymbol{p}}},t_{1})\\ \boldsymbol{ q}(\boldsymbol{x}_{\mathbf{1}},t_{2}) & \boldsymbol{q}(\boldsymbol{x}_{\mathbf{2}},t_{2}) & \cdots \, & \boldsymbol{q}(\boldsymbol{x}_{\boldsymbol{N}_{\boldsymbol{p}}},t_{2})\\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{q}(\boldsymbol{x}_{\mathbf{1}},t_{N_{T}}) & \boldsymbol{q}(\boldsymbol{x}_{\mathbf{2}},t_{N_{T}}) & \cdots \, & \boldsymbol{q}(\boldsymbol{x}_{\boldsymbol{N}_{\boldsymbol{p}}},t_{N_{T}})\end{array}\right],\end{eqnarray}$$

where $N_{p}$ is the number of grid points in the computational domain and $N_{T}$ is the number of snapshots.

Because of the high dimensionality of the input data $\unicode[STIX]{x1D64C}$ , we first reduce the dimension of the dynamical system using the snapshot POD method (Sirovich Reference Sirovich1986). Proper orthogonal decomposition has been widely applied to obtain optimal bases that represent the most energetic content of the system dynamics with as few basis functions as possible (Lumley Reference Lumley1967). The snapshot POD approach starts with a decomposition of the vector quantities $\boldsymbol{q}(\boldsymbol{x},t)$ into the mean flow, $\overline{\boldsymbol{q}}(\boldsymbol{x})$ , and fluctuations, $\boldsymbol{q}^{\prime }(\boldsymbol{x},t)$ . The latter can be further expanded into a combination of spatial modes $\unicode[STIX]{x1D753}_{i}(\boldsymbol{x})$ and their temporal coefficients $\boldsymbol{a}_{i}(t)$ for a defined number of $N$ modes as

(2.3) $$\begin{eqnarray}\boldsymbol{q}(\boldsymbol{x},t)=\overline{\boldsymbol{q}}(\boldsymbol{x})+\boldsymbol{q}^{\prime }(\boldsymbol{x},t)=\overline{\boldsymbol{q}}(\boldsymbol{x})+\mathop{\sum }_{i=1}^{N}\boldsymbol{a}_{i}(t)\unicode[STIX]{x1D753}_{i}(\boldsymbol{x}).\end{eqnarray}$$

To calculate the POD correlation matrix of the data set $\unicode[STIX]{x1D64C}$ some specific norm must be used. For an incompressible flow, a kinetic energy norm provides an optimal result, however, for a compressible flow, other norms can be employed (Rowley et al. Reference Rowley, Colonius and Murray2004). Hence, we define a vector $\unicode[STIX]{x1D6C8}=\left[\begin{array}{@{}ccccc@{}}\unicode[STIX]{x1D702}_{1} & \unicode[STIX]{x1D702}_{2} & \unicode[STIX]{x1D702}_{3} & \unicode[STIX]{x1D702}_{4} & \unicode[STIX]{x1D702}_{5}\end{array}\right]^{\text{T}}$ that determines which norm should be used to compute the correlation between two snapshots. For example, a pressure-based norm uses $\unicode[STIX]{x1D6C8}=\left[\begin{array}{@{}ccccc@{}}0 & 0 & 0 & 0 & 1\end{array}\right]^{\text{T}}$ and a kinetic energy norm uses $\unicode[STIX]{x1D6C8}=\left[\begin{array}{@{}ccccc@{}}0 & 1 & 1 & 1 & 0\end{array}\right]^{\text{T}}$ . The correlation between two snapshots is computed using the $L^{2}$ inner product. Therefore, considering $\boldsymbol{q}^{\prime }(\boldsymbol{x},t_{i})=\boldsymbol{q}_{i}^{\prime }$ and $\boldsymbol{q}^{\prime }(\boldsymbol{x},t_{j})=\boldsymbol{q}_{j}^{\prime }$ , the elements of the correlation matrix $\unicode[STIX]{x1D63E}$ are given by

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D60A}_{ij}=\int _{\unicode[STIX]{x1D6FA}}[\unicode[STIX]{x1D702}_{1}\unicode[STIX]{x1D70C}_{i}^{\prime }\unicode[STIX]{x1D70C}_{j}^{\prime }+\unicode[STIX]{x1D702}_{2}(\unicode[STIX]{x1D70C}u)_{i}^{\prime }(\unicode[STIX]{x1D70C}u)_{j}^{\prime }+\unicode[STIX]{x1D702}_{3}(\unicode[STIX]{x1D70C}v)_{i}^{\prime }(\unicode[STIX]{x1D70C}v)_{j}^{\prime }+\unicode[STIX]{x1D702}_{4}(\unicode[STIX]{x1D70C}w)_{i}^{\prime }(\unicode[STIX]{x1D70C}w)_{j}^{\prime }+\unicode[STIX]{x1D702}_{5}p_{i}^{\prime }p_{j}^{\prime }]\,\text{d}\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the fluid region of interest for the reconstruction and the matrix $\unicode[STIX]{x1D63E}$ is of size $N\times N$ . For the problems studied in this work we employ norms based on kinetic energy and pressure. Despite the changes in the POD modes computed for the different norms, we observed that the ROMs obtained by either kinetic energy or pressure norms were similar. In both cases, stable and accurate models could be reconstructed by the DNNs and further comments are provided in the results section.

In turbulent flows the POD temporal modes may contain contributions from several frequencies including high-frequency noise. For such cases in order to provide a better identification of the individual modes and to smooth out the temporal coefficients, we employ the SPOD described by Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016). The SPOD technique is able to filter the temporal modes while preserving the information of the FOM by a redistribution of the energy of the system to higher POD modes. The technique consists of applying a filter function to the POD correlation matrix which results in a matrix $\tilde{\unicode[STIX]{x1D63E}}$ with elements given as

(2.5) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D60A}}_{ij}=\mathop{\sum }_{k=-N_{f}/2}^{N_{f}/2}g_{k}\unicode[STIX]{x1D60A}_{i+k,j+k}.\end{eqnarray}$$

Here, $g_{k}$ is the filter function and $N_{f}$ is the size of the filter window. We consider a periodic temporal series and assume that the correlation matrix is also periodic (see Ribeiro & Wolf (Reference Ribeiro and Wolf2017) for details). Hence, $N_{f}=f_{snap}/N_{T}$ , where $f_{snap}$ is the number of snapshots used in the SPOD filter. Following this notation, if we apply 50 % of filter to the correlation matrix, it means that we are filtering 50 % of the total number of snapshots. Several filter functions can be applied to the POD correlation matrix and their effects are reported by Ribeiro & Wolf (Reference Ribeiro and Wolf2017). The filter function employed in this work is the box filter represented by $g_{k}=1/(2N_{f}+1)$ .

The temporal coefficients $\boldsymbol{a}_{i}=\left[\begin{array}{@{}ccc@{}}a_{i}(t_{1}) & \cdots \, & a_{i}(t_{N_{T}})\end{array}\right]^{\text{T}}$ and the POD eigenvalues $\unicode[STIX]{x1D706}_{i}$ are determined from the filtered correlation matrix $\tilde{\unicode[STIX]{x1D63E}}$ as

(2.6) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D63E}}\boldsymbol{a}_{i}=\unicode[STIX]{x1D706}_{i}\boldsymbol{a}_{i}.\end{eqnarray}$$

Singular value decomposition (SVD) can be employed to compute the eigenvalues $\unicode[STIX]{x1D706}_{i}$ and eigenvectors $\boldsymbol{a}_{i}$ of $\tilde{\unicode[STIX]{x1D63E}}$ since the matrix is real symmetric positive-definite. The spatial modes are obtained from the projection of the fluctuation quantities onto the temporal coefficients

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D753}_{i}(\boldsymbol{x})=\frac{1}{\unicode[STIX]{x1D706}_{i}}\mathop{\sum }_{j=1}^{N}\boldsymbol{a}_{i}(t_{j})\boldsymbol{q}^{\prime }(\boldsymbol{x},t_{j}).\end{eqnarray}$$

Finally, the temporal coefficients and spatial modes can be stored in matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D731}$ as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D63C}=\left[\begin{array}{@{}c@{}}\boldsymbol{ a}^{\text{T}}(t_{1})\\ \boldsymbol{ a}^{\text{T}}(t_{2})\\ \vdots \\ \boldsymbol{a}^{\text{T}}(t_{N_{T}})\\ \end{array}\right]=\left[\begin{array}{@{}cccc@{}}a_{1}(t_{1}) & a_{2}(t_{1}) & \cdots \, & a_{N}(t_{1})\\ a_{1}(t_{2}) & a_{2}(t_{2}) & \cdots \, & a_{N}(t_{2})\\ \vdots & \vdots & \ddots & \vdots \\ a_{1}(t_{N_{T}}) & a_{2}(t_{N_{T}}) & \cdots \, & a_{N}(t_{N_{T}})\\ \end{array}\right],\end{eqnarray}$$

and

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D731}=\left[\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D719}_{1}(\boldsymbol{x}_{\mathbf{1}}) & \unicode[STIX]{x1D719}_{1}(\boldsymbol{x}_{\mathbf{2}}) & \cdots \, & \unicode[STIX]{x1D719}_{1}(\boldsymbol{x}_{N_{p}})\\ \unicode[STIX]{x1D719}_{2}(\boldsymbol{x}_{\mathbf{1}}) & \unicode[STIX]{x1D719}_{2}(\boldsymbol{x}_{\mathbf{2}}) & \cdots \, & \unicode[STIX]{x1D719}_{2}(\boldsymbol{x}_{N_{p}})\\ \vdots & \vdots & \ddots & \vdots \\ \unicode[STIX]{x1D719}_{N}(\boldsymbol{x}_{\mathbf{1}}) & \unicode[STIX]{x1D719}_{N}(\boldsymbol{x}_{\mathbf{2}}) & \cdots \, & \unicode[STIX]{x1D719}_{N}(\boldsymbol{x}_{N_{p}})\\ \end{array}\right],\end{eqnarray}$$

where the temporal coefficients are the columns of $\unicode[STIX]{x1D63C}$ and the spatial modes are the rows of $\unicode[STIX]{x1D731}$ . Hence, the matrix of fluctuation quantities can be written as

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D64C}\,^{\prime }=\unicode[STIX]{x1D63C}\,\unicode[STIX]{x1D731}.\end{eqnarray}$$

Taking the time derivative of (2.3), we arrive at the following set of equations

(2.11) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{q}(t)}{\text{d}t}=\mathop{\sum }_{i=1}^{N}\unicode[STIX]{x1D753}_{i}(\boldsymbol{x})\frac{\text{d}\boldsymbol{a}_{i}(t)}{\text{d}t}.\end{eqnarray}$$

The last term of (2.11) represents the temporal evolution of coefficients $\boldsymbol{a}_{i}(t)$ associated with the $N$ modes retained in the SPOD modal basis. We can express this system of coupled ODEs as

(2.12) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{a}(t)}{\text{d}t}=\dot{\boldsymbol{a}}(t)=\boldsymbol{F}(\boldsymbol{a}(t)).\end{eqnarray}$$

Next, we compute the derivative $\dot{\boldsymbol{a}}(t)$ numerically using the data $\boldsymbol{a}(t)$ for each temporal mode. Different numerical schemes are tested for the computation of the temporal derivatives and a discussion of the application of explicit schemes and compact schemes will be provided in the results section. For now, let us consider that the numerical scheme employed for the temporal derivatives is a 10th-order accurate compact scheme (Lele Reference Lele1992) which provides high spectral resolution being non-dissipative and low-dispersive. The derivative $\dot{\boldsymbol{a}}(t)$ is then obtained as

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{1}\dot{\boldsymbol{a}}_{i-2}+\unicode[STIX]{x1D6FF}_{2}\dot{\boldsymbol{a}}_{i-1}+\dot{\boldsymbol{a}}_{i}+\unicode[STIX]{x1D6FF}_{2}\dot{\boldsymbol{a}}_{i+1}+\unicode[STIX]{x1D6FF}_{1}\dot{\boldsymbol{a}}_{i+2}=\unicode[STIX]{x1D6FF}_{3}\frac{\boldsymbol{a}_{i+3}-\boldsymbol{a}_{i-3}}{6h}+\unicode[STIX]{x1D6FF}_{4}\frac{\boldsymbol{a}_{i+2}-\boldsymbol{a}_{i-2}}{4h}+\unicode[STIX]{x1D6FF}_{5}\frac{\boldsymbol{a}_{i+1}-\boldsymbol{a}_{i-1}}{2h}.\end{eqnarray}$$

In the equation above, $h$ is the time step and the coefficients of the numerical scheme are set as $\unicode[STIX]{x1D6FF}_{1}=1/20$ , $\unicode[STIX]{x1D6FF}_{2}=1/2$ , $\unicode[STIX]{x1D6FF}_{3}=1/100$ , $\unicode[STIX]{x1D6FF}_{4}=101/150$ and $\unicode[STIX]{x1D6FF}_{5}=17/12$ for a 10th-order discretization. The system of (2.13) written for each temporal mode can be solved as a pentadiagonal linear system for the unknown derivatives $\dot{\boldsymbol{a}}(t)$ . The boundary data points can be computed from the interior points following Olson (Reference Olson2012).

The derivatives $\dot{\boldsymbol{a}}(t)$ are then arranged into a matrix $\dot{\unicode[STIX]{x1D63C}}$

(2.14) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D63C}}=\left[\begin{array}{@{}c@{}}\dot{\boldsymbol{a}}^{\text{T}}(t_{1})\\ \dot{\boldsymbol{a}}^{\text{T}}(t_{2})\\ \vdots \\ \dot{\boldsymbol{a}}^{\text{T}}(t_{N_{T}})\\ \end{array}\right]=\left[\begin{array}{@{}cccc@{}}{\dot{a}}_{1}(t_{1}) & {\dot{a}}_{2}(t_{1}) & \cdots \, & {\dot{a}}_{N}(t_{1})\\ \dot{a_{1}}(t_{2}) & \dot{a_{2}}(t_{2}) & \cdots \, & {\dot{a}}_{N}(t_{2})\\ \vdots & \vdots & \ddots & \vdots \\ {\dot{a}}_{1}(t_{N_{T}}) & {\dot{a}}_{2}(t_{N_{T}}) & \cdots \, & {\dot{a}}_{N}(t_{N_{T}})\\ \end{array}\right].\end{eqnarray}$$

2.2 Regression via DNN

Once the matrix of temporal derivatives is computed, we can set up a regression problem to find the weights $\unicode[STIX]{x1D652}$ and biases $\boldsymbol{b}$ that determine the function $F(\boldsymbol{a}(t))$ presented in (2.12),

(2.15) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D63C}}=\unicode[STIX]{x1D652}\unicode[STIX]{x1D723}(\unicode[STIX]{x1D63C})+\boldsymbol{b},\end{eqnarray}$$

where $\unicode[STIX]{x1D723}(\unicode[STIX]{x1D63C})$ is the matrix of features. In the SINDy algorithm, Brunton et al. (Reference Brunton, Proctor and Kutz2016) suggest that $\unicode[STIX]{x1D723}(\unicode[STIX]{x1D63C})$ may consist of constant, polynomial, exponential and trigonometric functions. However, in many cases it is difficult to know what set of features should be extracted from the data. Hence, we use machine learning to circumvent the problem of finding the functions which represent the dynamics of the problem. Therefore, the methodology can discover not only the weights $\unicode[STIX]{x1D652}$ and biases $\boldsymbol{b}$ but also the features $\unicode[STIX]{x1D723}(\unicode[STIX]{x1D63C})$ . The ‘learned’ features often result in a better performance when compared to those obtained using ‘engineered’ features. A learning algorithm can find a proper set of features in minutes or hours, depending on the task complexity. On the other hand, manually engineered features would require a great amount of human time and effort for complex tasks (Goodfellow, Bengio & Courville Reference Goodfellow, Bengio and Courville2016).

Figure 2.

An example of a deep feedforward network with $N$ inputs, several hidden layers (HL), and one output layer with $N$ outputs.

Deep learning methods are feature learning algorithms that can find a proper set of features using multiple layers, from higher layer features defined in terms of lower layer features. Automatically learning features at multiple processing layers allows the learning of complex functions through mapping the input to the output directly from a given data (Bengio Reference Bengio2009). In the present work, a DNN is used to learn the weights $\unicode[STIX]{x1D652}$ , biases $\boldsymbol{b}$ and features $\unicode[STIX]{x1D723}(\unicode[STIX]{x1D63C})$ of the dynamical systems investigated. Figure 2 shows a sample DNN architecture where the input $\unicode[STIX]{x1D653}$ of the DNN is the matrix $\unicode[STIX]{x1D63C}$ and the target $\unicode[STIX]{x1D654}$ is the matrix $\dot{\unicode[STIX]{x1D63C}}$ . The DNN calculation procedure is presented in the following algorithm form 1. In the current work, the open source machine learning framework Tensorflow (Abadi et al. Reference Abadi2015) is used for training the DNN.

Once we have learned the parameters $\unicode[STIX]{x1D652}^{[l]}$ and $\boldsymbol{b}^{[l]}$ of our model (2.15), we can use them to predict the temporal coefficients $a_{i}$ given the initial conditions $\boldsymbol{a}(t_{1})$ . The system of coupled ODEs (2.12) can now be given by

(2.16) $$\begin{eqnarray}\frac{\text{d}\hat{\boldsymbol{a}}(t)}{\text{d}t}=\unicode[STIX]{x1D652}^{[L]}\unicode[STIX]{x1D642}^{[L-1]}+\boldsymbol{b}^{[L]}=\boldsymbol{F}(\boldsymbol{a}(t)),\end{eqnarray}$$

where $L$ is the index of the last layer, $\unicode[STIX]{x1D642}^{[L-1]}$ is the activation function matrix of the penultimate layer, $\unicode[STIX]{x1D652}^{[L]}$ is the weight matrix of the last layer, and $\boldsymbol{b}^{[L]}$ is the bias vector of the last layer.

This system of coupled ODEs is integrated using an explicit five stage 4th-order Runge–Kutta scheme derived by Kennedy, Carpenter & Lewis (Reference Kennedy, Carpenter and Lewis1999). As we can see in algorithm 1, $\unicode[STIX]{x1D642}^{[L-1]}$ depends on the weights $\unicode[STIX]{x1D652}$ and biases $\boldsymbol{b}$ from previous layers. Thus for each stage of the 4th-order Runge–Kutta, we need to perform forward propagation to obtain $\boldsymbol{F}(\boldsymbol{a}(t))$ . As we have the temporal coefficients $\boldsymbol{a}(t)$ , one can reconstruct the flow field using (2.3). However, we are interested in using a ROM in circumstances other than simply reproducing the training data. The approach presented in this work allows us to predict the flow field beyond the training window because $\overline{\boldsymbol{q}}(\boldsymbol{x})$ and $\unicode[STIX]{x1D753}(\boldsymbol{x})$ depend only on the spatial coordinates $\boldsymbol{x}$ and they are calculated using only the training data.

2.3 Data-driven ROMs

As previously discussed in § 1, the construction of ROMs is an active area of research and different methodologies for the development of ROMs are available. For example, POD–Galerkin based methods are directly related to the physics of the problem through the projection of the Navier–Stokes equations into a system of ordinary differential equations. Galerkin projection methods require the treatment of the linear and nonlinear spatial terms appearing in the FOM. Another issue with these methods relates to their expensive application for dynamical systems with strong nonlinearities where one should employ, for example, hyper-reduction techniques (Chaturantabut & Sorensen Reference Chaturantabut and Sorensen2010; Carlberg et al. Reference Carlberg, Bou-Mosleh and Farhat2011; Zimmermann & Willcox Reference Zimmermann and Willcox2016).

In the current approach, we employ DNNs for the regression step in a context similar to the SINDy algorithm. Both the DNN and the original SINDy approaches are data-driven methods and, instead of having a direct connection with the physics of the problem, these techniques learn from data. An advantage of these methods is that the nonlinearities of the problem are considered in the temporal derivatives of the primitive variables, which are obtained from the FOM. Therefore, neither the SINDy nor the DNN method require the treatment of spatial derivatives.

In general, the resulting ROMs constructed via DNNs are not (directly) physically interpretable due to the nonlinearity of the matrix of features and their weights and biases, which are not sparse. We hope that interpretability of complex models can be improved by new topics of research such as ‘accountable machine learning’ (Navarro et al. Reference Navarro, Navarro, Rocha and Dahab2018). The physics of systems with complicated models can also be interpreted by the structure of the state space, e.g. fixed points, Lyapunov exponents, fractal dimensions, bifurcations and DNN models can also be useful for computing such measures. Learning from models simply uses the fact that many realizations of the physical system can be computed very quickly to explore parts of the state space which are not contained in the training data. In some cases, such as the flow past a cylinder, models constructed using sparse regression could directly lead to physically interpretable results (Brunton et al. Reference Brunton, Proctor and Kutz2016). As expected, the application of DNNs for the regression step adds a penalty cost compared to sparse regression. However it is shown in this work that, for the cases analysed, models obtained using neural networks present better long-term predictive capabilities compared to sparse regression. Both these issues will be discussed in the results section.

3 Hyperparameter optimization

The performance of the DNN described in algorithm 1 depends dramatically on the selection of hyperparameters such as the network depth, $L$ , number of hidden units for each layer, $n^{[l]}$ , regularization parameter, $\unicode[STIX]{x1D706}$ , and learning rate $\unicode[STIX]{x1D6FC}$ . Indeed, finding an optimal set of hyperparameters which minimizes the loss function over a hyperparameter space is a challenging task given the substantial number of free parameters involved.

The manual search, grid search, random search (Bergstra & Bengio Reference Bergstra and Bengio2012) and Bayesian optimization (Brochu, Cora & Freitas Reference Brochu, Cora and Freitas2010) are the most widely used procedures for the hyperparameter optimization. The manual search consists of a direct human trial and error procedure in the search for an optimal configuration of hyperparameters. This procedure is entirely based on prior experience of the user and there is a high probability that an optimal set of hyperparameters is not found. However, a manual search is still useful if the effect of a specific hyperparameter on the model performance can be monitored ‘on the fly’. In the grid search method, several combinations of hyperparameter values are tested in a range evenly spaced. This method is extremely time-consuming because the number of trials increases exponentially with the number of hyperparameters. In a random search, one randomly selects each hyperparameter from a defined range and evaluates the model performance. It is time-consuming when a high-dimensional hyperparameter space is analysed, however, Bergstra & Bengio (Reference Bergstra and Bengio2012) empirically show that a random search outperforms a grid search for the hyperparameter optimization both in terms of computational time and model performance.

One of the recent strategies to find an optimal set of hyperparameters is the Bayesian optimization. It is a technique that involves constructing a probabilistic surrogate model to the data in order to determine the most promising hyperparameters to evaluate. Snoek, Larochelle & Adams (Reference Snoek, Larochelle and Adams2012) showed that Bayesian optimization was able to find optimal hyperparameters for a three-layer convolutional neural network considerably faster than previous approaches and outperformed the state of the art performance at selecting the set of hyperparameters on the CIFAR-10 data set (Krizhevsky Reference Krizhevsky2009).

In this work we use two hyperparameter optimization strategies: random search and Bayesian optimization. For random search, the model generation procedure is presented in the following algorithm form 2. Likewise, Bayesian optimization has the same inputs as random search. To report the performance of each model from a set of candidates, we compute the mean absolute error (MAE) over the training data. The candidate models with lower MAE values are most likely those which will provide the ‘best’ models. It is possible, however, that the model with the lowest MAE suffers from overfitting.

The current metric does not assess the generalization of the model. However, it is the only one available since we cannot split our data into training and validation sets. The validation data set should provide an unbiased evaluation of the ROM. In our case, we employ the temporal coefficients of the POD modes for the training stage of the model. These modes are computed using a correlation of different snapshots and if we construct the ROM using POD modes obtained with data including the validation set, our model would use a biased set for the training stage. This occurs because the POD correlation matrix would be computed for snapshots of both the training and validation sets. In order to overcome this issue, we reconstruct the flow field of the candidate models with lowest MAEs and then compute $E_{r}=\Vert \boldsymbol{q}_{FOM}-\boldsymbol{q}_{ROM}\Vert$ over the validation set to select the best candidate model. This procedure could even be improved if $E_{r}$ were computed for all models. However, the computational cost would considerably increase if such metric were employed.

Table 1.

Effects of hyperparameters on model performance.

An alternative to the current MAE procedure is to use Akaike’s information criterion (AIC) (Akaike Reference Akaike1973) or the Bayesian information criterion (BIC) (Schwarz Reference Schwarz1978) as the model selection criteria. These methods try to balance the quality of the fit and model complexity and the main advantage is that there is no need for a validation set. The downside is that AIC and BIC impose a penalty for model complexity which is related to the number of parameters. This can be a problem when dealing with DNNs due to their large number of parameters.

In our experiments, we noticed that a few of the hyperparameters employed in algorithm 1 need to be tuned for the best performance of the DNN. These are listed in algorithm 2. The remaining hyperparameters are determined according to the following procedures in order to reduce the hyperparameter search space. For instance, hyperparameters related to the Adam optimization, such as the exponential decay rates $\unicode[STIX]{x1D6FD}_{1}$ and $\unicode[STIX]{x1D6FD}_{2}$ , and the constant for numerical stability $\unicode[STIX]{x1D716}$ , are set as described by Kingma & Ba (Reference Kingma and Ba2014). The learning rate $\unicode[STIX]{x1D6FC}$ and the number of iterations $n_{iter}$ are chosen via manual search. Further details regarding this process can be found in Goodfellow et al. (Reference Goodfellow, Bengio and Courville2016, p. 295). As can be observed in the results section, we use the same values of the previous hyperparameters for all flow simulations studied in this work and these values should serve as references for other flows. Following Goodfellow et al. (Reference Goodfellow, Bengio and Courville2016), we provide table 1 with information of the hyperparameters presented in algorithm 2 based on their effects on model performance. This table also presents a similar discussion for the hyperparameters chosen via manual search. We expect that it serves as a guideline for choosing the range of values to be explored for the hyperparameters presented in algorithm 2.

Here, we consider that the activation function is also a hyperparameter since it plays a key role in the DNN performance. There are several available activation functions such as sigmoid, hyperbolic tangent (tanh), rectified linear unit (ReLU) (Nair & Hinton Reference Nair and Hinton2010), exponential linear unit (ELU) (Clevert, Unterthiner & Hochreiter Reference Clevert, Unterthiner and Hochreiter2015), to name a few. For general regression problems, tanh and ELU are the most popular functions since they possess nonlinear properties and are continuously differentiable. Clevert et al. (Reference Clevert, Unterthiner and Hochreiter2015) show that the ELU function reduces the vanishing gradient effect since its positive part returns the identity. Thus, in the positive part, the derivative is unitary and it is not contractive. On the other hand, tanh is contractive almost everywhere. Furthermore, in our experiments the ELU function provided results with fewer iterations than a corresponding tanh network and, therefore, we use ELU as the activation function.

4 Results

In this section, the present method is first tested in the reconstruction of the dynamics of a damped cubic oscillator. Then, we evaluate the capability of the DNN-ROMs to reproduce the dynamics of a compressible flow past a cylinder including its noise generation. A comparison between the current DNN technique against sparse regression is also presented for the transient regime of an incompressible flow past a cylinder. Finally, the method is employed to create a ROM of a turbulent flow involving the dynamic stall of a plunging airfoil. Again, the DNN solution is compared against that obtained by sparse regression. In each example, we show the ability and limitations of the methods to identify the dynamics of the different nonlinear systems comparing the ROM and FOM solutions.

4.1 Nonlinear oscillator

For this first example, we consider a canonical problem in system identification (Brunton et al. Reference Brunton, Proctor and Kutz2016; Rudy et al. Reference Rudy, Kutz and Brunton2018): the two-dimensional nonlinear oscillator. The system dynamics are given by

(4.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left[\begin{array}{@{}c@{}}x_{1}\\ x_{2}\end{array}\right]=\left[\begin{array}{@{}cc@{}}-0.1 & 2\\ -2 & -0.1\end{array}\right]\left[\begin{array}{@{}c@{}}x_{1}^{3}\\[3.0pt] x_{2}^{3}\end{array}\right].\end{eqnarray}$$

with initial conditions $\left[\begin{array}{@{}cc@{}}x_{1} & x_{2}\end{array}\right]^{\text{T}}=\left[\begin{array}{@{}cc@{}}2 & 0\end{array}\right]^{\text{T}}$ .

We generate 4000 snapshots from $t=0$ to $t=40$ by integrating equation (4.1) using an explicit five stage 4th-order Runge–Kutta scheme (Kennedy et al. Reference Kennedy, Carpenter and Lewis1999) with a time step of $h=0.01$ . The training window spans the period $0\leqslant t\leqslant 10$ and the remaining data is used as the test set. The system reconstruction is obtained following the procedures defined in § 2. For this first example, it is not necessary to use POD due to the low dimensionality of the system. The parameters employed in algorithm 2 are presented in table 2 and the best set of hyperparameters obtained via random search is listed in table 3. Figure 3 presents the solutions obtained by the FOM (true model) and ROM for the damped cubic oscillator. Results show that the proposed algorithm accurately reproduces the system dynamics during the training window and beyond for the test set.

Figure 3.

Comparison between true model (solid black line) and ROM prediction (dashed cyan line) for nonlinear oscillator.

Table 2.

Model generation parameters for ROM of nonlinear oscillator.

Table 3.

Best set of hyperparameters for ROM of nonlinear oscillator.

4.2 Construction of ROMs for fluid flow simulations

For the following investigations of the flow past a cylinder and the dynamic stall of a plunging SD7003 airfoil, the system dynamics are modelled by the compressible Navier–Stokes equations in two dimensions and three dimensions, respectively. To simulate the airfoil undergoing a prescribed motion, the equations are solved in a non-inertial frame. In this form, source terms emerge from the grid curvature and frame movement. Here, all terms are solved in the full contravariant form to allow the use of a curvilinear coordinate system $\{\unicode[STIX]{x1D709}^{1},\unicode[STIX]{x1D709}^{2},\unicode[STIX]{x1D709}^{3}\}$ . For a frame of reference with varying linear velocity, the equations reduce to

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\sqrt{g}\unicode[STIX]{x1D70C})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{i}}(\sqrt{g}\unicode[STIX]{x1D70C}u^{i})=0, & \displaystyle\end{eqnarray}$$
(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\sqrt{g}\unicode[STIX]{x1D70C}u^{i})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{j}}[\sqrt{g}(\unicode[STIX]{x1D70C}u^{i}u^{j}-\unicode[STIX]{x1D70F}^{ij}+\unicode[STIX]{x1D628}^{ij}p)]+\left\{\begin{array}{@{}c@{}}i\\ jk\end{array}\right\}\sqrt{g}(\unicode[STIX]{x1D70C}u^{k}u^{j}+\unicode[STIX]{x1D628}^{jk}p-\unicode[STIX]{x1D70F}^{kj})=\sqrt{g}\unicode[STIX]{x1D70C}\ddot{h}^{i}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\sqrt{g}E)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{j}}\left\{\sqrt{g}\left[(E+p)u^{j}-\unicode[STIX]{x1D70F}^{ij}\unicode[STIX]{x1D628}_{ik}u^{k}-\frac{\unicode[STIX]{x1D707}}{Re\,Pr}\unicode[STIX]{x1D628}^{ij}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{i}}\right]\right\}=\unicode[STIX]{x1D70C}\sqrt{g}(h^{j}+u^{j})\unicode[STIX]{x1D628}_{ji}\ddot{h}^{i}.\end{eqnarray}$$

The set of equations above represent the continuity, momentum and energy equations. In order to close the system of equations the following relations are employed:

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle E=\frac{p}{\unicode[STIX]{x1D6FE}-1}+\frac{1}{2}\unicode[STIX]{x1D70C}u^{i}\unicode[STIX]{x1D628}_{ij}u^{j}+\frac{1}{2}\unicode[STIX]{x1D70C}{\dot{h}}^{i}\unicode[STIX]{x1D628}_{ij}{\dot{h}}^{i}, & \displaystyle\end{eqnarray}$$
(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}^{ij}=\frac{\unicode[STIX]{x1D707}}{Re}\left(\unicode[STIX]{x1D628}^{jk}u_{~|k}^{i}+\unicode[STIX]{x1D628}^{ik}u_{~|k}^{j}-\frac{2}{3}\unicode[STIX]{x1D628}^{ij}u_{~|k}^{k}\right), & \displaystyle\end{eqnarray}$$

and

(4.7) $$\begin{eqnarray}h=h_{o}\sin (kt).\end{eqnarray}$$

Here, $\unicode[STIX]{x1D70C}$ represents the density, $u^{i}$ the $i$ th component of the contravariant velocity vector and $p$ is the pressure. The term $h$ is the frame position (cross-stream motion of the plunging airfoil), $E$ is the total energy, $\unicode[STIX]{x1D707}$ is the dynamic viscosity and $T$ is the temperature. The dots represent temporal derivatives of the frame position, i.e., frame velocity and acceleration. The aforementioned terms are non-dimensionalised by freestream quantities such as density $\unicode[STIX]{x1D70C}_{\infty }$ and speed of sound $c_{\infty }$ . The length scales are made dimensionless by the cylinder diameter or airfoil chord. In the above equations, $\unicode[STIX]{x1D628}_{ij}$ and $\unicode[STIX]{x1D628}^{ij}$ are the covariant and contravariant metric tensors, $\sqrt{g}$ is the Jacobian of the covariant metric tensor, and $\big\{\begin{array}{@{}c@{}}i\\ jk\end{array}\big\}$ represents the Christoffel symbols of the second kind. Further details regarding the present formulation can be found in Aris (Reference Aris1989).

The numerical scheme for the spatial discretization of the equations is a high-resolution sixth-order accurate compact scheme (Nagarajan, Lele & Ferziger Reference Nagarajan, Lele and Ferziger2003) implemented on a staggered grid. A sixth-order compact interpolation scheme is also employed (Lele Reference Lele1992) to obtain fluid properties on the different nodes of the staggered configuration. Due to the non-dissipative characteristics of compact finite-difference schemes, numerical instabilities may arise from mesh non-uniformities, interpolations and boundary conditions and, hence, they have to be filtered to preserve stability of the numerical schemes. The high wavenumber compact filter presented by Lele (Reference Lele1992) is applied in flow regions far away from solid boundaries at prescribed time intervals in order to control numerical instabilities.

The time integration of the fluid equations is carried out by the fully implicit second-order scheme of Beam & Warming (Reference Beam and Warming1978) in the near-wall region in order to overcome the time step restriction due to the usual near-wall fine-grid numerical stiffness. A third-order Runge–Kutta scheme is used for time advancement of the equations in flow regions far away from solid boundaries. No-slip adiabatic wall boundary conditions are applied along the solid surfaces and characteristic plus sponge conditions are applied in the far field. For the study of dynamic stall, periodic boundary conditions are applied along the spanwise direction of the airfoil. In this case, we employ an implicit large eddy simulation (LES) to study the flow physics of deep dynamic stall over a plunging SD7003 airfoil profile, i.e., no subgrid scale model is used. The numerical tool has been previously validated for several simulations of unsteady compressible flows (Wolf Reference Wolf2011; Wolf, Azevedo & Lele Reference Wolf, Azevedo and Lele2012a ; Wolf et al. Reference Wolf, Lele, Jothiprasard and Cheung2012b ; D’Lelis et al. Reference D’Lelis, Wolf, Yeh and Taira2019).

4.2.1 Flow past a cylinder

In this case, the FOM is obtained by solving the compressible Navier–Stokes equations as detailed in § 4.2. The numerical simulations are conducted for Reynolds and Mach numbers $Re=150$ and $M=0.4$ , respectively. These dimensionless parameters are computed based on freestream quantities. The grid configuration consists of a body-fitted O-grid with $421\times 751$ points in the streamwise and wall-normal directions, respectively.

The flow is recorded for 1120 snapshots with dimensionless time steps of $h=0.05$ . The snapshots are collected after an initial transient period of the simulation is discarded. The ROM is obtained following the procedure described in § 2 using the pressure norm for the POD correlation matrix. It is important to mention that the use of a kinetic energy norm produced similar results for this case. The training data comprises the first 280 snapshots of the FOM and the remaining data is employed as the test set. The set of hyperparameters employed in algorithm 2 is listed in table 4 and the best set of these parameters obtained via random search is presented in table 5. For this case, Bayesian optimization was also able to produce accurate models but at a higher computational cost compared to random search. One can see that the best architecture of the neural network for this case has 10 layers. However, we also found several other stable and accurate models with six layers, for example. The fact that the best model was found with 10 layers is a coincidence.

Table 4.

Model generation parameters for ROM of compressible flow past a cylinder.

Table 5.

Best set of hyperparameters for ROM of compressible flow past a cylinder.

Figure 4.

Contours of $\boldsymbol{z}$ -vorticity,  $t=410$ .

Figure 5.

Contours of divergence of velocity, $t=410$ .

Figures 4 and 5 show contours of $z$ -vorticity and divergence of velocity, respectively, along the cylinder and wake regions at time $t=410$ , which is beyond the training window. The snapshots allow a comparison of the results between the FOM and ROM using 2 and 10 POD modes out of 280 modes. Hence, the flow is reconstructed using 0.7 % and 3.5 % of the total information available from the FOM. Although the current simulation is performed for a compressible flow at $M=0.4$ and $Re=150$ , we could verify that the POD spatial eigenfunctions are almost identical to those from Noack et al. (Reference Noack, Afanasiev, Morzynski, Tadmor and Thiele2003), which were obtained for an incompressible flow and $Re=100$ . A video comparing the FOM and ROM solutions is provided in the supplementary material (movie 1) available at https://doi.org/10.1017/jfm.2019.358. Reconstruction of the individual flow variables with 2 POD modes could recover between 50 % and 80 % of the total modal energy, depending on the variable. For example, density is reconstructed using 50 % of the total energy of the system dynamics while $y$ -momentum is reconstructed with 80 % of the total modal energy. Here we use the term ‘modal energy’ to refer to the ratio between the sum of $N$ POD eigenvalues used in a particular reconstruction over the entire range $N_{T}$ of eigenvalues available, $\sum _{i=1}^{N}\unicode[STIX]{x1D706}_{i}/\sum _{i=1}^{N_{T}}\unicode[STIX]{x1D706}_{i}$ . For 10 POD modes the reconstructions could recover 99 % of the energy for all variables and, therefore, should lead to an accurate flow representation. For this particular case, the SPOD technique is not required since the modes are almost monochromatic and, hence, do not require filtering.

One can observe from the figures that the computations of the flow using the ROM framework show good agreement compared to those obtained by the FOM. For the current Reynolds number, the flow develops a typical von Kármán vortex street along the cylinder wake. The periodical pattern of the vortex shedding can be observed in the contours of $z$ -vorticity. Noise generation also occurs in the current unsteady compressible flow simulation. In this case, pressure fluctuations along the cylinder surface are scattered to the far-field and can be observed in the contours of dilatation shown in figure 5. Both the near-field hydrodynamics and the far-field acoustics are recovered by the ROM. The reconstruction using 10 POD modes show excellent agreement with the FOM. When 2 POD modes are employed in the flow reconstruction, discrepancies between the ROM and FOM are evident from the figures. However, the main features of the flow are still recovered by the model. One should note that, despite the use of only two modes, the dynamical system is still stable beyond the training region.

In order to show a more qualitative evaluation of the model reconstructions, the density and $x$ -momentum fluctuation time histories are presented for the FOM and ROMs in figures  6 and 7, respectively. The figures on the left column show results for a probe located just behind the cylinder, close to the surface, at $(x,y)=(0.55,-0.06)$ . The cylinder has radius 0.5 and its centre is positioned in the origin of the Cartesian system. On the right column, results are obtained for a probe downstream of the cylinder wake, at $(x,y)=(1.1,-0.06)$ . Results are shown for both the training window period and beyond. When two POD modes are employed, the solutions show a less accurate representation of the dynamics observed in the FOM. The density reconstruction is that with the highest discrepancy and that is attributed to the lower energetic content achieved by the first two POD modes. One can notice that the ROM accurately reproduces the FOM results during and beyond the training window when 10 POD modes are employed in the reconstruction.

Figure 6.

Fluctuation time history of density.

Figure 7.

Fluctuation time history of $x$ -momentum.

In order to test the robustness of the method, we employ the DNN approach for the reconstruction of the transient regime of an incompressible cylinder flow. For this study, the two most energetic POD modes containing both the transient and limit cycle dynamics of the flow are obtained from Brunton et al. (Reference Brunton, Proctor and Kutz2016). These modes contain both the transient and limit cycle dynamics of the flow. Figure 8 shows the POD temporal modes used for training and testing the DNNs. The first 5000 temporal instants are used to train the models and one can see this data set represented by the black line to the left of the vertical line marking the end of the training window. To the right of the training window, the black line represents the test data which is the correct solution for the temporal dynamics. Figure 8 also shows the results obtained by the current DNN approach and by sparse regression. In the case of sparse regression, we employ the same model obtained by Brunton et al. (Reference Brunton, Proctor and Kutz2016) in table 11 of their supporting information. As one can observe, the DNN is able to recover accurately both the transient and limit cycle solutions of the test data for both POD modes. On the other hand, sparse regression reconstructs the transient portion of the dynamics but not the long-term prediction of the limit cycle. Several models obtained by the DNNs presented similar results compared to those shown in figure 8 and some had similar neural network architectures compared to that of table 5. These results show that the proposed DNN approach has good long-term predictive capabilities and can learn the transient features of the flow.

Figure 8.

Reconstruction of POD temporal modes for transient solution of incompressible flow past a cylinder.

4.2.2 Deep dynamic stall of plunging airfoil

The present study concerns a plunging SD7003 airfoil in deep dynamic stall. The flow conditions have a reduced frequency $k=\unicode[STIX]{x03C0}fL/U_{\infty }=0.5$ , where $f$ is the plunging frequency, $L$ is the chord of the airfoil and $U_{\infty }$ is the reference free stream velocity. The motion amplitude is set as $h_{o}/L=0.5$ with a static angle of attack $\unicode[STIX]{x1D6FC}_{0}=8^{\circ }$ . The chord Reynolds number based on the freestream velocity is $Re=60\,000$ and the freestream Mach number is $M=0.1$ . This flow condition is relevant for micro air vehicle applications and this case is selected based on the availability of published results from other high fidelity simulations (Visbal Reference Visbal2011) and particle image velocimetry data (Baik et al. Reference Baik, Rausch, Bernal and Ol2009; Ol et al. Reference Ol, Bernal, Kang and Shyy2009). In a previous work (D’Lelis et al. Reference D’Lelis, Wolf, Yeh and Taira2019), we performed a mesh refinement study and validation of the numerical solutions against the references above.

The present mesh configuration consists of a body-fitted O-grid with $441\times 300\times 64$ points in the streamwise, wall-normal and spanwise directions, respectively. The grid is generated with $70\,\%$ of the surface points located in the suction side of the airfoil to improve the capturing of finer flow scales developing along the turbulent region of the flow. Due to the favourable pressure gradients in the pressure side of the airfoil, the flow remains laminar along the entire cycle of the plunging motion. The trailing edge of the SD7003 airfoil is rounded in the current simulations with an arc of radius $r/L=0.0008$ . This procedure is required for retaining the smoothness of the metric terms computed by the high-order compact scheme. The spanwise domain is set as $z/c=0.4$ similarly to Visbal (Reference Visbal2011) and the dimensionless time step of the simulation is set as $\unicode[STIX]{x0394}t^{\ast }=((\unicode[STIX]{x0394}tU_{\infty })/L)=0.00008$ .

The plunge motion occurs with an effective angle of attack in the range of $-6^{\circ }\leqslant \unicode[STIX]{x1D6FC}\leqslant 22^{\circ }$ . Defining $\unicode[STIX]{x1D713}$ as the angular position in the plunging cycle, we say that at $\unicode[STIX]{x1D713}=0^{\circ }$ the airfoil has no velocity in the $y$ -direction and is at the top-most position of the plunging motion. At $\unicode[STIX]{x1D713}=90^{\circ }$ it has the highest velocity in the $y$ -direction downwards and, at $\unicode[STIX]{x1D713}=180^{\circ }$ , it has no velocity and is at the bottom-most position of the plunging motion. Finally, at $\unicode[STIX]{x1D713}=270^{\circ }$ it has the highest velocity in the $y$ -direction upwards. In summary, during the down-stroke instabilities begin to form in the suction side of the airfoil, growing and eventually breaking into finer structures. While this takes place, the dynamic stall vortex forms along the leading edge and is transported through the airfoil suction side increasing the overall lift and creating a nose-down pitching moment. As the leading-edge vortex (LEV) approximates the trailing edge, a trailing-edge vortex (TEV) forms pushing the LEV away from the airfoil. A more complete discussion of the flow dynamics can be found in Visbal (Reference Visbal2011) and D’Lelis et al. (Reference D’Lelis, Wolf, Yeh and Taira2019).

Figure 9 shows iso-surfaces of Q-criterion coloured by pressure and it is possible to observe the main flow features described above. In figures 9(a) and 9(b), it is possible to compare the 3-D solutions of the FOM and ROM, respectively, for the leading edge vortex formation. Figure 9(c,d) show a similar comparison for the instant where the trailing edge vortex forms. For both instants, one can observe that the ROM is able to reconstruct the larger scale features of the 3-D flow. A video comparing the FOM and ROM solutions for this flow configuration is provided together with the manuscript (movie 2). The ROM is trained using the first two cycles of the plunging motion and the solutions presented in figure 9 are computed for the fourth cycle, showing that the model is able to reproduce the 3-D flow dynamics beyond the training window. For this study the parameters employed in algorithm 2 are presented in table 6 and the best set of hyperparameters, obtained via random search, is listed in table 7. It is important to mention that, for this case, Bayesian optimization was unable to produce stable and accurate models.

Figure 9.

Iso-surfaces of Q criterion coloured by pressure at different instants of the plunge motion for the fourth cycle. The ROM was trained using only the first 2 cycles.

Table 6.

Model generation parameters for ROM of plunging airfoil in deep dynamic stall.

Table 7.

Best set of hyperparameters for ROM of plunging airfoil in deep dynamic stall (3-D flow).

The first 16 SPOD modes are employed to reduce the dimensionality of the input data in the 3-D flow reconstruction. Other POD reconstructions were tested with a different number of modes and it was observed that the first 16 modes contained the main features of the leading and trailing edge vortex formation. For example, adding 10 more modes did not improve significantly the solution and, beyond mode 26, the SPOD temporal modes presented complex behaviour, being composed of several frequencies and confusing the training stage of the DNNs. This issue could be improved by running the simulation for a longer period with a lower time step to improve convergence of the POD modes. One should note that the SPOD allows an energy shift of the modes to obtain frequency filtered temporal dynamics. In this sense, the high-frequency noise observed in the POD temporal modes is not discarded but is shifted to higher SPOD modes. This procedure allows a better pairing of POD modes if the coherent structures have periodicity. In the present study, due to the unsteady boundary conditions and the lack of symmetry in the flow, we do not expect pairing of the SPOD modes.

Figure 10.

Characteristics of POD and SPOD temporal modes for 3-D flow past airfoil under plunging motion.

For the present turbulent flow, the temporal modes of the standard snapshot POD are composed by several frequencies and, moreover, contain some noise that degrades the training of the DNNs. The SPOD provides the most energetic modes for specific frequency bands, allowing a better identification of the individual modes and smoothing out the temporal coefficients. Here, a box filter is applied to 25 % of the POD correlation matrix. We observed that filter values in the range of 20 % to 40 % provide good results for the model reconstructions. If a lower filter window is employed, noise can still be present in the temporal modes. On the other hand higher filter windows may considerably modify the temporal modes transforming them into quasi-sinusoidal signals. This is undesirable since the relevant dynamics of the flow can be lost and too many SPOD modes may be required for the construction of the ROM. An example of the impact of SPOD on the characteristics of temporal modes 2 and 16 can be seen in figure 10. In this example, the second POD mode exhibits high-frequency noise which is filtered by the SPOD. Meanwhile, POD mode 16 shows a complicated pattern due to the contribution of several frequencies. This temporal coefficient is considerably simplified by the application of SPOD as can be observed in the figure. It is clear that, different from the cylinder case, for the plunging airfoil the POD modes contain more complex dynamics which are composed of multiple Fourier modes.

In order to further evaluate the present DNN approach, the ROM reconstructions are also performed for spanwise-averaged solutions of the flow and the ROM is constructed following the procedures defined in § 2. We employ six plunging cycles (2244 snapshots) for computing the models, discarding initial flow transients. The training data contains the first two cycles (748 snapshots) of the FOM and the remaining data is used as the test set. The parameters employed in algorithm 2 are the same as those presented in table 6 and the best set of hyperparameters, obtained via random search, is listed in table 8. Here, we apply the box filter to 25 % of the POD correlation matrix and employ the first 10 SPOD modes for the ROM construction, which was sufficient to recover the most important features of the dynamic stall solution for the spanwise-averaged study.

Table 8.

Best set of hyperparameters for ROM of plunging airfoil in deep dynamic stall (spanwise-averaged flow).

Figures 11 and 12 present snapshots of density and $x$ -component of momentum, respectively, for the FOM and ROM. A video comparing these solutions for the spanwise-averaged flow is also provided together with the manuscript (movie 3). The flow features of the dynamic stall can be observed for different stages of the plunging cycle, given by different values of $\unicode[STIX]{x1D713}$ in the figures. It is possible to track the LEV over the suction side of the airfoil. The current ROM is able to recover the important dynamics of the leading edge vortex formation, its transport and ejection, besides the trailing edge vortex formation and ejection. One should note that the results presented in these figures are obtained for a plunging cycle beyond the training region. Therefore, the current ROM is capable of reproducing the flow dynamics for the test set. The high wavenumber features shown in the FOM are not present in the SPOD modes and, hence, they are not recovered by the ROM. If additional SPOD modes were employed in the model reconstruction, these features would appear. However, the overall cost of the simulations would increase considerably if accurate higher POD modes were required. For the present study, we observed that the fine scale flow dynamics are related to higher modes since most of the energetic content is related to the airfoil plunging motion represented by the first POD modes.

Figure 11.

Density contours at several phase angles for the fifth plunging cycle.

Figure 12.

Contours of $x$ -component of momentum at several phase angles for the fifth plunging cycle.

Fluctuation time histories are presented in figure 13 for density and $x$ -momentum. These properties are computed by the FOM and ROM at probe locations in the proximity of the leading and trailing edges, at $(x,y)=(0.07,0.07)$ and $(x,y)=(0.97,0.07)$ , respectively. It is important to mention that the airfoil leading edge is at the origin of the Cartesian system. The plots show the training and test regions demonstrating that the ROM is stable beyond the training set and can accurately reproduce the main dynamics of the flow. The fast oscillations resolved by the FOM are not represented due to the SPOD basis truncation. However, all the relevant features of the flow fluctuations are captured by the DNN model obtained with the first 10 SPOD modes.

Figure 13.

Fluctuation time histories computed by the FOM and ROM for probe locations in the proximity of the leading edge (a,c) and trailing edge (b,d).

In figure 14, one can observe a comparison between ROMs built using the current DNN approach and sparse regression. For this case, both the LASSO and Ridge regression techniques were tested with the original SINDy algorithm. Although the LASSO presented a higher computational cost, it provided better models compared to the Ridge regression. For this case, the computational cost of the DNN regression is, on average, 40 times higher than that of SINDy. However, as shown for SPOD modes 1 and 8, the SINDy algorithm was not able to provide stable models. The same can be said for the other SPOD modes. Sparse regression can learn the dynamics during the training window and also for a few cycles in the test set but its solution is not stable for long-time predictions. On the other hand, despite the higher cost, the DNN approach presents good long-time predictive capabilities with stable and accurate solutions beyond the training window.

Figure 14.

Reconstruction of SPOD temporal modes using the current DNN approach and the original SINDy algorithm with LASSO.

In order to evaluate the robustness of the present framework, we test the capability of the ROM in reconstructing the flow field using a time step different than that of the simulation. Figure 15 shows the reconstruction of SPOD temporal modes using different finite difference schemes and time steps. The derivatives of the temporal modes are computed via an explicit second-order scheme and by sixth- and tenth-order compact schemes. When the ROM is constructed using the same time step of the numerical simulation, the first 10 SPOD modes are highly resolved in time. For this case, all three schemes for computing the temporal derivatives provide accurate results and the random search is able to find a set of hyperparameters which results in an accurate ROM. However, when the DNN is trained using SPOD modes obtained by a time step five times larger than that of the simulation, the resolution of the higher POD temporal modes can become compromised. Then, it is important to use the compact schemes to obtain accurate representations of the derivatives. In this case, the random search could only find suitable ROMs using the compact schemes. As shown in figure 15, the models are still stable and reproduce the same dynamics observed with lower time steps. This shows that the current ROM framework is robust and can be constructed using a subset of the data of the FOM, with lower temporal resolution.

Figure 15.

Reconstruction of SPOD temporal modes using different finite difference schemes and time steps.

Finally, we compute the phase-averaged aerodynamic coefficients for the plunging SD7003 airfoil. Lift and drag coefficients are shown in figure 16 for the FOM and ROM using the six cycles computed in the training and test sets. Both the pressure and friction forces are accounted for in the calculation of the aerodynamic coefficients. However, for this case the pressure component dominates the aerodynamic loads. The effective incidence angles are depicted in the figures and the aerodynamic coefficients are calculated using the 10th-order compact scheme for the derivatives of the temporal modes and the first 10 SPOD modes. As a validation of the current solutions, the results from Visbal (Reference Visbal2011) are also shown. Results are presented for the ROMs computed using both the pressure and kinetic energy norms. The DNNs could find stable and accurate models for both POD norms and the aerodynamic coefficients obtained for the best models are compared in the figure. As one can see, the present ROM solutions show good comparisons to the current LES and Visbal (Reference Visbal2011). The total simulation cost of the FOM for this case was 100 000 h. On the other hand, the full cost of the ROM, including the optimization of the hyperparameters, plus the training and evaluation of 700 models, required approximately seven hours. Hence, the computational cost for the training procedure is around 100 models per hour in a single GPU. Once a ROM is chosen the cost of the simulation is reduced to a few seconds.

Figure 16.

Phase-averaged aerodynamic coefficients computed by Visbal (Reference Visbal2011), the current LES (FOM), and the DNN approach with POD norms based on pressure ( $P$ ) and kinetic energy ( $KE$ ).

5 Conclusions

We present a methodology for constructing ROMs combining flow modal decomposition and regression analysis via DNNs. The framework is implemented in a context similar to that of the SINDy algorithm recently proposed in the literature. The details of the methodology are described including algorithm charts which facilitate the understanding and implementation of the proposed framework. The source code can be downloaded from http://cces.unicamp.br/software/.

The method is tested for different problems involving nonlinear dynamical systems: the canonical nonlinear damped oscillator, the flow past a circular cylinder at a low Reynolds number and the turbulent flow past a plunging SD7003 airfoil under deep dynamic stall. For the previous two cases, the compressible Navier–Stokes equations are solved in full contravariant form and additional non-inertial terms are added to the equations to simulate the airfoil plunging motion. Numerical simulations are performed using a high-order compact finite difference flow solver. Then, the high fidelity results from the simulations are used as the input data for the construction of the ROMs.

To create stable and accurate ROMs of more complex flows, the application of the SPOD was found necessary. Hence, in order to filter high frequency content present in the temporal modes of the classical snapshot POD, we apply a filter function to the correlation matrix (SPOD approach). The proposed numerical framework allows the prediction of the flow field beyond the training window using larger time increments than those employed by the FOM, which demonstrates the robustness of the current ROMs constructed via DNNs. The resolution of the numerical schemes used for computation of the POD temporal mode derivatives is shown to have an important role when larger time increments are employed in the construction of the ROMs. In this case, the higher POD modes are composed by a broad range of frequencies and the accurate representation of the temporal derivatives is crucial for obtaining ROMs via regression analysis.

A discussion regarding the optimization of hyperparameters for obtaining the best ROMs via DNNs is provided together with a description of the effects of individual hyperparameters on model performance. Here, we test the random search and the Bayesian optimization which are the most widely used procedures for hyperparameter optimization. To report the performance of each model from a set of parameters, we compute the MAE over the training data and choose the candidate models with lower MAE values. Then, the best model is chosen based on the differences between ROM and FOM solutions over the validation set. Following this approach, the random search produced the best models at lower computational costs for all cases investigated in this work. It is worth mentioning that, in order to reduce the set of parameters for optimization, we first selected the hyperparameters related to the optimization step, such as the learning rate and the number of iterations, by manual search. We also chose the exponential linear unit activation function over the hyperbolic tangent function for the current problems since it provided results with fewer iterations.

Using 10 POD modes, 99 % of the flow energy is recovered in the cylinder flow study. The ROM obtained for this case shows an excellent agreement with the FOM. Even when two POD modes are employed in the reconstruction, the ROM is still stable beyond the training set and captures most of the dynamics of the vortex shedding and sound wave propagation. Both the DNN and the original SINDy approaches are tested for the reconstruction of the transient regime of an incompressible flow past a cylinder. In this case, regression is performed for the two most energetic POD modes of the flow, obtained from Brunton et al. (Reference Brunton, Proctor and Kutz2016). The DNN model is able to accurately recover both the transient and limit cycle solutions of the test data for both POD modes. On the other hand, sparse regression reconstructs the transient portion of the dynamics but not the long-term prediction of the limit cycle.

In the dynamic stall configuration, the flow is turbulent along the airfoil suction side, mostly during the downstroke motion. The complex flow dynamics of this case exhibit unpaired POD modes with high frequency noise. We show how the SPOD approach modifies the temporal modes for this study. Reduced-order models are constructed both for the full 3-D and the spanwise-averaged flow solutions. In both cases, the ROMs are able to capture the dynamics of the leading edge stall vortex, including its formation, transport and ejection, and of the trailing edge vortex. The total simulation cost of the FOM for this case was 100 000 h. On the other hand, the full cost of the ROM was approximately seven hours in a single GPU. This cost already includes the optimization of hyperparameters, plus the training and evaluation of 700 models. When the best model is chosen, the cost of the simulation is reduced to a few seconds. The data obtained by the DNN-ROMs were used to compute aerodynamic coefficients of the dynamic stall and good agreement was found compared to the LES used as the FOM. Again, we compare models obtained by DNNs and sparse regression and show that, despite the higher computational cost, the DNN approach presents good long-time predictive capabilities with stable and accurate solutions beyond the training window. On the other hand, the best model obtained from sparse regression could learn the dynamics during the training window and also for a few cycles in the test set but its solution was not stable for long-time predictions. We expect that, in a future work, the current methodology can be further improved for applications in flow control and extrapolation of flow configurations other than those used for training.

Acknowledgements

The authors acknowledge the financial support received from Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP, under grant nos. 2013/08293-7 and 2013/07375-0, and from Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, under grant nos. 304335/2018-5 and 407842/2018-7. The support provided by CNPq through a scholarship for the first author is likewise gratefully appreciated. We also thank CEPID-CeMEAI for providing the computational resources used in this work. We gratefully acknowledge B. D’Lelis and T. Ricciardi for providing the high fidelity dynamic stall and cylinder data employed in the current work.

Supplementary movies

Supplementary movies are available at https://doi.org/10.1017/jfm.2019.358.

References

Abadi, M. et al. 2016 TensorFlow: Large-scale machine learning on heterogeneous distributed systems. arXiv:1603.04467.Google Scholar
Akaike, H. 1973 Information Theory and an Extension of the Maximum Likelihood Principle, pp. 199–213. Springer.Google Scholar
Aris, R. 1989 Vectors, Tensors, and the Basic Equations of Fluid Mechanics. Dover.Google Scholar
Baik, Y. S., Rausch, J., Bernal, L. & Ol, M. V. 2009 Experimental investigation of pitching and plunging airfoils at Reynolds number between 1 × 104 and 6 × 104 . In AIAA Paper 2009-4030.Google Scholar
Balajewicz, M., Tezaur, I. & Dowell, E. 2016 Minimal subspace rotation on the Stiefel manifold for stabilization and enhancement of projection-based reduced order models for the compressible Navier–Stokes equations. J. Comput. Phys. 321, 224–241.Google Scholar
Beam, R. M. & Warming, R. F. 1978 An implicit factored scheme for the compressible Navier–Stokes equations. AIAA J. 16 (4), 393–402.Google Scholar
Bengio, Y. 2009 Learning deep architectures for AI. Found. Trends Mach. Learn. 2 (1), 1–127.Google Scholar
Bergstra, J. & Bengio, Y. 2012 Random search for hyper-parameter optimization. J. Mach. Learn. Res. 13, 281–305.Google Scholar
Brochu, E., Cora, V. M. & Freitas, N.2010 A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv:1012.2599.Google Scholar
Brunton, S. L. & Noack, B. R. 2015 Closed-loop turbulence control: progress and challenges. Appl. Mech. Rev. 67, 050801.Google Scholar
Brunton, S. L., Proctor, J. L. & Kutz, N. J. 2016 Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl Acad. Sci. USA 113 (15), 3932–3937.Google Scholar
Carlberg, K., Bou-Mosleh, C. & Farhat, C. 2011 Efficient non-linear model reduction via a least-squares Petrov–Galerkin projection and compressive tensor approximations. Intl J. Numer. Meth. Engng 86 (2), 155–181.Google Scholar
Cazemier, W., Verstappen, R. W. C. P. & Veldman, A. E. P. 1998 Proper orthogonal decomposition and low-dimensional models for driven cavity flows. Phys. Fluids 10 (7), 1685–1699.Google Scholar
Chaturantabut, S. & Sorensen, D. 2010 Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput. 32, 2737–2764.Google Scholar
Clainche, S. & Vega, J. M. 2017 Higher order dynamic mode decomposition to identify and extrapolate flow patterns. Phys. Fluids 29, 084102.Google Scholar
Clevert, D., Unterthiner, T. & Hochreiter, S.2015 Fast and accurate deep network learning by exponential linear units (ELUs). arXiv:1511.07289.Google Scholar
D’Lelis, B. O. R., Wolf, W. R., Yeh, C. & Taira, K.2019 High-fidelity simulation and flow control of a plunging airfoil under deep dynamic stall. AIAA Paper 2019-1395.Google Scholar
Goodfellow, I., Bengio, Y. & Courville, A. 2016 Deep Learning. MIT Press.Google Scholar
Green, M. A., Rowley, C. W. & Haller, G. 2007 Detecting vortex formation and shedding in cylinder wakes using Lagrangian coherent structures. J. Fluid Mech. 572, 111–120.Google Scholar
Haller, G. 2015 Lagrangian coherent structures. Annu. Rev. Fluid Mech. 47, 137–162.Google Scholar
Juang, J. N. & Pappa, R. S. 1985 An eigensystem realization algorithm for modal parameter identification and model reduction. J. Guid. Control Dyn. 8, 620–627.Google Scholar
Kennedy, C. A., Carpenter, M. H. & Lewis, R. M.1999 Low-storage, explicit Runge–Kutta schemes for the compressible Navier–Stokes equations. ICASE Report No. 99-22.Google Scholar
Kingma, D. P. & Ba, J. 2014 Adam: a method for stochastic optimization. In Proceedings of the 3rd International Conference on Learning Representations (ICLR). arXiv:1412.6980.Google Scholar
Krizhevsky, A.2009 Learning multiple layers of features from tiny images. Tech. Rep.Google Scholar
Kutz, J. 2017 Deep learning in fluid dynamics. J. Fluid Mech. 814, 1–4.Google Scholar
Lele, S. K. 1992 Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103 (1), 16–42.Google Scholar
Ling, J., Kurzawski, A. & Templeton, J. 2016 Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. J. Fluid Mech. 807, 155–166.Google Scholar
Ling, J. & Templeton, J. 2015 Evaluation of machine learning algorithms for prediction of regions of high Reynolds averaged Navier Stokes uncertainty. Phys. Fluids 27, 085103.Google Scholar
Luhar, M., Sharma, A. S. & McKeon, B. J. 2014 On the structure and origin of pressure fluctuations in wall turbulence: predictions based on the resolvent analysis. J. Fluid Mech. 751, 38–70.Google Scholar
Lumley, J. L. 1967 The structure of inhomogeneous turbulence. In Atmospheric Turbulence and Wave Propagation, pp. 166–178. Nauka.Google Scholar
Nagarajan, S., Lele, S. K. & Ferziger, J. H. 2003 A robust high-order compact method for large eddy simulation. J. Comput. Phys. 191 (2), 392–419.Google Scholar
Nair, V. & Hinton, G. E. 2010 Rectified linear units improve restricted Boltzmann machines. In Proceedings of the 27th International Conference on International Conference on Machine Learning, pp. 807–814. Omnipress.Google Scholar
Navarro, L. C., Navarro, A. K. W., Rocha, A. & Dahab, R. 2018 Connecting the dots: Toward accountable machine-learning printer attribution methods. J. Visual Commun. Image Represent. 53, 257–272.Google Scholar
Noack, B. R., Afanasiev, K., Morzynski, M., Tadmor, G. & Thiele, F. 2003 A hierarchy of low-dimensional models for the transient and post-transient cylinder wake. J. Fluid Mech. 497, 335–363.Google Scholar
Ol, M. V., Bernal, L., Kang, C. & Shyy, W. 2009 Shallow and deep dynamic stall for flapping low Reynolds number airfoils. Exp. Fluids 46 (5), 883–901.Google Scholar
Olson, B. J.2012 Large-eddy simulation of multi-material mixing and over-expanded nozzle flow. PhD thesis, Stanford University.Google Scholar
Östh, J., Noack, B. R., Krajnovi, S., Barros, D. & Borée, J. 2014 On the need for a nonlinear subscale turbulence term in POD models as exemplified for a high-Reynolds-number flow over an Ahmed body. J. Fluid Mech. 747, 518–544.Google Scholar
Pan, S. & Duraisamy, K.2018 Long-time predictive modeling of nonlinear dynamical systems using neural networks. arXiv:1805.12547.Google Scholar
Protas, B., Noack, B. R. & Östh, J. 2015 Optimal nonlinear eddy viscosity in Galerkin models of turbulent flows. J. Fluid Mech. 766, 337–367.Google Scholar
Ribeiro, J. H. M. & Wolf, W. R. 2017 Identification of coherent structures in the flow past a NACA0012 airfoil via proper orthogonal decomposition. Phys. Fluids 29 (8), 085104.Google Scholar
Rockwood, M. P., Taira, K. & Green, M. A. 2016 Detecting vortex formation and shedding in cylinder wakes using Lagrangian coherent structures. AIAA J. 55, 15–23.Google Scholar
Rowley, C. W., Colonius, T. & Murray, R. M. 2004 Model reduction for compressible flows using POD and Galerkin projection. Physica D: Nonlinear Phenomena 189 (12), 115–129.Google Scholar
Rudy, S. H., Brunton, S. L., Proctor, J. L. & Kutz, N. J. 2017 Data-driven discovery of partial differential equations. Sci. Adv. 3 (4), e1602614.Google Scholar
Rudy, S. H., Kutz, J. N. & Brunton, S. L.2018 Deep learning of dynamics and signal-noise decomposition with time-stepping constraints. arXiv:1808.02578.Google Scholar
San, O. & Maulik, R. 2018 Extreme learning machine for reduced order modeling of turbulent geophysical flows. Phys. Rev. E 97, 042322.Google Scholar
Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28.Google Scholar
Schwarz, G. 1978 Estimating the dimension of a model. Ann. Stat. 6, 461–464.Google Scholar
Sharma, A. S., Moarref, R., McKeon, B. J., Park, J. S., Graham, M. D. & Willis, A. P. 2016 Low-dimensional representations of exact coherent states of the Navier–Stokes equations from the resolvent model of wall turbulence. Phys. Rev. E 93, 021102.Google Scholar
Sieber, M., Paschereit, C. O. & Oberleithner, K. 2016 Spectral proper orthogonal decomposition. J. Fluid Mech. 792, 798–828.Google Scholar
Sirovich, L. 1986 Turbulence and the dynamics of coherent structures. Part I. Coherent structures. Q. Appl. Maths 45, 561–571.Google Scholar
Snoek, J., Larochelle, H. & Adams, R. P. 2012 Practical Bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems 25, pp. 2951–2959. Curran Associates.Google Scholar
Ströfer, C. M., Wu, J. L., Xiao, H. & Paterson, E. 2019 Data-driven, physics based feature extraction from fluid flow fields using convolutional neural networks. Commun. Comput. Phys. 25, 625–650.Google Scholar
Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W., Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S., Theofilis, V. & Ukeiley, L. S. 2017 Modal analysis of fluid flows: an overview. AIAA J. 55 (12), 4013–4041.Google Scholar
Towne, A., Schmidt, O. T. & Colonius, T. 2018 Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. J. Fluid Mech. 847, 821–867.Google Scholar
Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L. & Kutz, N. J. 2014 On dynamic mode decomposition: theory and applications. J. Comput. Dyn. 1, 391–421.Google Scholar
Visbal, M. R. 2011 Numerical investigation of deep dynamic stall of a plunging airfoil. AIAA J. 49, 2152–2170.Google Scholar
Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P. & Koumoutsakos, P. 2018 Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks. Proc. R. Soc. Lond. A 474 (2213).Google Scholar
Wan, Z. Y., Vlachas, P., Koumoutsakos, P. & Sapsis, T. 2018 Data-assisted reduced-order modeling of extreme events in complex dynamical systems. PLOS ONE 13 (5), 1–22.Google Scholar
Wang, J. X., Wu, J. L. & Xiao, H. 2017 Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data. Phys. Rev. Fluids 2, 034603.Google Scholar
Wolf, W. R.2011 Airfoil aeroacoustics: LES and acoustic analogy predictions. PhD thesis, Stanford University.Google Scholar
Wolf, W. R., Azevedo, J. L. F. & Lele, S. K. 2012a Convective effects and the role of quadrupole sources for aerofoil aeroacoustics. J. Fluid Mech. 708, 502–538.Google Scholar
Wolf, W. R., Lele, S. K., Jothiprasard, G. & Cheung, L.2012 Investigation of noise generated by a DU96 airfoil. AIAA Paper 2012-2055.Google Scholar
Xavier, G. & Yoshua, B. 2010 Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics (ed. Teh, Y. W. & Titterington, M.).Google Scholar
Zimmermann, R. & Willcox, K. 2016 An accelerated greedy missing point estimation procedure. J. Sci. Comput. 38, 2827–2850.Google Scholar
Figure 0

Figure 1.Schematic of the current method used for construction of ROMs via DNNs and flow modal decomposition.

Figure 1

Figure 2.An example of a deep feedforward network with $N$ inputs, several hidden layers (HL), and one output layer with $N$ outputs.

Figure 2

Table 1.Effects of hyperparameters on model performance.

Figure 3

Figure 3.Comparison between true model (solid black line) and ROM prediction (dashed cyan line) for nonlinear oscillator.

Figure 4

Table 2.Model generation parameters for ROM of nonlinear oscillator.

Figure 5

Table 3.Best set of hyperparameters for ROM of nonlinear oscillator.

Figure 6

Table 4.Model generation parameters for ROM of compressible flow past a cylinder.

Figure 7

Table 5.Best set of hyperparameters for ROM of compressible flow past a cylinder.

Figure 8

Figure 4.Contours of $\boldsymbol{z}$-vorticity, $t=410$.

Figure 9

Figure 5.Contours of divergence of velocity, $t=410$.

Figure 10

Figure 6.Fluctuation time history of density.

Figure 11

Figure 7.Fluctuation time history of $x$-momentum.

Figure 12

Figure 8.Reconstruction of POD temporal modes for transient solution of incompressible flow past a cylinder.

Figure 13

Figure 9.Iso-surfaces of Q criterion coloured by pressure at different instants of the plunge motion for the fourth cycle. The ROM was trained using only the first 2 cycles.

Figure 14

Table 6.Model generation parameters for ROM of plunging airfoil in deep dynamic stall.

Figure 15

Table 7.Best set of hyperparameters for ROM of plunging airfoil in deep dynamic stall (3-D flow).

Figure 16

Figure 10.Characteristics of POD and SPOD temporal modes for 3-D flow past airfoil under plunging motion.

Figure 17

Table 8.Best set of hyperparameters for ROM of plunging airfoil in deep dynamic stall (spanwise-averaged flow).

Figure 18

Figure 11.Density contours at several phase angles for the fifth plunging cycle.

Figure 19

Figure 12.Contours of $x$-component of momentum at several phase angles for the fifth plunging cycle.

Figure 20

Figure 13.Fluctuation time histories computed by the FOM and ROM for probe locations in the proximity of the leading edge (a,c) and trailing edge (b,d).

Figure 21

Figure 14.Reconstruction of SPOD temporal modes using the current DNN approach and the original SINDy algorithm with LASSO.

Figure 22

Figure 15.Reconstruction of SPOD temporal modes using different finite difference schemes and time steps.

Figure 23

Figure 16.Phase-averaged aerodynamic coefficients computed by Visbal (2011), the current LES (FOM), and the DNN approach with POD norms based on pressure ($P$) and kinetic energy ($KE$).

Lui and Wolf supplementary movie 1

Movie showing x-momentum for full and reduced order models (cylinder flow)

Download Lui and Wolf supplementary movie 1(Video)
Video 9.3 MB

Lui and Wolf supplementary movie 2

Movie showing Q-criterion colored by pressure for full and reduced order models (3D flow of plunging airfoil)

Download Lui and Wolf supplementary movie 2(Video)
Video 49.4 MB

Lui and Wolf supplementary movie 3

Movie showing x-momentum for full and reduced order models (spanwise-averaged flow of plunging airfoil)

Download Lui and Wolf supplementary movie 3(Video)
Video 12.1 MB