Prediction of Aerodynamic Flow Fields Using Convolutional Neural Networks

TL;DR

This paper introduces a CNN-based model that rapidly predicts flow fields around airfoils, significantly reducing computation time compared to traditional RANS simulations, enabling near-real-time aerodynamic analysis.

Contribution

The study presents a novel CNN architecture tailored for flow field prediction over airfoils, demonstrating high accuracy and speed in diverse conditions, surpassing conventional methods.

Findings

CNN predicts velocity and pressure fields orders of magnitude faster than RANS.

The model effectively generalizes across different airfoil shapes and flow conditions.

Robustness to noise and specific convolution techniques enhance prediction accuracy.

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.

Tables

Table 1. Table 1 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure in the wake region of the S805 airfoil and the entire flow field around it (separated decoder Fig.  3(a) ).

	$𝐔$	$𝐕$	$𝐏$
Error in the wake region	24.9%	10.15%	24.97%
Error in the entire flow	13.51%	11.92%	13.50%

Table 2. Table 2 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure in the wake region of the S805 airfoil and the entire flow field around it (shared decoder Fig.  3(b) ).

	$𝐔$	$𝐕$	$𝐏$
Error in the wake region	15.08%	7.98%	14.82%
Error in the entire flow	9.62%	8.65%	7.31%

Table 3. Table 3 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure (shared decoder).

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%

Table 4. Table 4 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure (separated decoder).

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%

Table 5. Table 5 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure (separated decoder).

Airfoil	AOA	Re	Variable	Error in the	Error in the	Error in the	Error in the
				wake region	wake region	entire flow	entire flow
		$\times {10}^6$		without GS	with GS	without GS	with GS
S809	$1^{\circ }$	1	U	11.43%	11.30 %	7.79%	9.92%
S809	$1^{\circ }$	1	V	15.53%	19.43 %	8.74%	8.88 %
S809	$1^{\circ }$	1	P	5.76%	6.82 %	7.36%	8.09 %
S814	${19}^{\circ }$	3	U	27.23%	9.49 %	13.20%	7.38 %
S814	${19}^{\circ }$	3	V	5.57%	3.09%	4.69%	2.94 %
S814	${19}^{\circ }$	3	P	12.93%	4.71%	5.71%	3.11 %
S805	$4^{\circ }$	3	U	13.40%	4.47%	6.44%	2.99 %
S805	$4^{\circ }$	3	V	3.96%	1.42%	5.60%	4.85 %
S805	$4^{\circ }$	3	P	8.48%	3.72%	5.37%	3.22 %

Table 6. Table 6 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure with and without GS (shared decoder).

Airfoil	AOA	Re	Variable	Error in the	Error in the	Error in the	Error in the
				wake region	wake region	entire flow	entire flow
		$\times {10}^6$		without GS	with GS	without GS	with GS
S809	$1^{\circ }$	1	U	12.25%	8.35 %	10.35%	5.57%
S809	$1^{\circ }$	1	V	24.27%	17.43 %	11.53%	6.52 %
S809	$1^{\circ }$	1	P	5.14%	3.18 %	8.40%	3.86 %
S814	${19}^{\circ }$	3	U	30.80%	19.89 %	13.13%	11.37 %
S814	${19}^{\circ }$	3	V	10.43%	3.52%	5.49%	3.04 %
S814	${19}^{\circ }$	3	P	13.84%	7.95%	5.70%	3.57 %
S805	$4^{\circ }$	3	U	15.68%	15.15%	6.50%	9.59 %
S805	$4^{\circ }$	3	V	4.20%	1.29%	7.63%	10.70 %
S805	$4^{\circ }$	3	P	8.58%	32.77%	4.03%	15.06 %

Table 7. Table 7 : APE at probe locations (shared decoder).

Airfoil	AOA	Re	Variable	Error at	Error at	Error at
		$\times {10}^6$		LE probe	TE probe	wake probe
S805	${12}^{\circ }$	3	U	2.36%	1.67%	31.66%
S805	${12}^{\circ }$	3	V	0.96%	0.74%	1.50%
S805	${12}^{\circ }$	3	P	9.35%	6.56%	22.44%
S809	$1^{\circ }$	3	U	0.08%	2.38%	5.45%
S809	$1^{\circ }$	3	V	0.73%	2.22%	1.93%
S809	$1^{\circ }$	3	P	0.23%	3.16%	1.07%
S814	${15}^{\circ }$	2	U	0.75%	1.98%	1.64%
S814	${15}^{\circ }$	2	V	0.76%	1.17%	0.18%
S814	${15}^{\circ }$	2	P	5.67%	6.56%	8.02%

Table 8. Table 8 : MAPE for the components of the velocity field ( U 𝑈 U and V 𝑉 V respectively) and pressure in the wake region of S807 and S819 airfoils respectively and the entire flow field around them.

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%

Equations

