Prediction of Aerodynamic Flow Fields Using Convolutional Neural Networks

Introduction

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 et al 2019) or replace (Guo et al 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; Jürgen 2015) have shown significant successes in learning from data, and have been successfully used in the development of novel computational approaches (Raissi et al 2018; Raissi and Karniadakis 2018; Raissi et al 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 et al 1998; Taylor et al 2010; Zuo et al 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 1 illustrates the simplified layout of a typical CNN, LeNet-5 (Lecun et al 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.

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

Guo et al 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 and Jaiman 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 and You 2017; Lee and You 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 2017 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 et al 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 et al 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 et al 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 et al 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

In this work, flow computations and analyses are performed using the OVERTURNS CFD code (Duraisamy 2005; Lakshminarayan and Baeder 2010). This code solves the compressible RANS equations using a preconditioned dual-time scheme (Pandya et al 2003). Iterative solutions are pursued using the implicit approximate factorization method (Pulliam and Chaussee 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) (van Leer 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 and Allmaras 1992) turbulence model and $\gamma-\overline{Re_{\theta t}}$ model (Medida and Baeder 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 et al 2012). Computations are performed using structured C-meshes with dimensions $394\times 124$ in the wrap-around and normal directions respectively. Figure 2 shows the airfoils and their near-body meshes.

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 et al 2012).

2 Convolutional Neural Networks

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.

3 Network 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 et al 2016. Encoder-decoder CNNs are most widely used for machine translation from a source language to a target language (Chollampatt and Tou Ng 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 3 illustrates the proposed CNN architecture designed in this work.

Guo et al 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 et al 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.

4 Geometry Representation

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 et al 1992; Carr et al 2001; Kazhdan et al 2006; Fuhrmann and Goesele 2014), or shape representation and classification (Zhang and Lu 2004; Ling and Jacobs 2007; Xu et al 2015; Fernando et al 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 et al 2016; Prantl et al 2017), we use the SDF sampled on a Cartesian grid. Guo et al 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.

5 Signed Distance Function

A mathematical definition of the signed distance function of a set of points X determines the minimum distance of each given point $\textbf{x}\in\textbf{X}$ from the boundary of an object $\partial\Omega$.

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 x from the object boundary points. The distance sign determines whether the given point is inside or outside of the object. Figure 4 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 5 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 et al 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.

6 Convolutional Encoder-Decoder Approach

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 et al 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.

7 Data Preparation

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. 6.

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 7 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\times 10^{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. 8 and the CFD data in Fig. 7 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 and Haralick 2000).

8 Network Training and Hyper-parameter Study

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. 2.6). 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:

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 et al 2015; Lee and You 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:

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

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 9 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.

The Absolute percent error (APE) or the unsigned percentage error is used as a metric for comparison:

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.

2 Numerical simulations

At a fixed Reynolds number ($Re=1\times 10^{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. 2.6 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,

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 10 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.

Table. 1 and 2 present the MAPE calculated in the wake region and the entire flow field around the S805 airfoil (see Fig.10), where the fluid flow characteristics are the angle of attack of $\alpha=12.5^{\circ}$ and the Reynolds number of $1\times 10^{6}$.

The results in table. 1 and 2, 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 11 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).

2.2 Shape, angle of attack, and Reynolds number variation

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.2.6 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,

where $\lambda_{\text{MSE}}=1$ and $\lambda_{\text{L2}}=10^{-5}$ are user defined parameters.

Figures. 12 and 13 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\times 10^{6})$ and $(\alpha=19^{\circ},~Re=3\times 10^{6})$.

Quantitative results are presented in Tables 3 and 4.

2.3 Shape, angle of attack, and Reynolds number variation with gradient sharpening

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 et al 2015; Lee and You 2018) in the loss functions combination and present the cost function over the training set as,

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. 14 and 15 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.

Visual comparisons of the predictions and the absolute difference with and without GS as illustrated in Figs. 14 and 15 are proofs of further gains and sharpness in the network predictions. The “absolute difference” 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. 5 and 6. 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 16 illustrates these three probes around different airfoils, S805, S809, and S814, respectively.

Table. 7 presents the APE (Eq. 7) at the probe locations (LE, TE, and wake region probe).

Figures. 17, 18 and 19 illustrate the flow-field predictions with gradient sharpening in the loss function and in comparison with the reference results from the OVERTURNS CFD code.

Figures. 20 and 21 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. 22, 23, and 24 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.

2.4 Prediction for unseen airfoil shapes

To further explore the predictive ability and accuracy of the trained network, three unseen geometries are considered as shown in Figure 25). 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 26, 27, 28 illustrate the prediction of the network on the unseen airfoils in comparison to CFD simulations.

Table. 8 provides a quantification of the results, and suggests good generalization properties of the network to an unseen shape.

Conclusions and Future Work

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 et al 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

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

The RANS equations are derived by ensemble-averaging the conservation equations of mass, momentum and energy. These equations, for compressible flow are given by:

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,

References