---
abstract: 'An approximation model based on convolutional neural networks (CNNs) is proposed for flow field predictions. The CNN is used to predict the velocity and pressure field in unseen flow conditions and geometries given the pixelated shape of the object. In particular, we consider Reynolds Averaged Navier-Stokes (RANS) flow solutions over airfoil shapes. The CNN can automatically detect essential features with minimal human supervision and shown to effectively estimate the velocity and pressure field orders of magnitude faster than the RANS solver, making it possible to study the impact of the airfoil shape and operating conditions on the aerodynamic forces and the flow field in near-real time. The use of specific convolution operations, parameter sharing, and robustness to noise are shown to enhance the predictive capabilities of CNN. We explore the network architecture and its effectiveness in predicting the flow field for different airfoil shapes, angles of attack, and Reynolds numbers.'
author:
- Saakaar Bhatnagar
- Yaser Afshar
- Shaowu Pan
- Karthik Duraisamy
- Shailendra Kaushik
date: 'Received: date / Revised version: date'
title: Prediction of Aerodynamic Flow Fields Using Convolutional Neural Networks
---

Introduction {#intro}
============

With advances in computing power and computational algorithms, simulation-based design and optimization has matured to a level that it plays a significant role in an industrial setting. In many practical engineering applications, however, the analysis of the flow field tends to be the most computationally intensive and time-consuming part of the process. These drawbacks make the design process tedious, time consuming, and costly, requiring a significant amount of user intervention in design explorations, thus proving to be a barrier between designers from the engineering process.

Data-driven methods have the potential to augment [@Duraisamy:2019] or replace [@Guo:2016] these expensive high-fidelity analyses with less expensive approximations. Learning representations from the data, especially in the presence of spatial and temporal dependencies, have traditionally been limited to hand-crafting of features by domain experts. Over the past few years, deep learning approaches [@Bengio:2009; @Schmidhuber:2015] have shown significant successes in learning from data, and have been successfully used in the development of novel computational approaches [@Raissi:2018; @Raissi:2018b; @Raissi:2019].

Deep learning presents a fast alternative solution as an efficient function approximation technique in high-dimensional spaces. Deep learning architectures such as deep neural networks (DNNs), routinely used in data mining, are well-suited for application on big, high-dimensional data sets, to extract multi-scale features.

Deep convolutional neural networks (CNN) belong to a class of DNNs, most commonly applied to the analysis of visual imagery. Previous works [@Lecun:1998; @Taylor:2010; @Zuo:2015] have illustrated the promise of CNNs to learn high-level features even when the data has strong spatial and temporal correlations. Increasing attention being received by CNNs in fluid mechanics partly originates from their potential benefit of flexibility in the shape representation and scalability for 3D and transient problems. Figure \[fig:fig1\] illustrates the simplified layout of a typical CNN, LeNet-5 [@Lecun:1998] applied to the handwritten digit recognition task.

The main advantage of a CNN is that it exploits the low dimensional high-level abstraction by convolution. The key idea of CNN is to learn the representation and then to use a fully connected standard layer to fit the relationship between the high-level representation and output.

State of the art in application of CNNs in fluid dynamics {#State_of_the_art}
---------------------------------------------------------

The use of deep neural networks in computational fluid dynamics recently has been explored in some rudimentary contexts.

[@Guo:2016] reported the analysis and prediction of non-uniform steady laminar flow fields around bluff body objects by employing a convolutional neural network (CNN). The authors reported a computational cost lower than that required for numerical simulations by GPU-accelerated CFD solver. Though this work was pioneering in the sense that it demonstrated generalization capabilities, and that CNNs can enable a rapid estimation of the flow field, emphasis was on qualitative estimates of the velocity field, rather than on precise aerodynamic characteristics.

[@Miyanawala:2017] used a CNN to predict aerodynamic force coefficients of bluff bodies at a low Reynolds number for different bluff body shapes. They presented a data-driven method using CNN and the stochastic gradient-descent for the model reduction of the Navier-Stokes equations in unsteady flow problems.

[@Lee:2017; @Lee:2018] used a generative adversarial network (GAN) to predict unsteady laminar vortex shedding over a circular cylinder. They presented the capability of successfully learning and predicting both spatial and temporal characteristics of the laminar vortex shedding phenomenon.

[@hennigh:2017b] presented an approach to use a DNN to compress both the computation time and memory usage of the Lattice Boltzmann flow simulations. The author employed convolutional autoencoders and residual connections in an entirely differentiable scheme to shorten the state size of simulation and learn the dynamics of this compressed form.

[@Tompson:2016] proposed a data-driven approach for calculating numerical solutions to the inviscid Euler equations for fluid flow. In this approach, an approximate inference of the sparse linear system used to enforce the Navier-Stokes incompressibility condition, the “pressure projection” step. This approach cannot guarantee an exact solution pressure projection step, but they showed that it empirically produces very stable divergence-free velocity fields whose runtime and accuracy is better than the Jacobi method while being orders of magnitude faster.

[@Zhang:2017] employed a CNN as feature extractor for a low dimensional surrogate modeling. They presented the potential of learning and predicting lift coefficients using the geometric information of airfoil and operating parameters like Reynolds number, Mach number, and angle of attack. However, the output is not the flow field around the airfoil but the pressure coefficients at several locations. It is unclear whether this model would have good performance in predicting the drag and pressure coefficient when producing the flow field at the same time.

The primary contribution of the present work is a framework that can be used to predict the flow field around different geometries under variable flow conditions. Towards this goal and following [@Guo:2016], we propose a framework with a general and flexible approximation model for near real-time prediction of non-uniform steady RANS flow in a domain based on convolutional neural networks. In this framework, the flow field can be extracted from simulation data by learning the relationship between an input feature extracted from geometry and the ground truth from a RANS simulation. Then without standard convergence requirements of the RANS solver, and its number of iterations and runtime, which are irrelevant to the prediction process, we can directly predict the flow behavior in a fraction of the time. In contrast to previous studies, the present work is focused on a more rigorous characterization of aerodynamic characteristics. The present study also improves on computational aspects. For instance, [@Guo:2016] use an separated decoder, whereas the present work employs shared-encoding and decoding layers, which are computationally efficient compared to the separated alternatives.

 Methodology {#Methodology}
============

CFD Simulation {#CFD}
--------------

In this work, flow computations and analyses are performed using the OVERTURNS CFD code [@Duraisamy:2005; @Lakshminarayan:2010]. This code solves the compressible RANS equations using a preconditioned dual-time scheme [@Pandya:2003]. Iterative solutions are pursued using the implicit approximate factorization method [@Pulliam:1981]. Low Mach preconditioning [@Turkel:1999] is used to improve both convergence properties and the accuracy of the spatial discretization. A third order Monotonic Upwind Scheme for Conservation Laws (MUSCL) [@VanLeer:1979] with Koren’s limiter [@Koren:1993] and Roe’s flux difference splitting [@Roe:1986] is used to compute the inviscid terms. Second order accurate central differencing is used for the viscous terms. The RANS closure is the SA [@Spalart:1992] turbulence model and $\gamma - \overline{Re_{\theta t}}$ model [@Medida:2011] is used to capture the effect of the flow transition. No-slip boundary conditions imposed on the airfoil surface. The governing equations are provided in the Appendix.

Simulations are performed over the S805 [@Somers:1997a], S809 [@Somers:1997b], and S814 [@Somers:2004] airfoils. S809 and S814 are among a family of airfoils which contain a region of pressure recovery along the upper surface which induces a smooth transition from laminar to turbulent flow (so-called “transition-ramp”). These airfoils are utilized in wind turbines [@Aranake:2012]. Computations are performed using structured C-meshes with dimensions $394 \times 124$ in the wrap-around and normal directions respectively. Figure \[fig:fig2\] shows the airfoils and their near-body meshes.

[0.45]{}

[0.45]{}

[0.45]{}

Simulations are performed at Reynolds numbers $0.5,~1,~2,~\text{and}~3 \times 10^6$, respectively, and a low Mach number of $0.2$ is selected to be representative of wind turbine conditions. At each Reynolds number, the simulation is performed for different airfoils with a sweep of angles of attack from $\alpha=0^{\circ}$ to $\alpha=20^{\circ}$. The OVERTURNS CFD code has been validated for relevant wind turbine applications in [@Aranake:2012].

Convolutional Neural Networks {#CNN}
-----------------------------

In this study, we consider the convolutional neural network to extract relevant features from fluid dynamics data and to predict the entire flow field in near real-time. The objective is a properly trained CNN which can construct the flow field around an airfoil in a non-uniform turbulence field, using only the shape of the airfoil and fluid flow characteristics of the free stream in the form of the angle of attack and Reynolds number. In this section, we describe the structure and components of the proposed CNN.

Network Structure {#Structure}
-----------------

To develop suitable CNN architectures for variable flow conditions and airfoil shapes, we build our model based on an encoder-decoder CNN, similar to the model proposed by [@Guo:2016]. Encoder-decoder CNNs are most widely used for machine translation from a source language to a target language [@Chollampatt:2018]. The encoder-decoder CNN has three main components: a stack of convolution layers, followed by a dense layer and subsequently another stack of convolution layers. Figure \[fig:fig3\] illustrates the proposed CNN architecture designed in this work.

[@Guo:2016] used a shared-encoder but separated decoder. We conjecture that the separated decoder may be a limiting performance factor. To address this issue, we designed shared-encoding and decoding layers in our configuration, which save computations compared to the separated alternatives. Explicitly, the weights of the layers of the decoder are shared where they are responsible for extracting high-level representations of pressure and different velocity components. This design provides the same accuracy of the separated decoders but, it is almost utilized fifty percent fewer parameters compared to the separated alternatives. Also, in the work of [@Guo:2016], the authors used only one low Reynolds number for all the experiments, but here, the architecture is trained with four high Reynolds numbers, three airfoils with different shapes and 21 different angles of attacks. In this architecture, we use three convolution layers both in the shared-encoding and decoding parts.

The inputs to the network are the airfoil shape and the free stream conditions of the fluid flow. We use the convolution layers to extract the geometry representation from the inputs. The decoding layers use this representation in convolution layers and generate the mapping from the extracted geometry representation to the pressure field and different components of the velocity. The network uses the Reynolds number, the angle of attack, and the shape of the airfoil in the form of $150 \times 150$ 2D array created for each data entry. The geometry representation has to be extracted from the RANS mesh and fed to the network with images. Using images in CNNs allows encoding specific properties into the architecture, and reducing the number of parameters in the network.

 Geometry Representation {#Geometry}
------------------------

A wide range of approaches are employed to capture shape details and to classify points into a learnable format. Among popular examples are methods like implicit functions in image reconstruction [@Hoppe:1992; @Carr:2001; @Kazhdan:2006; @Fuhrmann:2014], or shape representation and classification [@Zhang:2004a; @Ling:2007; @Xu:2015; @Fernando:2015]. In applications such as rendering and segmentation and in extracting structural information of different shapes, signed distance functions (SDF) are widely used. SDF provides a universal representation of different geometry shapes and represents a grid sampling of the minimum distance to the surface of an object. It also works efficiently with neural networks for shape learning. In this study, to capture shape details in different object representations, and following [@Guo:2016; @Prantl:2017], we use the SDF sampled on a Cartesian grid. [@Guo:2016] reported the effectiveness of SDF in representing the geometry shapes for CNNs. The authors empirically showed that the values of SDF on the Cartesian grid provide not only local geometry details but also contain additional information on the global geometry structure.

Signed Distance Function {#SDF}
------------------------

A mathematical definition of the signed distance function of a set of points $\textbf{X}$ determines the minimum distance of each given point $\textbf{x} \in \textbf{X}$ from the boundary of an object $\partial\Omega$. $$\begin{aligned}
\textrm{SDF}(\textbf{x}) = \left\{\begin{matrix}
d(\textbf{x}, \partial \Omega) & ~\textbf{x}\notin\Omega \\
0 & ~~~\textbf{x}\in\partial\Omega\\
-d(\textbf{x}, \partial \Omega) & ~\textbf{x}\in\Omega 
\end{matrix}\right.,
\label{eq:eq1}\end{aligned}$$ where, $\Omega$ denotes the object, and $d(\textbf{x}, \partial \Omega) = \min_{\textbf{x}_I \in \partial \Omega}{\left( |\textbf{x} - \textbf{x}_I |\right )}$ measures the shortest distance of each given point $\textbf{x}$ from the object boundary points. The distance sign determines whether the given point is inside or outside of the object. Figure \[fig:fig4\] illustrates the signed distance function contour plot for a S814 airfoil.

Here, the SDF has positive values at points which are outside of the airfoil, and it decreases as the point approaches the boundary of the airfoil where the SDF is zero, and it takes negative values inside the airfoil. Fast marching method [@Sethian:1996] and fast sweeping method [@Zhao:2005] are among the popular algorithms for calculating the signed distance function. To generate a signed distance function, we use the CFD input structured C-mesh information and define the points around the object (airfoil). Figure \[fig:fig5\] shows the C-mesh representation of an airfoil (S814) and its boundary points on a Cartesian grid.

We find the distance of Cartesian grid points from the object boundary points, using the fast marching method [@Sethian:1996]. To find out whether a given point is inside, outside, or just on the surface of the object, we search the boundary points and compute the scalar product between the normal vector at the nearest boundary points and the vector from the given point to the nearest one and judge the function sign from the scalar product value. For other non-convex objects, one can also use different approaches of crossing number or winding number method which are common in ray casting [@Foley:1995].

After pre-processing the CFD mesh files, we use the SDF as an input to feed the encoder-decoder architecture with multiple layers of convolutions. Convolution layers in the encoding-decoding part extract all the geometry features from the SDF.

Convolutional Encoder-Decoder Approach {#Encoder}
--------------------------------------

To learn all the geometry features from an input SDF, we compose the encoder and decoder with convolution layers and convolutional filters. Every convolutional layer is composed of 300 convolutional filters. Therefore, a convolution produces a set of 300 activation maps. Every convolution in our design is wrapped by a non-linear Swish activation function [@Ramachandran:2017]. Swish is defined as $x.\sigma(\beta x)$ where $\sigma(z)=(1+exp(-z))^{-1}$ is the sigmoid function and $\beta$ is either a constant or a trainable parameter. The resulting activation maps are the encoding of the input in a low dimensional space of parameters to learn. The decoding operation is a convolution as well, where the encoding architecture fixes the hyper-parameters of the decoding convolution. Compared to the encoding convolution layer, here a convolution layer has reversed forward and backward passes. This inverse operation is sometimes referred to “deconvolution”. The decoding operation unravels the high-level features encoded and transformed by the encoding layers and generates the mapping to the pressure field and different components of the velocity. When we use the CNN, neurons in the same feature map plane have identical weights so that the network can study concurrently, and it learns implicitly from the training data. The training phase of the CNN comprises the input function, the feed-forward process, and the back-propagation process.

Data Preparation {#Data}
----------------

In total, a set of 252 RANS simulations were performed. This data includes our CFD predictions for three different S805, S809, and S814 airfoils. The training data-set consists of $85$ percent of the full set, and the remaining data sets are used for testing, as shown in Fig. \[fig:fig6\].

The test points are chosen uniformly at random on the feature space, providing an unbiased evaluation of a model fit on the training data-set while tuning the model’s hyper-parameters.

Figure \[fig:fig7\] shows the x-component of the velocity field ($U$) around the S814 airfoil on the structured C-mesh. The simulation is performed at an angle of attack of $\alpha = 9^\circ$ and with the Reynolds number of $3\times10^6$.

The CFD data has to be interpolated onto a $150\times 150$ Cartesian grid which contains the SDF. A triangulation-based scattered data interpolation method [@Amidror:2002] is used. After the interpolation of the data to the Cartesian grid, the interior points masked and the velocity is set to zero. The comparison of the reconstructed data in Fig. \[fig:fig8\] and the CFD data in Fig. \[fig:fig7\] shows evidence of interpolation errors.

The interpolated data is normalized using the standard score normalization by subtracting the mean from the data and dividing the difference by the standard deviation of the data. Scaling the data causes each feature to contribute approximately proportionately to the training, and also results in a faster convergence of the network [@Aksoy:2000].

Network Training and Hyper-parameter Study {#Training}
------------------------------------------

The network learns different weights during the training phase to predict the flow fields. In each iteration, a batch of data undergoes the feed-forward process followed by a back-propagation (see Sec. \[Encoder\]). For a given set of input and ground truth data, the model minimizes a total loss function which is a combination of two specific loss functions and an L2 regularization as follows: $$\begin{aligned}
\label{eq:MSE}
&\text{MSE}_\text{shared} = \frac{1}{m (n_x-2) (n_y-2)}~\sum_{l=1}^{m} \sum_{j=2}^{n_y-1} \sum_{i=2}^{n_x-1}\\
&\nonumber \quad \big[(U^l_{{ij}_\text{truth}}-U^l_{{ij}_\text{pred}})^2~+~(V^l_{{ij}_\text{truth}}-V^l_{{ij}_\text{pred}})^2~+~\\
&\nonumber \quad ~(P^l_{{ij}_\text{truth}}-P^l_{{ij}_\text{pred}})^2 \big], \\
\label{eq:GSshared}
&\text{GS}_\text{shared} = \frac{1}{6m(n_x-2)(n_y-2)}~\sum_{l=1}^{m} \sum_{j=2}^{n_y-1} \sum_{i=2}^{n_x-1} \\
&\nonumber \quad [(\frac{\partial P^l}{\partial x}_{{ij}_\text{truth}} - \frac{\partial P^l}{\partial x}_{{ij}_\text{pred}})^2 ~+~(\frac{\partial P^l}{\partial y}_{{ij}_\text{truth}} - \frac{\partial P^l}{\partial y}_{{ij}_\text{pred}})^2 ~+~ \\
&\nonumber \quad ~ (\frac{\partial U^l}{\partial x}_{{ij}_\text{truth}} - 
\frac{\partial U^l}{\partial x}_{{ij}_\text{pred}})^2 ~+~(\frac{\partial U^l}{\partial y}_{{ij}_\text{truth}} -\frac{\partial U^l}{\partial y}_{{ij}_\text{pred}})^2 ~+~ \\
&\nonumber \quad ~ (\frac{\partial V^l}{\partial x}_{{ij}_\text{truth}} - \frac{\partial V^l}{\partial x}_{{ij}_\text{pred}})^2 ~+~(\frac{\partial V^l}{\partial y}_{{ij}_\text{truth}} -
\frac{\partial V^l}{\partial y}_{{ij}_\text{pred}})^2],\\
\label{eq:L2}
&\text{L2}_\text{regularization} = \frac{1}{2m}\sum_{l=1}^{L}\sum_{i=1}^{n_l}(\theta^l_{i})^2,\end{aligned}$$ where $U$, and $V$ are the x-component and y-component of the velocity field respectively, and $P$ is the scalar pressure field. $m$ is the batch size, $n_x$ is the number of grid points along the x-direction, $n_y$ is the number of grid points along the y-direction, and $L$ is the number of layers with trainable weights, and $n_l$ represents number of trainable weights in layer $l$. MSE is the mean squared error, and GS is gradient sharpening or gradient difference loss (GDL) [@Mathieu:2015; @Lee:2018]. In this paper, we use gradient sharpening based on a central difference operator. The network was trained for $30,000$ epochs with a batch size of $214$ data points, which took 33 GPU hours. For the separated decoder, the following loss functions are used: $$\begin{aligned}
\label{eq:MSEseparate}
&\text{MSE}_\text{separated} = \frac{1}{m (n_x-2) (n_y-2)}~\sum_{l=1}^{m} \sum_{j=2}^{n_y-1} \sum_{i=2}^{n_x-1}\\
&\nonumber \quad \big[(X^l_{{ij}_\text{truth}}-X^l_{{ij}_\text{pred}})^2],\\
\label{eq:GSseparate}
&\text{GS}_\text{separated} = \frac{1}{2m(n_x-2)(n_y-2)}~\sum_{l=1}^{m} \sum_{j=2}^{n_y-1} \sum_{i=2}^{n_x-1} \\
&\nonumber \quad [(\frac{\partial X^l}{\partial x}_{{ij}_\text{truth}} - \frac{\partial X^l}{\partial x}_{{ij}_\text{pred}})^2 ~+~(\frac{\partial X^l}{\partial y}_{{ij}_\text{truth}} - \frac{\partial X^l}{\partial y}_{{ij}_\text{pred}})^2],\end{aligned}$$ where $X$ stands for $U,~V$ or $P$.

Finding the optimal set of hyper-parameters for the network is an empirical task and is done by performing a grid search consisting of an interval of values of each hyper-parameter, and training many networks with several different combinations of these hyper-parameters. The resulting networks are compared based on generalization tendency and the difference between the truth and prediction.

Results and Discussion {#Results}
======================

We first show the capability of the designed network architecture to accurately estimate the velocity and pressure field around different airfoils given only the airfoil shape. Then, we quantitatively assess the error measurement followed by a sequence of results which demonstrate usability, accuracy and effectiveness of the network.

Figure \[fig:fig9\] illustrates the training and validation results from the network. It shows the working concept of the proposed structure, by incorporating the fluid flow characteristics and airfoil geometry. Results are presented at the epoch number with the lowest validation error.

Model validation {#Model}
----------------

The Absolute percent error (APE) or the unsigned percentage error is used as a metric for comparison: $$\begin{aligned}
\text{APE}=\frac{|\text{Prediction} - \text{Truth} |} {|\text{Truth}|} \times 100.
\label{eq:APE}\end{aligned}$$

The mean value of the absolute percent error (MAPE) is standard as a Loss function for regression problems. Here, model evaluation is done using MAPE due to the very intuitive interpretation regarding the relative error and its ease of use.

In this paper, the MAPE between the prediction and the truth is calculated in the wake region of an airfoil and the entire flow field around the airfoil. Here, the wake region of the airfoil is an area defined as $\{(x,y)|x\in\left[1.1, 1.5\right],y\in \left[-0.5, 0.5\right]\}$, and $\{(x,y)|x\in\left[-0.5, 1.5\right],y\in \left[-0.5, 0.5\right]\}$ is the entire flow field area around the airfoil. The predictions contain $2-3\%$ of points with an error value greater than $100\%$, which are treated as outliers and not included in the reported errors.

Numerical simulations {#Experiments}
---------------------

### Angle of attack variation {#Performance1}

At a fixed Reynolds number ($Re=1\times10^6$) and fixed airfoil shape (S805), we consider simulations with angles of attack of one-degree increments from $\alpha=0^{\circ}$ to $\alpha=20^{\circ}$. By using this small set of data (21 data points), we train the network with 50 filters instead of the aforementioned 300 filters in each layer (see Sec. \[Encoder\] for more details). The total loss function comprises only an MSE and with no regularization during training. Thus the cost function over the training set is presented as, $$\begin{aligned}
\text{Cost} = \lambda_\text{MSE} \times \text{MSE},
\label{eq:cost1}\end{aligned}$$ where $\lambda_\text{MSE}$ is a user defined parameter (here it is $\lambda_\text{MSE}=1$).

After the network training is complete, testing is performed on four unseen angles of attacks, $\alpha=2.5^\circ,~7.5^\circ,~12.5^\circ,~\text{and}~19.5^\circ$ respectively. Figure \[fig:fig10\] shows the comparison between the network prediction and the actual observation from the CFD simulation for the x-component of the velocity field around the S805 airfoil at an angle of attack of $\alpha = 12.5^\circ$. A visual comparison shows that the prediction is in agreement with the truth.

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

Table. \[tab:tab1separate\] and \[tab:tab1shared\] present the MAPE calculated in the wake region and the entire flow field around the S805 airfoil (see Fig.\[fig:fig10\]), where the fluid flow characteristics are the angle of attack of $\alpha = 12.5^\circ$ and the Reynolds number of $1\times10^6$.

[![width 1pt]{}l![width 1pt]{}c|c|c![width 1pt]{}]{}   & $\mathbf{U}$ & $\mathbf{V}$ & $\mathbf{P}$\
**Error in the wake region** & 24.9% & 10.15% & 24.97%\
**Error in the entire flow** & 13.51% & 11.92% & 13.50%\

[![width 1pt]{}l![width 1pt]{}c|c|c![width 1pt]{}]{}   & $\mathbf{U}$ & $\mathbf{V}$ & $\mathbf{P}$\
**Error in the wake region** & 15.08% & 7.98% & 14.82%\
**Error in the entire flow** & 9.62% & 8.65% & 7.31%\

The results in table. \[tab:tab1separate\] and \[tab:tab1shared\], illustrate that the errors in the wake region are generally similar to the errors in the entire flow field. This trend is true not only for this case but also in subsequent experiments. Figure \[fig:fig11\] shows the comparison between the CFD result and the network prediction of the x-component velocity profile of the airfoil wake at $x=1.1$ (downstream location from the leading edge).

[0.5]{}

[0.5]{}

### Shape, angle of attack, and Reynolds number variation {#Shape}

We train the network using $85$ percent of the 252 RANS simulation data-sets, with the variation of the airfoil shape, angle of attack and Reynolds number. Every convolutional layer is composed of 300 convolutional filters (see Sec.\[Encoder\] for more details). The total loss function during training comprises an MSE loss function with the L2 regularization. Thus, the cost function over the training set is presented as, $$\begin{aligned}
\text{Cost} = \lambda_\text{MSE} \times \text{MSE} + \lambda_\text{L2} \times \text{L2}_\text{regularization},
\label{eq:cost2}\end{aligned}$$ where $\lambda_\text{MSE}=1$ and $\lambda_\text{L2}=10^{-5}$ are user defined parameters.

Figures. \[fig:fig102\] and \[fig:fig103\] present the comparisons between the network predictions and observations for the x-component of the velocity field around the S809 and S814 airfoils at $(\alpha = 1^\circ,~Re = 1\times10^6)$ and $(\alpha = 19^\circ,~Re = 3\times10^6)$.

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

Quantitative results are presented in Tables \[tab:tab2separate\] and \[tab:tab2shared\].

[![width 1pt]{}c|c|c|c|c|c![width 1pt]{}]{} **Airfoil** & **AOA** & **Re** & **Variable** & **Error in the** & **Error in the**\
  &   & $\times 10^6$ &   & **wake region** & **entire flow**\
**S809** & $1^\circ$ & 1 & U & 12.25% & 10.35%\
**S809** & $1^\circ$ & 1 & V & 24.27% & 11.53%\
**S809** & $1^\circ$ & 1 & P & 5.14% & 8.40%\
**S814** & $19^\circ$ & 3 & U & 30.80% & 13.13%\
**S814** & $19^\circ$ & 3 & V & 10.43% & 5.49%\
**S814** & $19^\circ$ & 3 & P & 13.84% & 5.70%\

[![width 1pt]{}c|c|c|c|c|c![width 1pt]{}]{} **Airfoil** & **AOA** & **Re** & **Variable** & **Error in the** & **Error in the**\
  &   & $\times 10^6$ &   & **wake region** & **entire flow**\
**S809** & $1^\circ$ & 1 & U & 11.43% & 7.79%\
**S809** & $1^\circ$ & 1 & V & 15.53% & 8.74%\
**S809** & $1^\circ$ & 1 & P & 5.76% & 7.36%\
**S814** & $19^\circ$ & 3 & U & 27.23% & 13.20%\
**S814** & $19^\circ$ & 3 & V & 5.57% & 4.69%\
**S814** & $19^\circ$ & 3 & P & 12.93% & 5.71%\

### Shape, angle of attack, and Reynolds number variation with gradient sharpening {#Gradient}

To penalize the difference of the gradient in the loss function, and to address the lack of sharpness in predictions, we use gradient sharpening (GS) [@Mathieu:2015; @Lee:2018] in the loss functions combination and present the cost function over the training set as, $$\begin{aligned}
\text{Cost} = \lambda_\text{MSE} \times \text{MSE} + \lambda_\text{GS} \times \text{GS} +\lambda_\text{L2} \times \text{L2}_\text{regularization},
\label{eq:cost3}\end{aligned}$$ where $\lambda_\text{MSE},~\lambda_\text{GS}~\text{and}~\lambda_\text{L2}$ are the user defined parameters and their values are set via systematic experimentation, as $0.9,~0.1~\text{and}~10^{-5}$ respectively.

Figures. \[fig:fig104\] and \[fig:fig105\] present the comparisons between the network predictions with and without GS loss for the x-component of the velocity field around S809 and S814 airfoils respectively.

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

Visual comparisons of the predictions and the absolute difference with and without GS as illustrated in Figs. \[fig:fig104\] and \[fig:fig105\] are proofs of further gains and sharpness in the network predictions. The “absolute difference”\[absolute\] between the prediction and ground truth, for example, is defined as the absolute difference in the subtraction of each element in prediction from the corresponding element in ground truth. The MAPE for the components of the velocity field and pressure of the airfoils (S809 and S814 discussed above) are presented in table. \[tab:tab3separate\] and \[tab:tab3shared\]. The errors are reported in the wake region and the entire flow field around the airfoils with and without GS.

The predictions with GS in the loss function compared to not having it show significantly reduced errors in the wake region of the airfoil (twenty percent or more in the x-component of the velocity and pressure predictions) and obvious gains and sharpness in the entire flow field around the airfoil.

To further compare the accuracy of the network predictions, we use three probes around different airfoils in different flow conditions. These probes are leading edge probe (LE), trailing-edge probe (TE), and the probe at the wake region of an airfoil. Figure \[fig:fig200\] illustrates these three probes around different airfoils, S805, S809, and S814, respectively.

[0.45]{}

[0.45]{}

[0.45]{}

Table. \[tab:tab4\] presents the APE (Eq. \[eq:APE\]) at the probe locations (LE, TE, and wake region probe).

Figures. \[fig:fig17\], \[fig:fig18\] and \[fig:fig19\] illustrate the flow-field predictions with gradient sharpening in the loss function and in comparison with the reference results from the OVERTURNS CFD code.

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

[0.32]{}

Figures. \[fig:fig500\] and \[fig:fig501\] illustrate the x-component velocity profile of the airfoil wake at $x=1.1$ (downstream location from the leading edge). These predictions include GS in the loss function.

As a further comparison of the network prediction accuracy, we consider the pressure distribution on the upper and lower boundaries. Figures. \[fig:fig22\], \[fig:fig23\], and \[fig:fig24\] depicts the Ground truth vs. Predictions of the normalized pressure using the standard score normalization along the surface of the S805, S809, and S814 airfoils respectively. It is noteworthy that the surface with a one-pixel gap adjacent to the airfoil surface is used to obtain the pressure values. This change is due to the masking of the airfoil as an input during the training.

Overall, results are in good agreement with the ground truth simulation results in the entire range of angles of attacks and Reynolds numbers for the three different airfoils.

### Prediction for unseen airfoil shapes {#Unseen}

To further explore the predictive ability and accuracy of the trained network, three unseen geometries are considered as shown in Figure \[fig:fig25\]). The first one, denoted by “new airfoil” is an averaged shape of S809 and S814 airfoils. In addition, the S807 and S819 airfoils are also considered.

Figures \[fig:fig26\], \[fig:fig27\], \[fig:fig28\] illustrate the prediction of the network on the unseen airfoils in comparison to CFD simulations.

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[0.45]{}

[![width 1pt]{}c|c|c|c![width 1pt]{}]{} **Airfoil** & **Variable** & **Error in the** & **Error in the**\
  &   & **wake region** & **entire flow**\
**New** & U & 10.55% & 6.41%\
**New** & V & 5.71% & 8.59%\
**New** & P & 4.41% & 5.27%\
**S807** & U & 11.35% & 8.4%\
**S807** & V & 2.2% & 8.8%\
**S807** & P & 6.0% & 7.6%\
**S819** & U & 13.3% & 10.3%\
**S819** & V & 2.4% & 8.6%\
**S819** & P & 5.5% & 7.5%\

Table. \[tab:tab6\] provides a quantification of the results, and suggests good generalization properties of the network to an unseen shape.

Conclusions and Future Work {#Conclusion}
===========================

A flexible approximation model based on convolutional neural networks was developed for efficient prediction of aerodynamic flow fields. Shared-encoding and decoding was used and found to be computationally more efficient compared to separated alternatives. The use of convolution operations, parameter sharing and robustness to noise using the gradient sharpening were shown to enhance predictive capabilities. The Reynolds number, angle of attack, and the shape of the airfoil in the form of a signed distance function are used as inputs to the network and the outputs are the velocity and pressure fields.

The framework was utilized to predict the Reynolds Averaged Navier–Stokes flow field around different airfoil geometries under variable flow conditions. The network predictions on a single GPU were four orders of magnitude faster compared to the RANS solver, at mean square error levels of less than 10% over the entire flow field. Predictions were possible with a small number of training simulations, and accuracy improvements were demonstrated by employing gradient sharpening. Furthermore, the capability of the network was evaluated for unseen airfoil shapes.

The results illustrate that the CNNs can enable near real-time simulation-based design and optimization, opening avenues for an efficient design process. It is noteworthy that using three airfoil shapes in training, is a data limitation and reduces the general prediction behavior for unseen airfoil geometries from other families. Future work will seek to use a rich data set including multiple airfoil families in training and to augment the training data-sets to convert a set of input data into a broader set of slightly altered data [@Shijie:2017] using operations such as translation and rotation. This augmentation would effectively help the network from learning irrelevant patterns, and substantially boost the performance. Furthermore, exploring physical loss functions can be helpful in explicitly imposing physical constraints such as the conservation of mass and momentum to the networks.

Acknowledgements {#acknowledgements .unnumbered}
================

This work was supported by General Motors Corporation under a contract titled “Deep Learning and Reduced Order Modeling for Automotive Aerodynamics.” Computing resources were provided by the NSF via grant 1531752 MRI: Acquisition of Conflux, A Novel Platform for Data-Driven Computational Physics (Tech. Monitor: Stefan Robila).

Appendix: Governing equations {#appendix-governing-equations .unnumbered}
=============================

The RANS equations are derived by ensemble-averaging the conservation equations of mass, momentum and energy. These equations, for compressible flow are given by: $$\begin{aligned}
&\frac{\partial\bar\rho}{\partial t}+\frac{\partial\left(\bar\rho\hat u_i\right)}{\partial x_i} =0 \\
&\frac{\partial\left(\bar\rho\hat u_i\right)}{\partial t}+\frac{\partial\left(\bar\rho\hat u_i\hat u_j\right)}{\partial x_j} =-\frac{\partial\bar p}{\partial x_i}+\frac{\partial\bar\sigma_{ij}}{\partial x_j}+\frac{\partial\tau_{ij}}{\partial x_j} \\
&\frac{\partial\left(\bar\rho\hat E\right)}{\partial t}+\frac{\partial\left(\bar\rho\hat H\hat u_j\right)}{\partial x_j} = \frac{\partial}{\partial x_j}\left(\bar\sigma_{ij}\hat u_i+\overline{\sigma_{ij}u_i''}\right) - \\
&\nonumber \quad\quad \frac{\partial}{\partial x_j}\left(-\hat\kappa\frac{\partial\hat T}{\partial x_j}+c_P\overline{\rho u_j'' T''}-\hat u_i\tau_{ij}+\frac{1}{2}\overline{\rho u_i'' u_i'' u_j''}\right),\end{aligned}$$ where the overbar indicates conventional time-average mean, $u_i$ is the fluid velocity, $\rho$ is the density, $p$ is the pressure, $\tau_{ij}$ is the Reynolds stress term, $c_P$ is the heat capacity at constant pressure, and $\kappa$ is the kinetic energy of the fluctuating field (local turbulent kinetic energy). The density weighted time averaging (Favre averaging) of any quantity $\xi$, denoted by $\hat\xi$ is given as $\hat\xi=\overline{\rho\xi}/\bar\rho$, where, $$\begin{aligned}
&\hat H=\hat E+\frac{\bar p}{\bar\rho},\\
&\bar\sigma_{ij}=\mu_t\left(\frac{\partial\hat u_i}{\partial x_j}+\frac{\partial\hat u_j}{\partial x_i}-\frac{2}{3}\frac{\partial\hat u_k}{\partial x_k}\delta_{ij}\right),\\
&\tau_{ij}=-\overline{\rho u_i'' u_j''},\\
&k=\frac{\widehat{u_i''^2}+\widehat{v_i''^2}+\widehat{w_i''^2}}{2},\\
&\bar p = (\gamma-1)\bar\rho\left[\hat E - \frac{\hat u^2+\hat v^2+\hat w^2}{2} - k\right].\end{aligned}$$

To provide closure to the above equations, we use the model proposed by [@Spalart:1992]. In this closure, the Boussinesq hypothesis relates the Reynolds stress and the effect of turbulence as an eddy viscosity $\mu_t$. Employing the Boussinesq approach, and Reynolds Analogy a transport equation for a working variable $\tilde\nu$ is solved to estimate the eddy viscosity field at every iteration. $$\begin{aligned}
\label{eq:eq10}
&\frac{\partial \tilde{\nu}}{\partial t}+u_j\frac{\partial \tilde{\nu}}{\partial x_j} =  C_{b1}\left[1-f_{t2}\right]\tilde{S}\tilde{\nu} + \\
\nonumber &~\frac{1}{\sigma}\left\lbrace\nabla\cdot\left[\left(\nu+\tilde{\nu}\right)\nabla\tilde{\nu}\right]+C_{b2}\left|\nabla\tilde{\nu}\right|^2\right\rbrace
-\left[C_{w1}f_w-\frac{C_{b1}}{\kappa^2}f_{t2}\right]\left(\frac{\tilde{\nu}}{d}\right)^2.\end{aligned}$$ The turbulent eddy viscosity is computed as $\mu_t=\bar\rho\tilde\nu f_{v1}$, where, $$\begin{aligned}
&\nonumber f_{v1}=\frac{\chi^3}{\chi^3+C_{v1}^3},~\chi=\frac{\tilde{\nu}}{\nu},~\nu=\frac{\mu}{\bar\rho},\\
&\nonumber f_{t2}=C_{t3} \exp \left(-C_{t4}\chi^2\right),\\
&\nonumber \tilde{S}=S+\frac{\tilde{\nu}}{\kappa^2 d^2}f_{v2},\\
&\nonumber S=\sqrt{2\Omega_{ij}\Omega{ij}},~f_{v2}=1-\frac{\chi}{1+\chi f_{v1}},\\
&\nonumber f_w=g\left[\frac{1+C_{w3}^6}{g^6+C_{w3}^6}\right]^{1/6}, \\ 
&\nonumber g=r+C_{w2}(r^6-r),~
r=\frac{\tilde{\nu}}{\tilde{S}\kappa^2d^2}, \\
&\nonumber C_{w1}=\frac{C_{b1}}{\kappa^2}+\frac{1+C_{b2}}{\sigma},~ C_{b1}=0.1355,~\sigma=2/3,~C_{b2}=0.622,\\
&\nonumber \kappa=0.41,~C_{w2}=0.3,~
C_{w3}=2.0,~C_{v1}=7.1,~C_{t3}=1.2,~C_{t4}=0.5.\end{aligned}$$ The first term on the right hand side of this Eq. \[eq:eq10\] is the production term for $\tilde\nu$ while the second term represents dissipation.

[48]{} \[1\][\#1]{} \[1\][[\#1]{}]{} urlstyle \[1\][DOI \#1]{}

\[2\]\[\][[\#2](#2)]{}

Aksoy S, Haralick RM (2000) Feature normalization and likelihood-based similarity measures for image retrieval. Pattern Recognition Letters 22:563–582

Amidror I (2002) Scattered data interpolation methods for electronic imaging systems: a survey. J Electronic Imaging 11(2):157–176

A, [Lakshminarayan]{} V, [Duraisamy]{} K (2012) Assessment of transition model and cfd methodology for wind turbine flows. In: 42nd AIAA Fluid Dynamics Conference and Exhibit, American Institute of Aeronautics and Astronautics, p 2720

Bengio Y (2009) Learning deep architectures for ai. Foundations and Trends in Machine Learning 2(1):1–127

Carr JC, Beatson RK, Cherrie JB, Mitchell TJ, Fright WR, McCallum BC, Evans TR (2001) Reconstruction and representation of 3d objects with radial basis functions. In: Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, SIGGRAPH ’01, pp 67–76

S, [Tou Ng]{} H (2018) [A Multilayer Convolutional Encoder-Decoder Neural Network for Grammatical Error Correction]{}. arXiv e-prints arXiv:1801.08831,

Duraisamy K (2005) Studies in tip vortex formation, evolution and contro, dept of aerospace engineering, univ of maryland. PhD thesis, University of Maryland

Duraisamy K, Iaccarino G, Xiao H (2019) Turbulence modeling in the age of data. Annual Review of Fluid Mechanics 51(1):357–377

B, [Karaoglu]{} S, [Saha]{} SK (2015) [Object Class Detection and Classification using Multi Scale Gradient and Corner Point based Shape Descriptors]{}. arXiv e-prints arXiv:1505.00432,

Foley JD, van Dam A, Feiner SK, Hughes JF (1995) Computer Graphics: Principles and Practice in C, second edition edn. Addison-Wesley

Fuhrmann S, Goesele M (2014) Floating scale surface reconstruction. ACM Trans Graph 33(4):46:1–46:11

Guo X, Li W, Iorio F (2016) Convolutional neural networks for steady flow approximation. In: Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, ACM, KDD ’16, pp 481–490

O (2017) [Lat-Net: Compressing Lattice Boltzmann Flow Simulations using Deep Neural Networks]{}. arXiv e-prints arXiv:1705.09036,

Hoppe H, DeRose T, Duchamp T, McDonald J, Stuetzle W (1992) Surface reconstruction from unorganized points. SIGGRAPH Comput Graph 26(2):71–78

Jürgen S (2015) Deep learning in neural networks: An overview. Neural Networks 61:85–117

Kazhdan M, Bolitho M, Hoppe H (2006) [Poisson Surface Reconstruction]{}. In: Sheffer A, Polthier K (eds) Symposium on Geometry Processing, The Eurographics Association, pp 61–70

Koren B (1993) A robust upwind discretization method for advection, diffusion and source terms, Vieweg, pp 117–138

Lakshminarayan VK, Baeder JD (2010) Computational investigation of microscale coaxial-rotor aerodynamics in hover. Journal of Aircraft 47(3):940–955

Lecun Y, Bottou L, Bengio Y, Haffner P (1998) Gradient-based learning applied to document recognition. Proceedings of the IEEE 86(11):2278–2324

S, [You]{} D (2017) [Prediction of laminar vortex shedding over a cylinder using deep learning]{}. arXiv e-prints arXiv:1712.07854,

S, [You]{} D (2018) [Data-driven prediction of unsteady flow fields over a circular cylinder using deep learning]{}. arXiv e-prints arXiv:1804.06076,

van Leer B (1979) Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov’s method. Journal of Computational Physics 32(1):101–136

Ling H, Jacobs DW (2007) Shape classification using the inner-distance. IEEE Transactions on Pattern Analysis and Machine Intelligence 29(2):286–299

M, [Couprie]{} C, [LeCun]{} Y (2015) [Deep multi-scale video prediction beyond mean square error]{}. arXiv e-prints arXiv:1511.05440,

Medida S, Baeder J (2011) Application of the Correlation-based Gamma-Re Theta t Transition Model to the Spalart-Allmaras Turbulence Model, American Institute of Aeronautics and Astronautics

TP, [Jaiman]{} RK (2017) [An Efficient Deep Learning Technique for the Navier-Stokes Equations: Application to Unsteady Wake Flow Dynamics]{}. arXiv e-prints arXiv:1710.09099,

Pandya S, Venkateswaran S, Pulliam T (2003) Implementation of Preconditioned Dual-Time Procedures in OVERFLOW, American Institute of Aeronautics and Astronautics

L, [Bonev]{} B, [Thuerey]{} N (2017) [Generating Liquid Simulations with Deformation-aware Neural Networks]{}. arXiv e-prints arXiv:1704.07854,

Pulliam T, Chaussee D (1981) A diagonal form of an implicit approximate-factorization algorithm. Journal of Computational Physics 39(2):347–363

Raissi M, Karniadakis GE (2018) Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics 357:125–141

M, [Yazdani]{} A, [Karniadakis]{} GE (2018) [Hidden Fluid Mechanics: A Navier-Stokes Informed Deep Learning Framework for Assimilating Flow Visualization Data]{}. arXiv e-printse-prints p arXiv:1808.04327

Raissi M, Perdikaris P, Karniadakis G (2019) Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378:686–707

P, [Zoph]{} B, [Le]{} QV (2017) [Searching for Activation Functions]{}. arXiv e-prints arXiv:1710.05941,

Roe PL (1986) Characteristic-based schemes for the euler equations. Annual Review of Fluid Mechanics 18(1):337–365

Sethian JA (1996) A fast marching level set method for monotonically advancing fronts. Proc Natl Acad Sci USA 93:1591–1595

J, [Ping]{} W, [Peiyi]{} J, [Siping]{} H (2017) Research on data augmentation for image classification based on convolution neural networks. In: 2017 Chinese Automation Congress (CAC), pp 4165–4170

Somers D (1997) Design and experimental results for the s805 airfoil. Tech. Rep. NREL/SR-440-6917

Somers DM (1997) Design and experimental results for the s809 airfoil. Tech. Rep. NREL/SR-440-6918

Somers DM (2004) S814 and s815 airfoils: October 1991–july 1992. Tech. rep.

Spalart P, Allmaras S (1992) A one-equation turbulence model for aerodynamic flows, American Institute of Aeronautics and Astronautics

Taylor GW, Fergus R, LeCun Y, Bregler C (2010) Convolutional learning of spatio-temporal features. In: Computer Vision – ECCV 2010, pp 140–153

J, [Schlachter]{} K, [Sprechmann]{} P, [Perlin]{} K (2016) [Accelerating Eulerian Fluid Simulation With Convolutional Networks]{}. arXiv e-prints arXiv:1607.03597,

Turkel E (1999) Preconditioning techniques in computational fluid dynamics. Annual Review of Fluid Mechanics 31(1):385–416

K, [Kim]{} VG, [Huang]{} Q, [Kalogerakis]{} E (2015) [Data-Driven Shape Analysis and Processing]{}. arXiv e-prints arXiv:1502.06686,

Zhang D, Lu G (2004) Review of shape representation and description techniques. Pattern Recognition 37(1):1–19

Y, [Sung]{} WJ, [Mavris]{} D (2017) [Application of Convolutional Neural Network to Predict Airfoil Lift Coefficient]{}. arXiv e-prints arXiv:1712.10082,

Zhao H (2005) A fast sweeping method for [E]{}ikonal equations. Math Comput 74(250):603–627

Zuo Z, Shuai B, Wang G, Liu X, Wang X, Wang B, Chen Y (2015) Convolutional recurrent neural networks: Learning spatial dependencies for image representation. In: 2015 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pp 18–26