$$\text{SDF}(\text{x})=\{\begin{array}{cc}d(\text{x},\partial \Omega )& \text{x}\notin \Omega \\ 0& \text{x}\in \partial \Omega \\ -d(\text{x},\partial \Omega )& \text{x}\in \Omega \end{array},$$

$${\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}$$

$$\displaystyle\quad\big{[}(U^{l}_{{ij}_{\text{truth}}}-U^{l}_{{ij}_{\text{pred}}})^{2}~{}+~{}(V^{l}_{{ij}_{\text{truth}}}-V^{l}_{{ij}_{\text{pred}}})^{2}~{}+~{}$$

$$\displaystyle\quad~{}(P^{l}_{{ij}_{\text{truth}}}-P^{l}_{{ij}_{\text{pred}}})^{2}\big{]},$$

$${\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}$$

$$[({\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+$$

$$({\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+$$

$$({\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],$$

$${\text{L2}}_{\text{regularization}}=\frac{1}{2m}\sum _{l=1}^L\sum _{i=1}^{n_l}({\theta }_i^l{)}^2,$$

$${\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}$$

$$\displaystyle\quad\big{[}(X^{l}_{{ij}_{\text{truth}}}-X^{l}_{{ij}_{\text{pred}}})^{2}],$$

$${\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}$$

$$[({\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],$$

$$\text{APE}=\frac{|\text{Prediction}-\text{Truth}|}{|\text{Truth}|}\times 100.$$

$$\text{Cost}={\lambda }_{\text{MSE}}\times \text{MSE},$$

$$\text{Cost}={\lambda }_{\text{MSE}}\times \text{MSE}+{\lambda }_{\text{L2}}\times {\text{L2}}_{\text{regularization}},$$

$$\text{Cost}={\lambda }_{\text{MSE}}\times \text{MSE}+{\lambda }_{\text{GS}}\times \text{GS}+{\lambda }_{\text{L2}}\times {\text{L2}}_{\text{regularization}},$$

$$\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^{\prime \prime }}\right)-$$

$$\frac{\partial }{\partial x_j}\left(-\hat{\kappa }\frac{\partial \hat{T}}{\partial x_j}+c_P\overline{\rho u_j^{\prime \prime }{T}^{\prime \prime }}-{\hat{u}}_i{\tau }_{ij}+\frac{1}{2}\overline{\rho u_i^{\prime \prime }{u}_i^{\prime \prime }{u}_j^{\prime \prime }}\right),$$

$$\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^{\prime \prime }{u}_j^{\prime \prime }},$$

$$k=\frac{\widehat{u_i^{\prime \prime 2}}+\widehat{v_i^{\prime \prime 2}}+\widehat{w_i^{\prime \prime 2}}}{2},$$

$$\bar{p}=(\gamma -1)\bar{\rho }\left[\hat{E}-\frac{{\hat{u}}^2+{\hat{v}}^2+{\hat{w}}^2}{2}-k\right].$$

$$\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 }+$$

$$\frac{1}{\sigma }\left\{\nabla \cdot \left[\left(\nu +\tilde{\nu }\right)\nabla \tilde{\nu }\right]+C_{b2}{\left|\nabla \tilde{\nu }\right|}^2\right\}-\left[C_{w1}{f}_w-\frac{C_{b1}}{{\kappa }^2}{f}_{t2}\right]{\left(\frac{\tilde{\nu }}{d}\right)}^2.$$

$$f_{v1}=\frac{{\chi }^3}{{\chi }^3+C_{v1}^3},\chi =\frac{\tilde{\nu }}{\nu },\nu =\frac{\mu }{\bar{\rho }},$$

$$f_{t2}=C_{t3}\exp \left(-C_{t4}{\chi }^2\right),$$

$$\tilde{S}=S+\frac{\tilde{\nu }}{{\kappa }^2{d}^2}{f}_{v2},$$

$$S=\sqrt{2{\Omega }_{ij}\Omega ij},f_{v2}=1-\frac{\chi }{1+\chi f_{v1}},$$

$$f_w=g{\left[\frac{1+C_{w3}^6}{g^6+C_{w3}^6}\right]}^{1/6},$$

$$g=r+C_{w2}(r^6-r),r=\frac{\tilde{\nu }}{\tilde{S}{\kappa }^2{d}^2},$$

$$C_{w1}=\frac{C_{b1}}{{\kappa }^2}+\frac{1+C_{b2}}{\sigma },C_{b1}=0.1355,\sigma =2/3,C_{b2}=0.622,$$

$$\kappa =0.41,C_{w2}=0.3,C_{w3}=2.0,C_{v1}=7.1,C_{t3}=1.2,C_{t4}=\mathrm{0.5.}$$

Full text

††journal: Computational Mechanics

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

1.1 State of the art in application of CNNs in fluid dynamics

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

2 Methodology

2.1 CFD Simulation

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

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

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

2.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$.

[TABLE]

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.

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

2.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).

2.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:

[TABLE]

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:

[TABLE]

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.

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

3.1 Model validation

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

[TABLE]

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.

3.2 Numerical simulations

3.2.1 Angle of attack variation

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,

[TABLE]

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

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

[TABLE]

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.

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

[TABLE]

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.

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

4 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:

[TABLE]

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,

[TABLE]

To provide closure to the above equations, we use the model proposed by Spalart and Allmaras (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.

[TABLE]

The turbulent eddy viscosity is computed as $\mu_{t}=\bar{\rho}\tilde{\nu}f_{v1}$, where,

[TABLE]

The first term on the right hand side of this Eq. 19 is the production term for $\tilde{\nu}$ while the second term represents dissipation.