1. Introduction
Surrogate models play an important role in product development, optimization problems, real time processing, or sub modelling. In many cases, surrogate models are required to replace CFD simulations due to their high computational costs. A surrogate model can be built via a fit to sample points that can either be simulated or measured. Prevalent models comprise linear regression, Bayesian models, random forests, and neural networks (NN), among others [
1].
According to the approximation theorem of Pinkus, a single hidden layer NN with enough neurons can approximate any function and its derivatives [
2]. This feature makes NN suitable for surrogate modeling as it avoids limitations in regression accuracy to complex functions that are related to other surrogate methods, such as linear regression. Consequently, several investigations show a superiority of NN when compared to other regression models [
3,
4,
5]. A vast variety of different architectures and algorithms of NN have been developed and applied, comprising feedforward NN, recurrent NN, convolutional NN, reinforcement learning, or generative adversarial networks, among others [
6]. A drawback of NNs is the need for large high-quality datasets in the training process [
7].
Recently, physics-informed neural networks (PINN) have become a focus of research activities. A PINN is a NN that is trained to respect given laws of physics, typically partial differential equations (PDE), as well as boundary conditions (BC). In the training process, PINNs create sample data by solving the PDE under the given BC. Hence, a lack of simulation or measurement data can be compensated. Since flow measurements and CFD simulations are often time consuming and expensive, the use of PINNs could improve and shorten the surrogate modelling process.
The concept of PINN has first been described by Lagaris et al. [
8] who showed the applicability of PINNs to solve several ordinary differential equations (ODEs), coupled ODEs and PDEs. Raissi et al. [
9,
10,
11] solved the Burgers equation, the Schrödinger equation, and the Allen-Cahn equation with a PINN for time dependent one-dimensional problems. The results exhibited excellent agreement between the PINN and the exact solution. The publications presented above demonstrated the suitability of PINNs to solve PDEs.
PINNs were also applied to the field of fluid dynamics to solve the Navier-Stokes (NS) equations. Here, investigations cover the application of feedforward NNs to two- and three dimensional problems [
12], near wall flows [
13], thermal flows [
14,
15], cavitation flows [
16], reconstruction of incomplete data [
17,
18], boundary condition enforcement and adaptive activation functions [
19], convolutional PINNs [
20,
21], and different fomulations of the NS equations [
22], among others.
In most publications PINNs were used to predict flow fields for Reynolds numbers below 300 [
12,
14,
19,
22,
23]. However, many technical problems are subject to high Reynolds numbers and corresponding turbulent flows. Turbulent flows result in highly complex flow fields in both space, and time. Hence, solving the Naver-Stokes equations with a PINN for highly turbulent flows can be a challenging task. As for many applications only stationary information is required, PINNs predicting the Reynolds-averaged Navier-Stokes (RANS) equations incorporating a turbulence model are suitable. A first approach has been done by Hennigh et al. [
24] who applied the mixing length turbulence model. Eivazi et al. [
25] trained a PINN to statisfy the RANS equations without a transport-equation-based model to predict a turbulent boundary layer, the high Reynolds flow in a sub domain on an airfoil, as well as a flow over a periodic hill. Xu et al. [
26] used a similar approach and trained a PINN to predict the turbulent viscosity without any imposed equation. They trained the model to CFD data and used it to explore missing flow data inside of a finite area downstream of a backward-facing step. Pioch et al. [
27] applied the mixing length model, the
k-
model, an equation-free pseudo-Reynolds stress model and an equation-free
model to a backward facing step flow. A favourable agreement with data from a direct numerical simulation was found when three or five lines of labeled training data were used.
In the research presented above, PINNs were trained to predict a two- or three-dimensional solution for a single geometry. Wandel et al. presented a convolutional PINN that takes the vector potential as an input. The model was trained with different 3D shapes, such as boxes, cylinders, or spinning balls. They reported plausible flow fields for geometries not contained in the training dataset, such as a fish or multiple boxes [
28]. Hennigh et al. presented the NVIDIA SimNet framework that can read standard tesselation language (stl) files and comes with the zero-equation mixing length turbulence model. They applied the framework on the NVSwitch with parameterized geometry and used a modified Fourier feature PINN to find the optimal parameter configuration [
24]. Arthurs and King presented an active learning algorithm and trained a PINN to predict the flow and pressure field in a tube for two shape parameters and reported low errors [
29]. These investigations demonstrate the general suitability of PINNs as surrogate models for flows around variable geometries. To further assess the capabilities and limitations of PINNs as multidimensional surrogate models of variable geometries, investigations of complex flows at high Reynolds numbers are needed.
In this work, we contribute to the question how multidimensional surrogate models of Reynolds-averaged flows could be generated using a PINN and how accurately the flow field is predicted. As an example, we present a PINN-based surrogate model of the DU99W350 airfoil at a high Reynolds number of
using the Python library DeepXDE [
30]. The model takes the two-dimensional coordinates as well as the angle of attack as input parameters and predicts the mean pressure and velocity fields around the airfoil. To investigate the applicability of PINNs to complex fluid dynamics, the model was designed to predict the developing flow separation on the suction surface. In the training process, the RANS equations were solved. We tested the standard
k-
model by Wilcox [
31] and a pseudo-Reynolds stress turbulence modeling approach [
25,
27]. We evaluate the accuracy of the model when compared to CFD simulations and further discuss the limitations and potential improvements of the presented PINN methodology. The scope of this work was to investigate the capability of a simple feed-forward network architecture to serve as a surrogate model of a complex Reynolds-averaged flow around a parameterized geometry at an elevated Reynolds number.
3. Methodology
To yield the multidimensional surrogate model, a feedforward NN was trained to predict the Reynolds-averaged velocity and pressure field. The input layer of the NN consisted of three neurons that represent the three input dimensions coordinates x, y, and angle of attack (AoA)
. In case the neural network is trained without any turbulence model, the output layer comprises three neurons and displays the Reynolds-averaged predictions of u, v, and p. The PINN then reads:
where
are the trainable parameters (weights and biases) of the NN.
Figure 2 shows the corresponding architecture. As sketched, the output of the NN is fed into three different losses that build a so-called composed loss function.
In the composed loss function, losses for boundary conditions, simulated training values, and a system of partial differential equations (PDEs) were added. The composed loss function
reads:
with the adjustable network weights and biases
, the loss on the boundary conditions
, the loss for the training data
, and the loss for the partial differential equations
. The individual loss terms were calculated using the mean squared error:
where
,
, and
represent the number of points for which the boundary conditions, labeled training data, and PDEs are trained.
and
are the given data for point
n on the boundaries and for the labeled data coordinates, respectively.
and
represent the output of the PINN at the corresponding points and
is the residual of the
kth equation at point
n.
In this work, the continuity equation and the two-dimensional Reynolds-averaged Navier-Stokes equations for an incompressible fluid were used as PDEs. The corresponding residuals read:
with Reynolds-averaged velocities and pressure
,
, and
, constant density
, and constant kinematic viscosity
. The RANS equations contain the Reynolds stress terms:
Here,
represents the Reynolds stresses. For traditional turbulence modeling as applied in computational fluid dynamics (CFD), the Reynolds stresses are estimated using the Boussinesq hypothesis:
where
k is the turbulent kinetic energy,
is the Kronecker delta, and
is the turbulent viscosity. This approach is based on the gradient diffusion hypothesis which states that the turbulent transport is aligned with the negative gradient of the mean flow. The isotropic turbulent viscosity
then serves as a diffusion coefficient, supposing an analogy between viscous stresses and turbulent stresses.
In this work, we considered two different turbulence-modeling approaches. The first model is the standard
k-
model of Wilcox [
31]. This model is frequently applied for CFD calculations. For the
k-
model, the turbulent viscosity is calculated using:
with the specific dissipation rate
. The 2D stationary equations for
k and
read:
with
, as recommended by Wilcox [
31]. For the composed loss function, the transport equations for
k and
are to be included in the system of PDEs. The corresponding PINN reads:
The second considered model is a pseudo-Reynolds stress model as applied by Eivazi et al. [
25] and Pioch et al. [
27]. For this model, the Reynolds stress terms are calculated without any specific modeling using only two pseudo-turbulent velocity fluctuations,
and
:
Using this approach, in the training process the NN is penalized to learn an averaged flow field as well as pseudo-Reynolds variables that together satisfy the RANS equations. The corresponding PINN reads:
To calculate the losses associated with the boundary conditions, a slip wall condition was defined for the upper and lower channel walls. For the airfoil surface, we used a no slip wall condition. At the left boundary of the channel, we deployed a velocity inlet boundary condition. No pressure outlet condition was defined.
Setup and training of the multidimensional PINN was done with a custom Python code using the tensorflow-based library DeepXDE [
30]. The model was trained with the described RANS equations and boundary conditions as well as simulation data. We used a PointSetBC to train the labeled reference data and anchor points to train the PDEs.
To define the training points inside the multidimensional training domain, the coordinates of two different CFD grids were used. The grids were generated with for different AoA. To obtain the labeled training data, two-dimensional simulations were carried out with Comsol Multiphysics 5.6 on a hybrid mesh consisting of unstructured triangular finite elements and layers of quadrilateral elements to account for the boundary layer. The mesh featured 15 inflation layers with a growth rate of 1.1, a maximal base size of 0.01 m with a maximal growth rate of 1.08, as well as a grid refinement on the airfoil surface with a maximum cell size of 0.002 m. The trial functions were 2nd order accurate. Further, a hybrid wall function, the k- turbulence model and a convergence limit of were defined. The maximal value was 2.9. The pressure values were downscaled by a factor of . To train the PINN, a batch of random data points was drawn from the complete labeled data set ahead of the training process. The coordinates for the training of the PDEs were derived from a coarser grid without any further selection. Some illustrations concerning the training coordinates are given in the Appendix. The advantage of this method is the concentration of training points in important regions such as the boundary layer. The training of the PDEs as well as the boundary conditions was conducted for AoA of 10.0°, 12.5°, 15.0°, and 17.5°. Labeled training data were considered for AoA of 10.0 and 15.0 only.
Figure 3 illustrates the applied boundary conditions as well as the considered training data. As seen, the training domain was built by a multidimensional hypercube and the training of the PINN was conducted at discrete AoA which represent slices of the hypercube. This procedure facilitated the generation of training coordinates for specific AoA using a CFD meshing algorithm. The PINN could then be trained on these coordinates using the PointSet boundary conditions for the labeled data and the anchor points for the PDEs. This method was applied as the automatic generation of appropriately positioned training coordinates for n-dimensional problems is challenging. In this case, definition of the n-dimensional geometry as well as a suitable training point coordinate generator were necessary. The method applied here allows, in contrast, for point generation using available CFD meshers discretizing two- or three-dimensional geometries which can be obtained using conventional computer aided design (CAD) software. If we assume that a three-dimensional geometry then was to be trained for one parameter in the fourth dimension, such as AoA, then the training points could be generated for the geometry at discrete AoA.
For the training of the neural network, the loss terms were weighted using:
with
,
, and
as constant weight factors for boundary conditions, labeled training data, and PDEs, respectively. In the training procedure, the weights and biases
of the NN were adjusted using the optimizers Adam and L-BFGS to minimize the weighted composed loss function. The model was trained for 30,000 epochs using the Adam optimizer. Afterwards, the model was further trained with the L-BFGS optimizer under the predefined default settings. In the Appendix, a supplementary illustration concerning the corresponding losses for PDEs, BCs, and labeled data is given.
Table 1 lists the relevant training parameters.
We trained a NN consisting of five hidden layers with 128 neurons each. The training was repeated ten times to account for the random initialization of the weights and biases. For the further evaluation, the model with the lowest normalized mean squared error (NMSE) was selected.
4. Sensitivity analyses
In this section, results of several preliminary sensitivity analyses are presented. For these investigations, we limited the training to an AoA of 15.0°. The number of training points , , and were 800, 1600, and 1558, respectively. The scope of this section was to identify important training settings and to evaluate the accuracy of the corresponding predictions. Here, we present evaluations of the Reynolds number and the corresponding scaling method, the amount and distribution of training data, and a comparison between the k- model and the pseudo-Reynolds stress model.
For the investigations concerning the Reynolds numbers, we tested the accuracy of the training for
. These Reynolds numbers were obtained following two different procedures. For the first approach, a fluid with a density
of 1000
and a dynamic viscosity
of 0.001
was defined. The inlet velocity u was then varied between and 0.0005 m/s and 5 m/s. For the second approach, the inlet velocity was set to 1 m/s and the viscosity was varied between 0.2
and 0.0002
. The results of this analysis are shown in
Figure 4. The training runs with different inlet velocities exhibited severe deviations from the reference solutions for lower Reynolds numbers. The training using a variable viscosity showed excellent agreement between the predictions and the reference data over the whole range of Reynolds numbers. However, the absolute error grew non-linearly with the associated Reynolds number. This finding is in agreement with the results of Sun et al. [
19] who, however, covered lower Reynolds numbers. The results demonstrate the necessity for an appropriate scaling of the investigated problem. For a successful training, input and output dimensions of the neural net needed to match. This was also found by Laubscher and Rousseau [
14] who trained the non-dimensional Navier-Stokes equations.
For the investigation concerning the amount of labeled training data, we used different numbers
of training points:
. The total number of points
on boundaries and labeled data was consequently set to
.
Figure 5 exhibits the corresponding results. As illustrated, the prediction accuracy positively correlated with the amount of training data. This finding emphasizes the requirement of introducing enough labeled training data for successfully training the PINN and is in agreement with the observations of Pioch et al. [
27]. As
, the reduction of the error stagnated. In the subsequent analysis,
was, consequently, set to 800 for every trained AoA that was supplied with labeled training data.
In a third step, a comparison between three different labeled training data sampling strategies was performed. For this analysis, the prediction accuracy obtained using 800 randomly selected coordinates from a CFD grid as well as from a uniform grid with a spacing of 0.005 m was compared with a distribution of training data along several vertical and horizontal lines. The vertical lines were positioned at
and the horizontal lines were positioned at
. The number of points
on the lines was 767. Some illustrations concerning the training coordinates are given in the Appendix.
Figure 6 shows the results of this analysis. From the results it is evident that the randomly selected training data drawn from a CFD grid outperformed the other approaches.
To assess the effect of different turbulence models, we compared the standard
k-
model and the pseudo-Reynolds stress model, as described above.
Figure 7 illustrates the obtained predictions. The pseudo-Reynolds stress model exhibited a favorable agreement with the reference CFD data while the predictions obtained with the
k-
model were subject to significant deviations. This result demonstrates the importance of an appropriate turbulence modeling approach. A potential explanation for the different outcomes of the two models is given in
Section 5.6 and
Section 5.7.
As the described investigations revealed the superiority of a scaling using a variable viscosity, a random distribution of data points drawn from a CFD grid, the pseudo-Reynolds stress model, as well as the requirement to set at least 800 points of labeled training data per AoA into the domain, the training of the multidimensional surrogate model for variable angles of attack was set up accordingly.
We have further tested the application of a Fourier-feature NN as well as different weighting factors , , and . However, these investigations did not show significant effects on improving the performance and, hence, are not covered here.
6. Summary and conclusions
In this work, we trained a fully-connected feedforward PINN to predict the Reynolds-averaged flow field around the DU99W350 airfoil for arbitrary AoA between 10.0° and 17.5°. The PINN was trained with simulated data randomly drawn from a CFD grid for angles of attack of 10.0° and 15.0°. Additionally, satisfaction of the boundary conditions as well as the RANS equations were trained on the coordinates of a second CFD grid.
A sensitivity analysis was conducted investigating the Reynolds number, the amount and distribution of labeled training data, as well as the turbulence model. The error grew non-linearly with increasing Reynolds number and the PINN yielded accurate solutions over a range of Reynolds numbers only in case the problem was scaled varying the viscosity. It was shown that a minimum of labeled training data for supervised learning was required to obtain satisfying solutions at the Reynolds number of . Furthermore, the sampling of training data from the coordinates of a CFD grid exhibited superior results when compared with other non-adaptive sampling strategies. The k- model yielded unfavorable predictions while the equation-free pseudo-Reynolds stress model agreed favorably with the reference data.
The predictions of the PINN trained for variable AoAs exhibited good agreement with the reference data for interpolations between the training data as well as extrapolations outside of the labeled training data set. The predictions exhibited excellent agreement even within the separating boundary layer on the suction surface of the airfoil. However, the results also revealed challenges in predicting high gradient regions accurately. This was evident for the stagnation point at the leading edge of the airfoil.
The results of this work demonstrate the capability of feedforward PINNs to predict the Reynolds-averaged flow field around variable geometries at high Reynolds numbers.
Figure 1.
Features of the simulated reference flow field around the DU99W350 airfoil. 1: Stagnation point; 2: high velocity region on the suction surface; 3: boundary layer on the airfoil surface; 4: separated shear layer; 5: wake and recirculative vortices; 6: constant free flow; 7: high velocity region on the pressure surface.
Figure 1.
Features of the simulated reference flow field around the DU99W350 airfoil. 1: Stagnation point; 2: high velocity region on the suction surface; 3: boundary layer on the airfoil surface; 4: separated shear layer; 5: wake and recirculative vortices; 6: constant free flow; 7: high velocity region on the pressure surface.
Figure 2.
Architecture of the deployed neural network.
Figure 2.
Architecture of the deployed neural network.
Figure 3.
Illustration of the multidimensional hypercube training domain and the corresponding boundary conditions and training data.
Figure 3.
Illustration of the multidimensional hypercube training domain and the corresponding boundary conditions and training data.
Figure 4.
Results of the sensitivity analysis concerning the training accuracy for different Reynolds numbers. The relative deviations are presented on a logarithmic scale.
Figure 4.
Results of the sensitivity analysis concerning the training accuracy for different Reynolds numbers. The relative deviations are presented on a logarithmic scale.
Figure 5.
Results of the sensitivity analysis concerning the amount of labeled training data.
Figure 5.
Results of the sensitivity analysis concerning the amount of labeled training data.
Figure 6.
Results of the sensitivity analysis concerning the distribution method of labeled training data.
Figure 6.
Results of the sensitivity analysis concerning the distribution method of labeled training data.
Figure 7.
Results of the sensitivity analysis concerning different turbulence models. The absolute deviations are presented on a logarithmic scale.
Figure 7.
Results of the sensitivity analysis concerning different turbulence models. The absolute deviations are presented on a logarithmic scale.
Figure 8.
Results and absolute error of the axial velocity around the DU99W350 airfoil at an angle of attack of 15.0° using a PINN in comparison with the reference data. Note that for this angle of attack, reference data were provided for the training of the PINN.
Figure 8.
Results and absolute error of the axial velocity around the DU99W350 airfoil at an angle of attack of 15.0° using a PINN in comparison with the reference data. Note that for this angle of attack, reference data were provided for the training of the PINN.
Figure 9.
Results and absolute error of velocity and pressure around the DU99W350 airfoil at an angle of attack of 12.5° using a PINN in comparison with the reference data. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 9.
Results and absolute error of velocity and pressure around the DU99W350 airfoil at an angle of attack of 12.5° using a PINN in comparison with the reference data. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 10.
Results and residual of the axial velocity on the suction surface of the DU99W350 airfoil at an angle of attack of 12.5° using a PINN in comparison with the reference data. Left: axial velocity; right: residual between PINN and CFD results. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 10.
Results and residual of the axial velocity on the suction surface of the DU99W350 airfoil at an angle of attack of 12.5° using a PINN in comparison with the reference data. Left: axial velocity; right: residual between PINN and CFD results. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 11.
Results and absolute error of velocity and pressure around the DU99W350 airfoil at an angle of attack of 17.5° using a PINN in comparison with the reference data. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 11.
Results and absolute error of velocity and pressure around the DU99W350 airfoil at an angle of attack of 17.5° using a PINN in comparison with the reference data. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 12.
Results and residual of the axial velocity on the suction surface of the DU99W350 airfoil at an angle of attack of 17.5° using a PINN in comparison with the reference data. Left: axial velocity; right: residual between PINN and CFD results. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 12.
Results and residual of the axial velocity on the suction surface of the DU99W350 airfoil at an angle of attack of 17.5° using a PINN in comparison with the reference data. Left: axial velocity; right: residual between PINN and CFD results. Note that for this angle of attack, no reference data were provided for the training of the PINN.
Figure 13.
Lift and drag results for the DU99W350 airfoil as predicted by the PINN in comparison with the CFD results.
Figure 13.
Lift and drag results for the DU99W350 airfoil as predicted by the PINN in comparison with the CFD results.
Table 1.
PINN training parameters.
Table 1.
PINN training parameters.
Parameter |
Value/setting |
Architecture |
five hidden layer |
|
with 128 neurons each |
Optimizer |
Adam, L-BFGS-B |
Epochs |
30,000 (Adam) |
Learning rate |
(10,000 epochs) |
|
and (20,000 epochs) |
activation function |
tanh |
Number training points
|
3200 |
Number training points
|
800 at 10.0° and 15.0° each |
Number training points
|
6277 |
Loss terms weighting factors |
, ,
|