Next Article in Journal
Dynamic Process Analysis and Voltage Stabilization Control of Multi-Load Wireless Power Supply System
Next Article in Special Issue
Transition Prediction in Incompressible Boundary Layer with Finite-Amplitude Streaks
Previous Article in Journal
Application of System Dynamic Modelling for Evaluation of Carbon Mitigation Strategies in Cement Industries: A Comparative Overview of the Current State of the Art
Previous Article in Special Issue
Tracking Turbulent Coherent Structures by Means of Neural Networks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development and Validation of a Machine Learned Turbulence Model †

1
Department of Mechanical Engineering, Mississippi State University, Starkville, MS 39762, USA
2
Center for Advanced Vehicular Systems, Mississippi State University, Starkville, MS 39762, USA
3
DoD High Performance Computing Modernization Program PET/GDIT, Vicksburg, MS 39180, USA
4
Engineer Research and Development Center (ERDC), Vicksburg, MS 39180, USA
*
Author to whom correspondence should be addressed.
DOD DISTRIBUTION STATEMENT A. Approved for Public Release: Distribution Unlimited.
Submission received: 28 January 2021 / Revised: 23 February 2021 / Accepted: 25 February 2021 / Published: 8 March 2021
(This article belongs to the Special Issue Machine-Learning Methods for Complex Flows)

Abstract

:
A stand-alone machine learned turbulence model is developed and applied for the solution of steady and unsteady boundary layer equations, and issues and constraints associated with the model are investigated. The results demonstrate that an accurately trained machine learned model can provide grid convergent, smooth solutions, work in extrapolation mode, and converge to a correct solution from ill-posed flow conditions. The accuracy of the machine learned response surface depends on the choice of flow variables, and training approach to minimize the overlap in the datasets. For the former, grouping flow variables into a problem relevant parameter for input features is desirable. For the latter, incorporation of physics-based constraints during training is helpful. Data clustering is also identified to be a useful tool as it avoids skewness of the model towards a dominant flow feature.

1. Introduction

Engineering applications encounter complex flow regimes involving turbulent and transition (from laminar to turbulent) flows, which encompass a wide range of length and time scales that increase with the Reynolds number (Re). Direct Numerical Simulations (DNS) require grid resolutions small enough to resolve the entire range of turbulent scales and are beyond the current computational capability. Availability of high-fidelity DNS and experimental datasets is fueling the emergence of machine learning tools to improve accuracy, convergence and speed-up of turbulent flow predictions [1,2]. Machine learning tools depend on neural networks to identify the correlation between input and output features, and have been used in different ways for turbulent flow predictions, such as direct field estimation, estimation of turbulence modeling uncertainty, or advance turbulence modeling.
In the direct field estimation approach, the entire flow field is predicted using a ML approach, i.e., the flow field is the desired output feature. For example, Milano and Koumoutsakos [3] estimated mean flow profile in the turbulent buffer-layer region using Burger’s equation and channel flow DNS results as training datasets. Hocevar et al. [4] predicted the scalar concentration spectra in an airfoil wake using experimental datasets as training datasets. Jin et al. [5] estimated unsteady velocity distribution in the near wake (within four diameters) of a 2D circular cylinder using laminar solutions as training datasets and surface pressure distribution as input feature. Obiols-Sales et al. [6] developed Computational Fluid Dynamics (CFD) Network CFDNet—a physical simulation and deep learning coupled framework to speed-up the convergence of CFD simulations. For this purpose, a convolution neural network was trained to predict the primary variables of the flow. The neural network was used as intermediate step in the flow prediction, i.e., CFD simulations are solved as a warmup preprocessing step, then the neural network is used to infer the steady state solution, and following that CFD simulations are performed to correct the solution to satisfy the desired convergence constraints. The method was applied for range of canonical flows such as channel flow, flow over ellipse, airfoil, cylinder, where the results were encouraging.
For the turbulence model uncertainty assessment, the desired output feature is the error in a CFD solution due to turbulence modeling. For example, Edeling et al. [7] used velocity data from several boundary layer flow experiments with variable pressure gradients to evaluate the k-ε model coefficient ranges. Then, the variation in model coefficients were used to estimate the uncertainty in Reynolds averaged Navier Stokes (RANS) solution. Ling and Tempelton [8] compared several flow features predicted by k-ε RANS with DNS/LES results for canonical flows (e.g., duct flow, flow over a wavy wall, square cylinder etc.) and estimated errors in RANS predictions due to νT < 0, νT isotropy, and stress non-linearity.
Simulations on coarser grids require a model for the turbulent stresses (τ). The stress terms account for the effect of unresolved (or subgrid) turbulent flow on the mean (or resolved) flow for RANS (or LES) computations. The primary question for turbulence modeling is how turbulent stresses are correlated with flow parameters or variables? A review of the literature as summarized in Appendix A Table A1 and Table A2 shows that machine learning has been used to either augment physics-based models to improve their predictive capability [9,10,11,12,13,14,15,16,17,18,19,20,21,22] or develop standalone turbulence models [23,24,25,26,27,28,29,30]. The details of the input and output features and test and validation cases used in the studies are listed in the Tables, and the salient points of the studies are discussed below.
Parish and Duraisamy [9] analyzed DNS of plane channel flow to estimate the response surface of TKE production multiplier (β) as a function of four turbulent flow variables. The model was used to argument k-ω model. In a follow-on study [10] experimental data for wind turbine airfoils was used to adjust the Spalart–Allmaras (SA) RANS model νT production as a function of five turbulent flow parameters. The β function was reconstructed using an artificial neural network to minimize the difference between the experimental data and SA model results. The models were used for aposteriori tests, where it showed significant improvement over the standard RANS models and was found to be robust even for unseen geometries and flow conditions. He et al. [11] developed a similar approach, wherein adjoint equations were derived for SA model solution error due to νT production. The SA model predictions were compared with the experimental data at selected locations during runtime. Then, a solution of β was obtained using the adjoint equation to minimize discrepancy between the predictions and experimental data. The approach was applied for several canonical test cases and encouraging results were reported. Yang and Xiao [12] extended the work of Duraisamy et al. [9,10] to train a correction term for the transition model time-scale. Their model was trained using DNS datasets for flow over an airfoil at different angle of attack using both random forest and artificial neural network and using six different flow variables as input features. The trained correction term was implemented in a 4-Equation transition model, and applied for aposteriori tests involving unseen flow conditions (both interpolation and extrapolation mode) and geometry. The study reported good agreement for the transition location, validating the efficacy of such models.
Ling et al. [13] used six different DNS/LES canonical flow results to obtain coefficients of a non-linear stress formulation consisting of ten (10) terms involving non-linear combinations of the of rate-of-strain (S) and rotation tensors (Ω). The model coefficients were trained using a deep neural network as a function of five invariants of S and Ω tensors. The trained model coefficient map was applied for both apriori and aposteriori tests, including unseen geometry. The study reported that the ML model performs better than both the linear and non-linear physics-based models and performed reasonably well for unseen geometries and flow conditions.
Xiao and coworkers [14,15] trained a response surface of the errors in k-ε RANS turbulent stress predictions as a function of ten flow features using random forest regression. The model was validated for apriori tests for two sets of cases, one where the test flow and training flow were similar, and second for unseen geometry and flow conditions. It was reported that the model performed better for the former case. Wu et al. [15] investigated metrics to quantify the similarity between the training flows and the test flow that can be used to provide guidelines for selecting suitable training flows to improve the prediction of such models. Wang et al. [16] extended the above approach for compressible flows. A model was trained using DNS datasets with 47 flow features obtained using combination of S, Ω, ∇k, and ∇T, inspired by Ling et al. [13]. The model was validated for apriori tests for flat-plate boundary layer simulations. The study identified that the machine learned model predictions depend significantly on the closeness to the training dataset.
Wu et al. [17] developed a model to address the ill-conditioned solutions predicted by machine learned model during aposteriori tests (i.e., small errors in the modeled Reynolds stresses results in large errors in velocity predictions). For this purpose, a model was trained to account for the errors in k-ε RANS model stress predictions (both linear and non-linear terms) using seven turbulent flow features. The model was applied as an aposteriori test for flows involving slightly different geometry than the training case. Yin et al. [18] investigated the role of the feature selection and grid topology on the unsmoothness of the solution and large prediction errors reported in the above [17] study. They trained a model using 47 input features, inspired by Ling et al. [13]. The model was applied for apriori and aposteriori tests, wherein for the latter machine learned turbulent stresses were frozen. The study concluded that unsmoothness of the solution was primarily due to grid topology issues which results in discontinuities in the input features.
Yang et al. [19] used neural networks to train a wall-model for LES. The model was trained using three sets of input features: (1) wall parallel velocity (u||) and d; (2) u||, d and grid aspect ratio; and (3) u||, d, grid aspect ratio and ∇p. The model was coupled with a Lagrangian dynamic Smagorinsky model and applied for channel flow over a wide range of flow conditions Reτ = 103 to 1010. The study reported that inclusion of additional flow features such as grid aspect ratio and pressure gradient does not show significant improvement in the predictions.
Weatheritt and Sandberg [20] used symbolic regression to derive an analytic formulation for turbulent stress anisotropy. The model was trained using hybrid RANS/LES solutions and a regression map for the model coefficients were obtained as a function of rate-of-strain and rotation tensor magnitudes. The anisotropy formulation was used along with the k-ω model, and the model showed very encouraging results for both apriori and aposteriori tests including unseen geometries. Jian et al. [21] used a deep neural network to train a regression map of the model coefficients for an algebraic RANS model. The model was trained using a single parameter |S|k/ε as the input feature. The model was validated for both apriori and aposteriori tests, including extrapolation mode, i.e., Reτ larger than those in the training dataset. The study reported that the ML model performed better than the non-linear physics-based models due to its capability to better capture the large stress-strain misalignment and strong stress anisotropy in the near-wall region. Xie et al. [22] used neural networks to train model coefficients of a mixed subgrid stress/thermal flux model. The model trained compressible isotropic turbulence flow using six flow features, and it was reported that the machine learned model performs better than the physics-based models.
Comparatively limited efforts have been made to develop standalone machine learned turbulence models. Schmelzer et al. [23] used a symbolic regression approach to infer the model coefficients of an algebraic RANS model. The model was trained using DNS datasets using invariants of S and Ω tensors. The model was applied for unseen flow conditions, where the machine learned model performed better than the k-ω RANS. Fang et al. [24] developed a model for turbulent shear stress τ u v using DNS of channel flow. The model used a deep neural network trained using mean flow gradients and Reτ as key input parameters, and non-slip boundary condition and spatial non-locality were enforced during training. The model was applied for aposteriori tests involving unseen flow conditions, and it was reported that the model worked better than the model proposed by Ling et al. [16] due to the use of Reτ and boundary condition enforcement. Zhu et al. [25] developed a regression map of RANS turbulent eddy viscosity using a radial basis function network using SA RANS solutions for flow over NACA 0012 and RAE2822 airfoils at different angles of attack using eight input features based on flow, gradient and wall distance. The flow domain was separated into near-wall, wake, and far-field regions, and a different model was trained in each region. The study reported that using wall-distance as weights helped during training, but training using log-transformation of dataset did not help. The model was applied for aposteriori tests for the training geometry but for unseen angles of attack (flow condition), and good predictions were reported for the lift/drag coefficients, and skin friction distributions.
King et al. [26] developed a subgrid stress model using DNS of isotropic and sheared turbulence using velocity, pressure, S and grid parameters as input features. The model was applied for apriori tests for the isotropic and sheared turbulence test cases on coarse grid, where it performed better than the dynamic LES model. Gamahara and Hattori [27] used a feed-forward neural network to train a regression map of LES turbulent stresses. For this purpose, DNS results for plane channel flow for Reτ = 180 to 800 were filtered on up to four times (in each direction) coarser grids, and the filtered flow field was used as input features. The study used four different sets of input features involving S, Ω, wall-distance, and velocity gradients, and models were validated for apriori tests for unseen flow conditions. The study reported that the best model was predicted when u and d were used as input features. Further, it was reported that the machine learned models are more accurate than the similarity models, because of their ability to learn the non-linear functional relation between the resolved flow field and the subgrid stresses better than those prescribed in the physics-based model. Zhou et al. [28] used a similar approach to develop LES subgrid-scale model. The model was trained for isotropic decaying turbulence using gradients of filtered velocity and filter width (Δ). Yuan et al. [29] also developed an LES model using deconvolution neural network with isotropic decaying turbulence datasets for training. The model was trained using filtered velocity as the input parameter and validated for both apriori and aposteriori tests. Maulik et al. [30] used an artificial neural network to train a regression map of subgrid term for 2D turbulence using primary variables and its gradients as input features. The model was validated for both apriori and aposteriori tests. The ML model provided good predictions; however, it was reported that some measure of aposteriori error must be considered during optimal model selection for greater accuracy.
Some recent studies [31,32,33] have investigated the development of machine learning models for turbulent flow predictions by incorporating physics-based constrained during the training. These models provide an important direction to the applicability of the machine learning approach, but thus far they have used only for direct field estimation.
Overall, a review of the literature shows that
(1)
Machine learning has been primarily used for RANS model augmentation, where either turbulence production is adjusted, or nonlinear stress components are added to linear eddy viscosity term. Limited effort has been made to develop a stand-alone model, except for some recent effort focusing on modeling of subgrid stresses for LES.
(2)
Studies have used wide range input flow features for machine learned model training. There is a consensus that combining flow features into physically relevant flow feature is desirable, as this helps incorporate physics in machine-learning. Use of a large numbers of input features have been found to be helpful to some extent as it allows output features to be uniquely identified in the different flow regimes. However, they introduce additional sources of inaccuracy. For example, unsmooth solutions have been reported due to inaccurate calculation in the input features involving higher-order formulations of the derivative terms.
(3)
Machine learned models have been applied for both apriori and aposteriori tests for both unseen geometry and flow conditions, including Re extrapolation mode. The model in general perform well for unseen flow conditions, but issues have been reported for unseen geometries. In general, the machine learned models are most accurate when the test flow has similar complexity to the training flow.
(4)
Studies have reported issues during training due to overlap in the output features in different flow regimes. It has been tacked by using more input features, as discussed above, and by segregating the flow domain into regions with similar flow characteristics, such as near-wall, wake and far-field regions, and train separate models in each region.
In summary, there are open issues such as, “What is the best way to use machine learning for turbulence model development?” Should machine learning be used to augment an existing model or to develop a stand-alone model. The model augmentation approach builds on a baseline physics-based model; thus, it has some inherent robustness, especially when the model is used in extrapolation mode or for unseen geometry. However, this approach undercuts the advantage of machine learning. If it is expected that neural networks can accurately learn the errors in a RANS model and provide a universal model for the errors, then there is no reason to believe that same approach cannot provide a model for Reynolds stresses. Second, none of the studies in the literature have validated the machine learned model in a similar fashion as the physics-based model, i.e., how does it perform when started from ill-posed flow conditions, i.e., when simulation is started from an arbitrary initial condition, or does it provide grid convergence, or can independently query output in contiguous regions result in kinks in the solution. Apart from the above questions, there are additional issues regarding the best practices for optimizing machine learning approach itself [31].
The objective of this paper is to investigate the predictive capability of a stand-alone machine learned turbulence model and shed light on some of the above issues and constraints. To achieve these objectives, a DNS/LES database is curated for incompressible boundary layer solution available in the literature (i.e., channel flow and flat-plate boundary layer solution at zero pressure gradient), and a DNS database has been generated for oscillating channel flow. The datasets are used to train a ML response surface for the turbulent stress. The model is validated for apriori and aposteriori tests, and predictions are compared with DNS and RANS model results. The preliminary results for the channel flow case have been presented in ASME conference paper [32], and those of oscillating channel flow case has been presented in SC20 conference paper [33]. The results from the above publications have been further refined and are presented herein.
The novel aspects of this study include: development of a stand-alone RANS model, which has not received same level of attention as the RANS model augmentation; the effect of physics-based constraints during the training of a model is investigated, which has not been done before; the machine learned model is applied for an unsteady flow, thus far none of the studies have applied and tested machine learned model for unsteady flows; and the ability of the machine learned model to adapt to ill-posed initial flow conditions, as expected in a typical CFD simulation, and their ability to provide grid independent solution as expected for a RANS model are investigated.
The following section provides an overview of the DNS/LES datasets curated or procured in this research. Section 2 provides an overview of the machine learning approach. Section 3 and Section 4 focuses on development and validation of the machine learned model for boundary layer flows and oscillating channel flow case, respectively. Finally, some key conclusions and future work are discussed in Section 5.

2. Machine Learning Approach

The basic framework of the deep learning neural network is described in Figure 1, inspired from LeCun et al. [34]. The network consists of various layers, including the first input layer comprising of input features, the last output layer comprising of the required output features, and multiple hidden layers comprising of features units obtained using linear combination of feature units from the previous layer. Each layer, excluding the input layer, has input features (z) and output features (y). Let’s say for the pth layer with m feature units, the input and output features are z1 to zm and y1 to ym, respectively. The input feature units for the pth layer is obtained from the linear combination of the n output feature units in the layer above or (p − 1)th layer, as below:
[ z 1 ( p ) .. .. z m ( p ) ] = [ w 11 w 1 n w m 1 w m n ] w i j ( p ) [ y 1 ( p 1 ) .. .. y n ( p 1 ) ]
where superscript (p) is used to represent the pth layer, i varies from 1 to m and j varies from 1 to n, and w i j ( p ) are the unknown weights that need to be estimated. Thus, there are m neurons that connect the (p − 1)th and pth layer. Note that for the first hidden layer z vector are the input features, and for the last output layer y vector are the output features. Also note that some of the features in the hidden layers can be dropped depending on the threshold of weights permitted. Usually in deep learning applications, the input features in the input layer are scaled to vary in-between 0 and 1, and similarly the weights are positive and normalized such that w i j ( p ) = 1 . For each hidden and output layers, the output features are obtained from the input features on that layer as
y i ( p ) = f ( z i ( p ) )
where f ( .. ) is a pre-defined non-linear analytic activation function. The most common functions are
f ( z ) = { max ( 0 , z )   : R e L U z   f o r   z 0 ,   ( e z 1 )   f o r   z < 0   : E L U   e z e z e z + e z   :   Hyperbolic   tangent 1 1 + e z   :   Logistic
Rectilinear Linear Units (ReLU) function provides a linear dependency on the input features, exponential Linear Units (ELU) are same as ReLU for non-negative inputs, but use exponential modulation for negative inputs, hyperbolic tangent and logistic functions modulate both the positive and negative inputs. As pointed out by Lee et al. [35], ReLU is the most commonly used function for training deep neural networks because they do not suffer from vanishing gradients and they are fast to train. However, they ignore the negative values and hence lose information, which is referred to as “dying ReLU”. ELU on the other hand can capture information for the negative inputs, so they used ELU for training their physics-informed neural network.
The unknown weight matrix in each layer is estimated using backpropagation approach, as described below. The network is first initialized with constant zero weight for each layer and the error in the network prediction are obtained by comparing them with training dataset or true values (T) using a user specified analytic cost function, such as L2 norm
E = i = 1 , m ( y i T i ) 2
For the above defined analytic cost function, the derivatives of the errors in the output layer with respect to the output feature is
E y i = 2 ( y i T i )
and the derivative of the errors with respect to the input features is
E z i = E y i y i z i = f ( z i ) E y i
where, f is a known analytic function from Equation (3). Since the input features of a layer are linearly related to the output features of the previous layer as shown in Equation (1), the derivative of the errors with respect to the weights are computed as
E w i j ( p ) = E z i ( p ) z i ( p ) w i j ( p ) = y j ( p 1 ) E z i ( p )
Lastly, the weights in each layer are adjusted as
w i j ( p ) | i t + 1 = w i j ( p ) | i t α E w i j ( p )
where α is learning rate, and subscript it represents the iteration level. Commonly available ML softwares provide optimizers which dictate the learning rate. For example, adaptive moment estimation (ADAM) optimizer [36] available in Keras application programming interface (API) requires a user specified initial learning rate, but the rate is adaptively adjusted during training. Note that in deep neural networks “iterations” refer to running a subset of the data (batch) forward and backwards through the network, whereas “epoch” refers to running all the training data forward and backward through the network. Thus, epochs are not same as iterations, unless the entire dataset is the batch size. Further note that the terms on RHS is known from Equations (5) and (6) for the output layer. For the hidden layers a backpropagation approach is used, where the derivative information is computed based on information from the layer ahead starting from the output layer, Equation (5), as below:
E y j ( p 1 ) = w i j ( p ) E z i ( p )
and then Equations (6) and (7) are used to obtain E w i j ( p 1 ) .
In the physics-informed machine learning (PIML) approach [37] the cost function is modified to include the residual in the governing equations ( ); thus, Equation (4) is modified to
E = i = 1 , m ( y i T i ) 2 + i = 1 , m ( i ) 2
Note that Equation (5) remains changed.

3. Test Cases and Database for Model Training

Two different text cases have been considered in this study, steady and unsteady channel flow. For the former, the mean flow solution describes the inner layer of a flat-plate boundary layer (with zero pressure gradient) at a fixed location on the plate. This case is a fundament test case for turbulence model development and several DNS studies are available for a range of flow conditions. The flow pattern for this case reduces to a one-dimensional (1D) problem. For the second test case, the inner boundary layer undergoes unsteadiness due to prescribed pressure gradient pulse. The mean flow pattern for this case reduces to an unsteady 1D problem, which adds another level of complexity to the first test case.

3.1. Plane Channel Flow

3.1.1. Governing Equation

This test case focuses on the simulation of the inner layer of the turbulent boundary layer under zero-pressure gradient at a fixed location on a flat-plate. The governing equations for such flow condition can be derived from the incompressible Navier–Stokes equations under the assumption that the flow is 2D and steady, and streamwise gradients are negligible compared to those along the wall normal direction (refer to Appendix B for derivation) as below:
u t = F ( τ w ) + ν 2 u y 2 + y ( τ u v )
where y is the coordinate direction normal to the wall and u is the ensembled averaged streamwise velocity. Note that the correlation between streamwise and wall-normal turbulent velocity fluctuations ( τ u v = u v ¯ , which are the shear stresses) is the only unknown quantity that needs to be modeled. Also note that although the above equations are valid for flows with zero pressure gradient, the above equation includes a body force term ( F ), which can be misconstrued as pressure-gradient term. This term is added because the simulations are performed for flow between two flat-plates; thus, a body force is required to balance the momentum loss due to wall friction and achieve a steady state. The body force term F ( τ w ) is a function of wall shear stress τ w expected at the simulated flat-plate location. The wall shear stress can be expressed in terms of friction velocity u τ , and is a fixed input parameter:
τ w = u τ 2
and
F ( τ w ) = u τ 2 / H
where H is the half channel height. Note that the simulated flow conditions in a channel flow case can be changed simply by changing the wall friction value (or the applied body force term). Further note that a time-derivative term in the governing equation is a pseudo time derivative, and used as a residual term and provides a measure of the solution convergence. A steady state solution is achieved as the time derivative term approaches zero.
Commonly used linear URANS models express the turbulent shear stress as [38]
τ u v = ν T u y
where, ν T is an unknown turbulent eddy viscosity. The one-equation URANS model requires solution of an additional transport equation for turbulent kinetic energy (k)
k t = ν T ( u y ) 2 Production 0.1643 k 3 / 2 Dissipation + y { ( ν + ν T ) k y } Diffusion
and turbulent eddy viscosity is obtained as
ν T = 0.3 k
The turbulent length scale is prescribed as
= κ d ( 1 e y + / 26 )
where von-Karman constant κ = 0.41 and d is distance from the wall. Refer to Warsi [38] for details for the modeling. Similar to the streamwise velocity equation (Equation (11)), the time derivative can be perceived as a residual term for steady simulations, and a steady state is achieved as term approaches zero. The time derivate term provides the time-accurate solution for URANS simulations.

3.1.2. DNS/LES Database

A DNS/LES database is curated to train a ML response surface for τ u v . As summarized in Table A3, the database includes, 21 channel flow DNS cases with Reynolds numbers ranging from Reτ = 109 to 5200 [39,40,41,42,43,44,45]; 18 flat-plate boundary layer with zero-pressure gradient DNS cases with Reynolds numbers ranging from Reθ = 670 to 6500 [46,47,48]; and 14 flat-plate boundary layer with zero-pressure gradient LES cases with Reθ = 670 to 11,000 [49,50]. The database contained around 20,000 data points of which 3.4% were in the sub-layer, 5.7% in buffer layer and 90.9% in the log-layer or the outer layer. Note that all the datasets are for flat-plate boundary layer with zero-pressure gradient, where the channel flow cases represent only the inner-layer.

3.1.3. ML Model Training and Refinement Using Apriori Tests

ML model is developed using the cost function in Equation (4), which is referred to as data-driven machine learned (DDML) model, and cost function including the residual in the governing equations, i.e., Equation (10), which is referred to as physics-informed machine learned (PIML) model. For the PIML model training, the residual of the governing equation is obtained using the integral form of the governing equation (Equation (11)) at steady state as below:
i = ν u i y + τ u v i u τ 2 / H d
where subscript ‘i’ denotes the solution at the epoch level.
The DDML model is trained used a multilayer perceptron (MLP) neural network consisting of two dense hidden layers, each with 512 neurons with ReLU activations, and a linear activation on the output layer. ADAM optimizer is used during training with a mean absolute error loss function, generating 276,000 trainable parameters. To mitigate over-fitting during training, each layer has a 20% dropout [51]. The PIML model is trained using an eight-layer deep neural network with each hidden layer having 20 neurons and a hyperbolic tangent activation function. The model is optimized using an L-BFGS quasi-Newton full-batch gradient-based optimization method and iterated for 2400 epochs. For the model training, the dataset is not separated into training and test sets; but rather 70% of the dataset is chosen randomly. Thus, it is possible that there could be an overlap in the datasets used for training and apriori tests. A parametric study was also performed by choosing 80% and 100% of datasets, which did not show much effect on the model predictions in apriori tests. Further note that the number of layers and neurons used for training DDML and PIML are different. For DDML training, a parametric study was performed using more layers and different number of neurons. The tests revealed that increasing the number of layers increases the training time, but does not necessarily improve the training accuracy. Increasing the number of neurons in lieu of layers was found to be computationally efficient and helped in reducing the training error. The numbers of layers and neurons used in the study are found to provide optimal model in terms of computational cost and accuracy. The optimal width and depth of neural network architecture also depends on the amount of data available for training the networks [52]. A shallow depth network is partially due to comparable small training dataset. For the PIML training a similar test was not performed and the training set-up was same as that used by [37]. Overall, although the numbers of layers and neurons differed between the PIML and DDML model training, but the neural network architecture for both provided the needed accuracy.
The machine learned turbulence model seeks to obtain a response surface of the shear stress τ u v , which is the output feature. The input features are the flow parameters. Referring to Figure 1, if a training batch uses N number of data points and each data point has M flow parameters, then the total number of input features x i : i = N × M , and total number of output features y m : m = N . The batch size was 128 for this case. Note that even though the M flow parameters at each data point are related; however, during the training they are considered independent of each other. The PIML model may leverage the correlation between the input features at a datapoint, as the feature set is expected to satisfy the governing equations.
Considering all the possible flow variables; the response surface of the shear stress τ u v is expected to have following functional form:
τ u v = f ( u , U c ,   d ,   ν , u τ , u y )
Instead of providing the flow features independently, they are grouped or non-dimensionalized into physically relevant parameter. The velocity derivative is a key term, which captures the rate-of-strain in the flow. It is commonly accepted that turbulent stresses can be modeled by extending Stokes’ law of friction (for viscous stresses), except that the linear assumption is debatable [38]. Thus, the velocity derivative normalized by the centerline velocity U c and channel height H is used as an input parameter. Figure 2 provides an overview of the correlation of τ u v with respect to key input features u y . As evident, the dataset shows quite an overlap especially for 5 u y 30 . Thus, it is clear that this parameter alone is not sufficient to train a reliable regression map.
The flow viscosity ν is one of the fundamental bulk fluid property and dictates one of the key non-dimensional parameters for the boundary layer flow, which is the Reynolds number (Re). Re is usually defined based on global properties, such as channel height H and centerline velocity U c , Re = U c H/ν. This bulk flow parameter has little significance to the local stresses; thus, is not a suitable input parameter for machine learning. Rather, Reynolds number based on wall distance R e d = U c d ν , or based on local velocity and wall distance R e l = u d ν seem to be better choice, as they also incorporate local parameters. Note that for channel flow datasets, both the global parameters ν and u τ are related via near-wall velocity gradient, and do not provide any additional information about the flow. Thus, when training a model just using the channel flow datasets u τ is not required. However, when the datasets include both flat-plate boundary layer and channel flow datasets, then u τ provides additional identifying information (as discussed below).
A preliminary apriori analysis is performed using just the channel flow DNS datasets to evaluate which Re formulation ( R e d or R e l ) is better, and also to investigated how inclusion of higher-order terms of rate-of-strain affects the model. For this purpose, the DDML model was trained using the following four sets of input parameters:
τ u v = { f { R e d = U c d ν , u y } f { R e l = u d ν , u y } f { R e l = u d ν , u y , ( u y ) 2 } f { R e l = u d ν , u y } , ( u y ) 2 weighted
As shown in Figure 3a, the stresses obtained using input feature set ( R e d , u y ) is around 8% under predictive whereas those using obtained using ( R e l , u y ) compare very well and is only 2% over predictive. An accurate prediction by the latter is not surprising, as local Reynolds number R e l = u + y + is a physical quantity that separates the flow regime, i.e., R e l < 25 is sub-layer; 25 R e l < 418 is buffer-layer; and R e l 418 is log-layer. The advantage of R e l over R e d is also evident in Figure 3b,c, where the former shows a better collapse in dataset in the buffer layer. Further, including ( u y ) 2 as an input feature deteriorates the model performance and the stresses are overpredicted by 9%. However, using ( u y ) 2 as a weighting parameter improves the results and the L2 norm error for this case is estimated to be 2.9 × 10−4, slightly better than the L2 norm error of 5.2 × 10−4 obtained using the feature set ( R e l , u y ) .
When the entire channel and flat-plate database is considered, we need an additional parameter to demarcate between the inner and outer boundary layer. The inner and outer layer demarcation can be achieved by using a combination of y + = d u τ ν and u y , where large y + and u y 0 corresponds to the start of the outer layer. Thus, a model was trained using the following set of input parameters:
τ u v =   f { R e l , y + , u y }
In addition, different weighting functions, ML-1 through ML-4 as shown below, were used during training of the DDML response surface as the data show significant overlap (as pointed out in Figure 2) at the intersection of the inner and outer layers:
  • ML-1: No weighting
  • ML-2: Weighted using curvature of the profiles
  • ML-3: τ u v levels were expanded to separate out the curves
  • ML-4: Not-weighted for τ u v   10 3 but weighted for τ u v >   10 3
Note that ( u y ) 2 weighting was also considered. But, similar to the earlier test this weighting did not show significant improvement over ML-1. This is expected as neural networks should be able to learn the dependency on higher-order forms of the input parameters all by itself. Thus, such weighting is not considered further.
An apriori analysis for a range of channel flow conditions in Figure 4 shows that the DDML response surface obtained without any weighting (ML-1) is consistently under predictive. Those obtained using profile curvature as a weighting function (ML-2) results in a waviness in the profile and is not recommended. DDML using τ u v weighting improves the results, and ML-4 provides somewhat better results than ML-3. The model shows some limitations for the prediction of peak stress for Reτ = 1000, which needs to be further investigated. The PIML model shows the best predictions among all the models for the entire range of Reτ, and compares very well with DNS, including those for Reτ = 1000. In addition, PIML removes the trial-and-error approach to train the model.

3.2. Oscillating Plane Channel Flow

Oscillating channel flow is a canonical test case used in the literature to understand turbulence flow physics and validate the predictive capability of LES models. Scotti and Piomelli [53] performed DNS and LES of oscillating channel flow to study the effect of pressure gradients on the modulation of the viscous sublayer, turbulent stresses and the topology of the coherent structures. In this test case, an unsteady pressure pulse (body force term) is applied along the streamwise direction, which introduces periodic unsteadiness in the flow.

3.2.1. Governing Equation

The governing equation for the mean flow for this case can be obtained similar to the inner boundary layer equations as discussed in Appendix B as below:
u t = F ( t ) + ν 2 u y 2 + τ u v y
where
F ( t ) = d P 0 d x [ 1 + α 0 cos ( ω t ) ]
is the prescribed pressure pulse. The above equation is derived under the assumption that mean flow is 2D in nature and wall-normal gradient is more dominant than the streamwise gradients. Both of these assumptions are valid for this case. Because DNS is performed using periodic domain, the mean flow variation along the streamwise direction is assumed negligible, and the velocity changes occur in time. Thus, this test case represents how boundary layer flow, at a fixed location on a flat plate, varies in time as the free-stream pressure gradient in varied. Further, note that the d P 0 d x term corresponds to the F ( τ w ) forcing required to obtain steady state solution for flow between two flat-plates, which is the baseline flow. The term d P 0 d x α 0 cos ( ω t ) accounts for the variation of the free-stream pressure gradients. Similar to the steady inner-boundary layer case, the closure of the above equation requires modeling of shear stress τ u v . The physics-based one-equation URANS model described in Section 3.1.1 is valid for this case as well.

3.2.2. DNS Database

The DNS for this case is performed using an in-house parallel pseudo-spectral solver, which discretizes the incompressible Navier–Stokes equations using Fast Fourier Transform (FFT) along the homogenous streamwise and spanwise directions and Chebyshev polynomials in the wall-normal direction. The solver is parallelized using a hybrid OpenMP/MPI approach, scales up to 16K processors on up to 1 billion grid points, and has been extensively validated for LES and DNS studies [54,55].
DNS were performed for three different (high, medium, and low) streamwise pressure pulses, as used in [53]. The simulations were performed using as domain size of 3π × 2 × π along the streamwise, wall normal and spanwise directions, respectively, on a grid consisting of 192 × 129 × 192 points. The domain size is consistent with those used in [52], but the grids are finer. The details of the simulation set-up are provided in Table A4. The simulations were performed using periodic boundary condition along streamwise and spanwise directions, and no-slip wall at y = ±1.
The DNS results were validated in a previous study [56] for the prediction of the alternating (AC) and mean (DC) components of the mean flow against results presented in [53]. The AC components were obtained by decomposing normalized-planar-averaged velocity profile at every one-eight cycle using Fast Fourier Transform at each wall normal location. The DNS results were found to be in good agreement with the available benchmark results. The low frequency pulse resulted in re-laminarization and transition behavior, which is a challenging case for RANS predictions. The high frequency case was primarily driven by the pressure gradient, and turbulence levels were small. The medium frequency case provides a compromise between the two extremes and is used in this study.
The variation of the wall shear stress and streamwise and wall-normal velocities for the medium frequency case is shown in Figure 5. The positive pressure gradient generates adverse flow conditions and decelerates the flow (and decreases wall shear stress magnitude) and the negative pressure gradient generates favorable flow conditions which accelerates the flow (and increases wall shear stress). The peak wall shear (and velocity) is observed around 3/4th cycle (t = 0.75). The wall normal velocity shows that the ejection events are subdued during the first part of the oscillation cycle (descending leg), i.e., at t = 0.25 the ejection events are limited within the quarter channel and diminished at the mid-cycle (t = 0.5). The ejection events are enhanced during the ascending leg of the pressure pulse. The ejection events are very prominent for t = 0. Thus, shear stress is expected to be generated during the last part of the oscillation cycle.

3.2.3. ML Model Training Using Apriori Tests

For training of the ML turbulence model, the 3D DNS dataset was processed to obtain an unsteady 1D dataset. Due to the use of a periodic boundary condition along the streamwise direction, the simulation domain is expected to move in time [57], and as demonstrated in Figure 6a. Further, the results at any time step represents multiple (every grid point in the streamwise-spanwise plane or 192 × 192) realizations of the turbulent flow field expected at that oscillation phase. The solutions were averaged over these realizations to obtain mean 1D flow along the wall normal (y) direction. The 1D solutions every 100th time step or 1/400th pressure oscillation cycle were used to generate a 2D y-time map of the mean flow as shown in Figure 6b. Solutions were collected over three pressure oscillation cycles resulting in 129 × 1201 (or around 155,000) data points.
For this case, k is also added as an unknown turbulence quantity, and the machine learned model seeks to obtain regression map with two outputs ( τ u v ,   k ). Based on the experience of the boundary layer case, the input flow parameters for this case are identified to be the time varying and steady mean flow quantities as below:
( τ u v ,   k ) O u t p u t = f ( R e l = u d ν ; F ( t ) ; y + ( t ) = d u τ ( t ) ν ; y 0 + = d u τ 0 ν ; U y ) I n p u t s
u τ ( t ) = | τ w ( t ) | ρ
where u τ ( t ) is the local friction velocity and u τ 0 is the baseline wall friction corresponding to d P 0 d x .
Only the DDML model is trained for this case, as temporal derivatives could not be computed during training. The model was trained using a three-layer-deep neural network with 512 neurons in each hidden layer with ReLU activation functions. The final layer was a linear fully connected layer. The model was trained for 200 epochs using an ADAM optimizer with a batch size of 128 and a learning rate of 0.01. During the training, the L2 norm error dropped by an order of magnitude in the first 25 epochs, and the error drop stalled thereafter.
The oscillating channel flow involves a wide range of turbulence regimes unlike the boundary layer case which only has sub-, buffer- and log-layer regimes. As shown in Figure 7 in this case, the turbulence is most prominent toward the end of the pressure-oscillation cycle, and varies significantly along the wall-normal direction. Thus, training a ML model using the entire dataset may be skewed toward the more prominent features. For example, for the boundary layer case the models perform much better in log-layer than in the buffer-layer, as 91% of the datasets were in log-layer region. The Birch clustering algorithm sklearn was used to cluster the datapoints into unique flow regimes. The clusters were automatically determined using a non-dimensional threshold of 0.3, which resulted in 412 clusters with (min, max, mean, std) of points to be (1, 10,700, 376, 937), as shown in Figure 7. Note that the clustering preserves the periodicity of data.
Several ML models were trained by randomly sampling various percentage of points within each cluster from 0.25% to 100% of the total dataset. As shown in Figure 8b–d, training using just one datapoint per cluster (0.25% of total) results in errors of 35%. The error level reduced to 28% when four points are used per cluster (or 1% of total). The error levels were around 12% when 10% of the datasets (either 10 points or 10% of each cluster) was used. A similar error levels were obtained when model was trained using all the datapoints. In addition, training using 10% of the data points required around eight times smaller computational cost compared to those using all of the datapoints. This indicates that an intelligently sampled dataset can generate reliable machine learned models at a lower computational cost.

4. Aposteriori Tests of the ML Model

For aposteriori tests, the machine learned turbulent shear stress (and TKE) regression map obtained in the pre-processing step above was coupled with the parallel pseudo-spectral solver used for DNS of oscillating channel flow (discussed in Section 3.2.2). Figure 9 provides a schematic diagram demonstrating the coupling of machine learned turbulence model with the CFD solver. As shown, the stress (or TKE) regression map is queried every time step using the flow predictions at that iteration level, and the stresses are used to advance the solution to the next time step. Both the test cases considered in this study are 1D in nature, i.e., mean streamwise velocity varies only along the wall-normal direction; thus, the simulations were performed using only two points in the streamwise and spanwise direction.

4.1. Steady Plane Channel Flow

The ML trained turbulent shear stress response surface (generated in previous section) was used in the solution of the boundary layer equation (Equation (11)). The equations were solved using a domain y = [−1,1] with no-slip boundary condition (u = 0) at y = ±1. The simulations were performed using 49 (coarse), 65 (medium), and 97 (fine) grid points. Two sets of simulations were performed for this case. For the first set, the simulations were started from a converged channel flow RANS solution obtained using one-equation model, and the second set of simulations were started from an ill-conditioned velocity field. The DDML and PIML model simulations were performed using input feature set ( R e l , y + , u y ) (as shown in Equation (21)). Further for the DDML simulations ML-3 and ML-4 weighting which provided the best predictions for apriori tests were considered.
Figure 10a shows solution convergence on all the three grids obtained using DDML when solution is started from RANS solution. The convergence time history shows a kink in residual early-on in the simulation, but eventually solution shows steady convergence. The convergence is faster for the fine grid and slowest for the coarse grid. On the coarse grid the residual converges to around 2 × 106. Whereas for both medium and fine grid simulations, the residual keep dropping even below 1.2 × 10−6 after 30 thousand (K) pseudo timestep iterations. The PIML results in Figure 10b shows a similar convergence.
Figure 11 compares the results obtained using DDML and PIML (after 30K pseudo timesteps) with DNS and RANS results. Among the DDML models, the model obtained using ML-4 weighting performed better than those obtained using ML-3. This is in contrast to the apriori tests, where ML-3 performed better than ML-4. Note that the one-equation model underpredicts the velocity at the center compared to the DNS, and both the DDML and PIML models resolve the issue. The DDML model predictions improve with grid refinement, whereas the largest improvement is obtained in between coarse and medium grids. However, the results show oscillation near the peak, which is clearly evident in the turbulent and shear stress predictions and more prominent for the fine grid predictions. The oscillations suggest that the “query output” jumps between the database curves. For the boundary layer equations, one of the physical constraints is the shear stress must satisfy C1 continuity, i.e., both stress and its derivative are continuous in space (i.e., along y direction). Node/point-based queries are definitely not satisfying this constraint, which is probably the cause of the oscillations.
To resolve the stress oscillation issue in the DDML model predictions, solution smoothness was enforced by implementing region-based query, i.e., use averaged τ u v output for the input parameter sets at the node and its two neighboring nodes, i.e.,
τ u v | j = 1 8 [ 2 τ u v ( R e l , j , u y | j , y j + ) + τ u v ( R e l , j 1 , u y | j , y j + ) + τ u v ( R e l , j + 1 , u y | j , y j + ) + τ u v ( R e l , j , u y | j 1 , y j + ) + τ u v ( R e l , j , u y | j + 1 , y j + ) + τ u v ( R e l , j 1 , u y | j 1 , y j + ) + τ u v ( R e l , j + 1 ; u y | j + 1 ; y j + ) ]
As shown in Figure 11, the above averaging approach helps get rid of oscillations, and the results improve with grid refinement consistent with grid convergence.
The PIML model predictions show only marginal improvement with grid refinement, and the results do not show oscillations similar to the DDML model. Further, when averaging was applied it significantly deteriorated the performance of the model. In general, the PIML model performs better than the DDML especially in the lower log-layer region, i.e., y+ ~ 30 to 40.
In the next test set, the simulations were started from ill-conditioned initial condition profile:
u = U c ( 1 d 8 )
This profile is much steeper than that of the turbulent boundary layer, as shown in Figure 12 (It = 0 plot), and the associated set of input features are not expected to be present in the DNS/LES database. Thus, the model operates in extrapolation mode. As shown in Figure 12a, the DDML model predictions converge to the correct sub-layer and log-layer profiles, but shows significant differences in the buffer- and lower log-layer (20 ≤ y+ < 80). A good prediction of the sub- and log-layer is very encouraging, and suggests that the extrapolation issue can be addressed by incorporating physics constraints during query process, such that the query recognizes that the inputs are out of bound of the database and provides a better educated guess of the output. As shown in Figure 12b, the PIML model converges slowly (as shown in the convergence plot in Figure 10b) to the DNS profile. This suggest that a well-trained machine learned model can converge to the correct solution.

4.2. Oscillating Plane Channel Flow

Three sets of simulations were performed for this case using 65 grid points in the wall-normal direction. The simulations were performed using the DDML model trained using input feature set shown in Equation (24), and using 10% of the datasets (i.e., either 10 points or 10% of each cluster). For set #1, the DDML model was used in apriori mode, i.e., simulations were performed using one-equation URANS model, and the local flow predictions were used to query the regression map. For set #2, the DDML model was used in aposteriori mode, where the simulations were started from fully converged channel flow velocity profile corresponding to Reτ = 350 obtained using RANS. For set #3, both URANS and DDML simulations were started from an ill-posed initial flow condition i.e., u t = 0 = U c ( 1 d 8 ) .
The DDML model predictions from set #1 and #2 are compared with DNS and URANS predictions in Figure 13. The URANS model performs quite well for the mean flow, but predicts significantly diffused shear stress and the peak values are 9% under predicted. The largest error is obtained for the turbulent kinetic energy for which the peak values are underpredicted by 60%. The DDML model predictions (in apriori mode) show significantly better shear stress predictions than the URANS model. The predictions have errors on the order of 3–5% for the mean velocity and shear-stress, and peak TKE are overpredicted by 15%. Since the DDML model uses the velocities and derivatives predicted by the URANS model, the improved prediction by the former can be attributed to its ability to learn the non-linear correlation between the turbulent stresses and rate-of-strain. The DDML model also works very well in the aposteriori mode, and both the mean velocity, shear stress and k compare within 8% of the DNS. Note that the simulations used a (wall-normal) grid and time step size two times smaller and 10 times larger, respectively, compared to those of DNS, thus the predictions are considered reasonably accurate.
As shown in Figure 14, for the simulations started from ill-posed initial flow conditions, the URANS solutions show large differences with the DNS during the early part of the simulations, but the solution slowly recovers, and the solutions for the second and third cycles are very similar to those for the earlier case. This is because the flow is primarily driven by the pressure gradient, and the flow adapts to the pressure variations. The DDML model adjusts to ill-posed initial flow condition much faster than the URANS model, and results compare very well with the DNS, and the predictions errors are similar to those of set #2.

5. Conclusions and Future Work

This study investigates the ability of neural networks to train a stand-alone turbulence model and the effects of input parameter selection and training approach on the accuracy of the machine learned regression map. To achieve these objectives, a DNS/LES database is curated and/or developed for steady and unsteady boundary layer flows, for which the mean flow simplifies to a 1D steady and unsteady problem, respectively, and closure of the governing equations require modeling of turbulent shear stress. The database was used to train data driven and physics-informed machine learned turbulence model. For the latter, the residual in the governing equation solution was incorporated in the cost function during the model training. The model was validated for apriori and aposteriori tests, including ill-posed flow condition.
Overall, the results demonstrate that machine learning can help develop a stand-alone turbulence model. Moreover, an accurately trained model can provide grid convergent, smooth solutions, which works well in extrapolation mode, and converge to a correct solution from ill-posed flow conditions. The accuracy of the machine learned response surface depends on
  • The choice of input parameters. Feature engineering was used to find the optimal input features for the neural network training. It was identified that grouping flow variables into a problem relevant parameter improves the accuracy of the model. For example, a model trained using Re based on local flow velocity and wall distance is more accurate compared to the model trained using Re based on global flow. Furthermore, higher order functions of an input variable, such as square of the rate-of-strain along with rate-of-strain, does not help in improving the accuracy of the map. However, they may be used as weighting function to reduce the overlap in the datasets; and
  • How the database is weighted to minimize the overlap between the datasets. This requires a trial-and-error method to come up with an appropriate weighting function. A better way to improve the accuracy of the regression surface is to include physical constraints to the loss function during training, which is referred to as the PIML approach. However, it is not straightforward to incorporate physical constrains during the training due to issues in calculation of the derivates, such as temporal derivatives, for unsteady problem. Data clustering is also identified to be a useful tool to improve accuracy of the machine learned model and reduce computational cost, as it avoids skewness of the model towards a dominant flow feature.
Herein, machine learning was applied for cases which are very similar to the training datasets, which limits the applicability of the model, as well as does not sufficiently challenge the robustness of the machine learning approach. The ongoing work is focusing on generation of a larger database encompassing steady and unsteady boundary layer flows, including separated flow regimes. A model trained using such a database will help in development of a more generic turbulence model for boundary layer flows. It is expected that data clustering will be very helpful for training such as model due to the presence of wide range of turbulence characteristics. On a final note, the machine learned models were found to be extremely computationally expensive compared to the physics-based model (the former was around two order of magnitude more expensive). This was because of the added cost associated with ML query every iteration for every grid point. This research primarily focused on the accuracy of the model and efforts were not made to improve the efficiency of the ML model query. However, practical application of ML model would require investigation of approaches to improve the computational efficiency of run-time ML query.

Author Contributions

Conceptualization, S.B.; methodology, S.B., G.W.B., W.B. and I.D.D.; software, S.B., G.W.B. and W.B.; validation, S.B. and G.W.B.; formal analysis, S.B.; investigation, S.B. and G.W.B.; resources, I.D.D.; data curation, S.B.; writing—original draft preparation, S.B.; writing—review and editing, G.W.B., W.B. and I.D.D.; visualization, S.B. and G.W.B.; supervision, I.D.D.; project administration, I.D.D.; funding acquisition, S.B., G.W.B. and W.B. All authors have read and agreed to the published version of the manuscript.

Funding

Effort at Mississippi State University was sponsored by the Engineering Research & Development Center under Cooperative Agreement number W912HZ-17-2-0014. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Engineering Research & Development Center or the US Government. This material is also based upon work supported by, or in part by, the Department of Defense (DoD) High Performance Computing Modernization Program (HPCMP) under User Productivity Enhancement, Technology Transfer, and Training (PET) contract #47QFSA18K0111, TO# ID04180146.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Symbols
QSecond invariant of rate-of-strain tensor
uLocal flow velocity vector
uTurbulent velocity fluctuation vector
U c Free-stream (global) velocity (or centerline channel velocity)
HHalf channel height
uτFriction velocity
pPressure
kTurbulent kinetic energy
εDissipation
νKinematic molecular viscosity
νTTurbulent eddy viscosity
dDistance from the wall
ωSpecific dissipation, ε/k
y+Wall distance normalized by friction velocity, d u τ / ν
τwWall shear stress
Gradient operator
SRate-of-strain tensor
ΩRotation tensor
|S|Magnitude of rate-of-strain tensor, 2 S : S
|Ω|Magnitude of rotation tensor, 2 Ω : Ω
PProduction of turbulent kinetic energy
ReReynolds number based on global flow variables, such as U c and geometry length
RetTurbulent Re based on distance from the wall, k d / ν
RedRe based on distance from wall, d U c
RelReynolds number based on the local velocity and distance from the wall, u d / ν
τ Shear stress tensor
τuvTurbulent shear stress component in x-y plane, u v ¯
MLMachine Learning
· Dot product
: Double dot product
Terminology
DatabaseCurated DNS/LES datasets for ML training
Response surfaceOutput from ML training
Input featuresInput flow variables for ML
Output featuresFlow variable for which the response surface is generated
Query inputsInput variables to query the response surface
Query outputOutput variables obtained from the query of the response surface
Unseen case (flow)Geometry (or flow condition) not used during ML training
Apriori testML model is applied as a post-processing step
Aposteriori testML model is coupled with CFD solver and its prediction is used during runtime

Appendix A

Table A1. Literature review of machine learning approach for turbulence model augmentation.
Table A1. Literature review of machine learning approach for turbulence model augmentation.
ReferenceResponse SurfaceTurbulence ModelInput FeaturesTraining FlowsValidation CaseComments
Parish and Duraisamy [9]TKE production multiplier: β (ηi) (Neural network)k-ω RANS:
ν T ( u ¯ y ) 2 β ( y ) α k ω + y [ ( ν + σ k ν T ) k y ]
4 features: ηi = | S | k ε , k d ν ,   P , y + DNS: Plane channel flow, Reτ = 180, 550, 950, 4200Plane channel flow,
Reτ = 2000
  • Inversion process as an optimization problem to minimize the difference between the experimental data and RANS results.
  • Correction term was applied for aposteriori tests.
  • Robust for unseen geometries and flow conditions
Singh et al. [10]Turbulence production multiplier: β (ηi). (Neural network)Spalart-Allmaras RANS:
D ν ˜ D t = β ( x ) P ( ν ˜ , U ) P r o d u c t i o n
D ( ν ˜ , U ) D e s t r u c t i o n + T ( ν ˜ , U ) T r a n s p o r t
5 features: ηi = | Ω | ,   χ = ν T ν ,   | S | | Ω | ,   τ τ w ,   P ε Experiment: lift coefficient (CL) and surface pressure (CP) for wind turbine airfoils S805, S809, S814, Re = 106, 2 × 106, 3 × 106CL and CP for S809,
α = 14°, Re = 2 × 106
He et al. [11]Adjoint equations for solution error.
β distribution to minimize error.
VelocitiesExperiment: At several cross-sections in the flow. Cylinder flow, Re = 2 × 104; Round jet, Re = 6000; Hump flow, Re = 9.4 × 105, Wall mounted cube, Re = 105
  • Data assimilation approach.
  • Model worked well for wide range of applications.
Yang and Xiao [12]Correction term for the Transition time-scale correction, β; (Random forest, Neural Network)Transition model timescale for first mode: τnt1= βτnt1d, streamline curvature, d 2 | Ω | | U | , Q, ∇pDNS: NLF(1)-0416 airfoil, α = 0° and 4°NLF(1)-0416 airfoil, α = 2° and 6°; NACA 0012, α = 3°
  • Extended work of Duraisamy et al. for transition flow
  • The model was applied for aposteriori test both in interpolation, extrapolation and for unseen case.
Ling et al. [13]Coefficients of non-linear stress terms: g n ( λ 1 , λ 2 , λ 3 , λ 4 , λ 5 )
(Deep Neural Network)
σ = n = 1 10 g n T ( n ) ; T ( n ) = S ; S · Ω Ω · S ; ( S · S ) * ; ( Ω · Ω ) * ;   Ω · S · S S · S · Ω ; ( Ω · Ω · S S · Ω · Ω ) * ; Ω · S · Ω · Ω Ω · Ω · S · Ω ; S · Ω · S · S S · S · Ω · S ;
( S · S · Ω · Ω ) * ;
Ω · S · S · Ω · Ω Ω · Ω · S · S · Ω
*Anisotropic component
Invariants of S and Ω:
λ 1 = S : S ;   λ 2 = Ω : Ω ;  
λ 3 = ( S · S ) : S ;
λ 4 = ( Ω · Ω ) : S ;  
λ 5 = ( Ω · Ω ) : ( S · S )
DNS/LES: Duct flow, Re = 3500; Channel flow, Reτ = 590; normal (Re = 5000) and inclined (Re = 3000) jet in crossflow; Square cylinder, Re = 2000; converging-diverging channel, Reτ = 600Duct flow Re = 2000
Wavy channel, Re = 6850
  • Model applied for both apriori and aposteriori tests.
  • Machine learned model performed better than linear and non-linear physics based model, and performed well for unseen geometry and flow conditions.
Wang et al. [14]Stress prediction error:
f n = Δ τ ( q n ) ; Δτ = τRANSτDNS/LES
(Random forest regression)
k-ε RANS model:
τRANS + Δτ( q n )
10 features q n : Q, k, k d / ν , ( u . ) p , ( u . ) k , k/ε, p · p , u u : u , streamline curvature etc.DNS: Duct flow, Re = 2200, 2600 and 2900Duct flow, Re = 3500
  • Apriori estimate of the Reynolds stress
  • Model works well when trained using same geometry but different flow conditions.
  • Identify a quantitative measure for apriori estimation of prediction confidence in data-driven turbulence modeling
DNS: Periodic hill, Re = 1400, 5600Periodic hill, Re = 10,595
DNS: Wavy channel Re = 360
LES: Curved backward facing step, Re = 13,200
Wu et al. [15]DNS: Periodic hill, Re = 1400, 2800, 5600; Curved backward facing step, Re = 13200; Converging-diverging channel, Re = 11,300; Backward facing step, Re = 4900; Wavy channel, Re = 360Periodic hill, Re = 10,595
Wang et al. [16]k-ω RANS model:
τRANS + Δτ( q n )
47 features q ( n ) based on combination of S, Ω, ∇k and ∇TDNS: Flat-plate boundary layer for Ma = 2.5, 6 and 7.8, Reτ ~ 400Flat-plate boundary layer, Ma = 8
  • Extended model for compressible flow.
  • Model performed well when trained using datasets with close wall thermal characteristics
Wu et al. [17]Eddy viscosity for linear stress: v T ; Non-linear anisotropic stress: b
(Random forest regression)
σ   = v T S + b S,Ω, ∇k,∇p, k d ν , k, k/εDNS and RANS: Duct flow, Re = 2200
LES and RANS: Periodic hill—Re = 5600
Duct flow, Re = 3500, 1.25 × 105
(Shallower) Periodic hill, Re = 5600
  • Focus on addressing ill-conditioning issue reported in above studies
  • Aposteriori test for different flow conditions, and slightly different geometry
Yin et al. [18]Stress prediction error:
f n = Δ τ ( q n ) ; Δτ = τRANSτDNS(Neural Network)
k-ω RANS model:
τRANS + Δτ( q n )
47 features q ( n ) based on combination of S, Ω, ∇p and ∇kDNS: Periodic hill with different steepness, L = 3.858 α +5.142, α = 0.8, 1.2, Re = 5600 Periodic hill, α = 0.5, 1, 1.5
(Re = 5600)
  • Investigated cause of ML prediction unsmoothness and error
  • Effect of input features and grid topology were investigated
  • Validation for apriori and aposterirori (frozen stresses) tests.
  • Grid significantly effects the flow features and accuracy
Yang et al. [19]τw = f(ηi)
(Fastforward Neural Network)
Wall-modeling for Lagrangian dynamic Smagorinsky (LES) modelηi: wall parallel velocity (u||), d, grid aspect ratio and ∇pDNS: Channel flow Reτ = 1000Channel flow Re = 1000 to 1010
  • Inducing grid aspect ratio and pressure gradient did not improve results.
Weatheritt and Sandberg [20]Analytic function of anisotropic stress coefficients: βi; (Symbolic regression) σ = −2νTS+2k a ;   ν T k-ω model
a = β 1 ( I 1 , I 2 ) S + 2 ( I 1 , I 2 ) { S · Ω Ω · S } + 3 ( I 1 , I 2 ) ( S · S ) *
I 1 = S : S
I 2 = Ω : Ω
Hybrid RANS/LES: Duct flows (Re = 104, Ar = 3.3), and diffuser flow (Re = 104, Ar = 1)Duct (Re = 104, Ar = 3.3 to 1), and diffuser (Re = 5000, 104, Ar = 1.7)
  • Model tested for both apriori and aposteriori tests
  • Framework is a viable methodology for RANS
Jian et al. [21]Model coefficients: Cμ, bmn, Cmn, dmn
(Deep neural network)
RANS: σ = 2 C μ k 2 ε S + b m n k 3 ε 2 ( S · S ) *
+ c m n k 3 ε 2 { S · Ω Ω · S } + d m n k 3 ε 2 ( Ω · Ω ) *
| S | k ε DNS: Plane channel flow, Reτ = 1000, 1990, 2020, 4100 Plane channel flow, Reτ = 650, 1000, 5200
  • Validated for both apriori and aposteriori tests and in extrapolation mode.
  • ML model better captures the stress-strain misalignment in the near-wall region.
Xie et al. [22]Model coefficients C1 and C2 for mixed SGS model (Neural networks)LES, subgrid stresses and heat flux
τ = C 1 Δ 2 | S | S + C 2 Δ 2 ( u · u T )
Vorticity magnitude, velocity divergence, |   u · u T | , |S|, ∇TDNS: Compressible isotropic turbulence,
Reλ = 260, Ma = 0.4, 0.6, 0.8, 1.02
Compressible isotropic turbulence, coarser grids
  • The ML model performed better than physics-based models for aposteriori tests
Table A2. Literature review of machine learning approach for stand-alone turbulence model.
Table A2. Literature review of machine learning approach for stand-alone turbulence model.
ReferenceResponse SurfaceTurbulence ModelInput FeaturesTraining FlowsValidation CaseComments
Schmelzer et al. [23]Analytic formulation of anisotropic stress;
(Symbolic regression)
RANS: σ = f (S, S · Ω Ω · S , ( S · S ) * , S : S , Ω : Ω )S, Ω, k,τDNS: Periodic hill, Re = 1.1 × 104; Converging-diverging channel, 1.26 × 104; Curved backward-facing step 1.37 × 104Periodic hill, Re = 3.7 × 104
  • Infer algebraic stress model using symbolic regression
Fang et al. [24]Shear stress: τuv
(Deep neural network)
RANS, τuvdu/dy, Reτ, near-wall van-Driest damping, spatial non-localityDNS: Channel flow, Reτ = 550, 1000, 2000, 5200Channel flow (unseen data)
  • Combination of Reynolds number injection and boundary condition enforcement improves results
Zhu et al. [25]Turbulent eddy viscosity: νT; (Radial basis function neural network) RANS: τ = νT SU, ρ, d, d2 |Ω|, velocity direction, vorticity, Entropy, strain-rateSA RANS: NACA0012 α = 0, 10, 15, Ma = 0.15, Re = 3 × 106; RAE2822 α = 2.8, Ma = 0.73 and 0.75, Re = 6.2–6.5 × 106Airfoil flow different α
  • Different maps for near-wall, wake and far-field regions
  • Weight based on wall distance, worked
  • Training using log transformation, did not work.
King et al. [26]Stress tensor τLES, subgrid stresses τU, p, filter width, resolved rate-of-strainDNS: Isotropic and sheared turbulence Isotropic and sheared turbulence on coarse grids
  • Applied for apriori tests
  • Performed better than dynamic model
Gambara and Hattori [27]Stress tensor τ
(Feedforward neural network)
LES, subgrid stresses τDifferent sets: S, d; S, Ω, d;u, d; ∇uDNS: Channel flow, Reτ = 180, 400, 600 and 800Channel flows for unseen flow conditions.
  • Model validated for apriori tests
  • u, d input parameter was most accurate.
  • Model was more accurate than similarity model
Zhou et al. [28]u, Δ (filter width)DNS: Isotropic decaying turbulence, Reλ = 129 and 302 Isotropic decaying turbulence, Reλ = 205
  • Validated for apriori and aposteriori tests
Yuan et al. [29]Stress tensor τ
(Deconvolutional neural network)
LES, subgrid stresses τFiltered velocityDNS: Isotropic decaying turbulence, Reλ = 252Isotropic decaying turbulence, Reλ = 252
  • Validated for apriori and aposteriori tests
  • ML model performed better than physics-based models.
Maulik et al. [30]Subgrid term π
(Artificial neural network)
LES, Subgrid term πVorticity, streamfunction, rate-of-strain, vorticity gradientDNS: Decaying 2D turbulence, Re = 3.2 × 104, 6.4 × 104Decaying 2D turbulence
  • Prediction error must be considered during optimal model selection for greater accuracy
Table A3. Channel and flat-plate boundary layer database curated for ML model training.
Table A3. Channel and flat-plate boundary layer database curated for ML model training.
CaseReferenceFlow Conditions#PointsDistribution of Data Points
ReτRecSublayer,
y+ < 6
Buffer Layer,
6 ≤ y+ ≤ 40
Log-Layer,
y+ > 40
Channel Flow (DNS)
1Iwamoto et al. [39]109.4191865132229
2191.83345.56561939
3150.182681.08273132138
4297.95788.151932440128
5395.767988.022572845183
6642.5413843.31931627149
7Alamo and Jimenez [40]186.343406.974971328
8Moser et al. [41]180.563298.596172553
9392.247896.97129142391
10587.1912,485.42129111998
11Lee and Moser [42]541.23211,365.96192841142
12Alamo and Jimenez [40]546.7411,476.1129121997
13933.9620,962.511931322157
14Abe et al. [43]1016.3623,433.92241530179
15Lee and Moser [42]997.422,534.12562028207
161990.6448,563.23842130332
17Hoyas and Jimenez [44]2004.348,683.87317817291
18Bernardini et al. [45]994.722,292.11921322157
192017.448621.83841930335
204072.6105,702.45121828466
21Lee and Moser [42]5180.73137,679.27681332722
Flat-plate (DNS)
Schlatter and Orlu [46]ReθReτ#PointsSublayer,
y+ < 6
Buffer layer,
6 ≤ y+ ≤ 40
Log-layer,
y+ > 40
22670252.25505131319481
231000359.37941320480
241410492.21151320480
252000671.12405131321479
263030974.18491421478
2732701043.42721421478
2836301145.16991421478
2939701244.77421422477
3040601271.53501422478
31Jimenesz et al. [47]1100445.46853451019316
321551577.78201020315
331968690.41221021314
34Sillero et al. [48]40001306.93735351019506
3540601271.53501422499
3645001437.06601019506
3750001571.19521419502
3860001847.65441019502
3965001989.47201019502
Flat-plate (LES)
ReθReτ#PointsSublayer,
y+ ≤ 6
Buffer layer,
7 ≤ y+ ≤ 40
Log-layer,
y+ >40
40Schlatter et al. [49]670257.19643851014361
411000359.5164914362
421410491.74861015360
432150721.53411014361
442560839.55761016359
4536601162.27231116358
4641001286.70141116358
47Eitel-Amor et al. [50]50001367.35865121015487
4860001561.0621015487
4970001750.51981016486
5080001937.31131016486
5190002118.08611016486
52100002299.21191016486
53110002478.99011018486
Total19,919670 (3.4%)1134 (5.7%)18,115 (90.9%)
Table A4. Flow parameters for oscillating channel flow DNS. The non-dimensional quantities are highlighted in yellow.
Table A4. Flow parameters for oscillating channel flow DNS. The non-dimensional quantities are highlighted in yellow.
Flow ParametersHigh FrequencyMed. FrequencyLow Frequency
Baseline flow Reτ,0350
Baseline flow Rec,07250
Baseline flow centerline velocity Uc1
Half channel height H1
Kinematic viscosity ν1.38 × 10−4
Domain size3π × 2 × π
Grid192 × 129 × 192
Baseline flow P 0 / x u τ , 0 2 = 0.002331
Density ρ1
Baseline flow uτ,00.048276
α200508
α   d P 0 / d x 0.46620.116550.01865
Non-dimensional pulse frequency
ω + = ω / u τ , 0 2
0.040.010.0016
Pulse frequency ω0.675650.168910.02703
Boundary layer thickness l s = 2 ν / ω 0.20210.40421.0106
l s + = 2 u τ , 0 2 / ν ω = 2 / ω + 7.07114.14235.355
R e s = U o 2 / ω ν 100200500
R e s / l s + = U o / u τ 10 2
U o / U c 0.03296
Time step size (Δt)0.00023250.000930.000969
Timesteps per period (2π/ωΔt)4000040,000240,000
Pressure pulse d P d x = d P 0 d x [ 1 + α cos ( ω t ) ]

Appendix B. Simplification of Navier-Stokes Equation

The mass conservation equation for incompressible Navier-Stokes equation for three-dimensional flows in Cartesian coordinate is:
u x + v y + w z
where, u, v and w are the mean velocities along streamwise (x), wall-normal (y) and spanwise (z) directions, respectively. For 2D flows, the velocity and the gradients along the spanwise direction are assumed to be negligible, i.e., w = 0 and ( ) z = 0 . In addition, assuming that the flow gradients along the streamwise direction are negligible compared to those along the wall-normal direction. The mass conservation equation simplifies to:
v y = 0
For wall-bounded flows the above equation results in a solution:
v = 0
For above 2D flows, the only momentum equation which is of interest is those along the streamwise direction, i.e.,
u t + u u x + v u y + w u z = 1 ρ p x + ν ( 2 u x 2 + 2 u y 2 + 2 u z 2 ) + ( τ u u x + τ u v y + τ u w z )
Using the assumptions for 2D flows, only the time derivative term survives on the right-hand side. Among the second derivative terms, 2 u y 2 is more dominant compared to 2 u x 2 based on the gradient assumption, thus later is neglected. τ u u ,   τ u v and τ u w are the unknow turbulent stresses due to correlation between streamwise turbulent fluctuations, streamwise and wall-normal fluctuations, and streamwise and spanwise fluctuations, respectively. The term involving τ u w can be dropped because of the 2D assumptions. The term involving τ u u can be also dropped assuming that both τ u u and τ u v have similar magnitude, and gradients along the streamwise direction are negligible compared to those along the wall-normal direction. Lastly, p is the pressure in the flow, and only the streamwise pressure gradient appears in the equation. For flat-plate boundary layer flows, the streamwise pressure gradients inside the boundary layer is driven by the pressure gradients outside the boundary layer. Thus, for inner boundary layer flows the pressure gradients can be specified as a forcing term ( F ( t ) ) to account for the free-stream pressure-gradients. Thus, the unsteady simplified boundary layer equation is:
u t = F ( t ) + ν 2 u y 2 + τ u v y

References

  1. Brunton, S.L.; Noack, B.R.; Koumoutsakos, P. Machine Learning for Fluid Mechanics. Annu. Rev. Fluid Mech. 2020, 52, 477–508. [Google Scholar] [CrossRef] [Green Version]
  2. Duraisamy, K.; Iaccarino, G.; Xiao, H. Turbulence Modeling in the Age of Data. Annu. Rev. Fluid Mech. 2019, 51, 357–377. [Google Scholar] [CrossRef] [Green Version]
  3. Milano, M.; Koumoutsakos, P. Neural network modeling for near wall turbulent flow. J. Comput. Phys. 2002, 182, 1–26. [Google Scholar] [CrossRef] [Green Version]
  4. Hocevar, M.; Sirok, B.; Grabec, I. A turbulent wake estimation using radial basis function neural networks. Flow Turbul. Combust. 2005, 74, 291–308. [Google Scholar] [CrossRef]
  5. Jin, X.W.; Cheng, P.; Chen, W.L.; Li, H. Prediction model of velocity field around circular cylinder over various Reynolds numbers by fusion convolutional neural networks based on pressure on the cylinder. Phys. Fluids 2018, 30, 047105. [Google Scholar] [CrossRef]
  6. Obiols-Sales, O.; Vishnu, A.; Chandramowlishwaran, A. CFDNet: A Deep Learning-Based Accelerator for Fluid Simulations. In Proceedings of the 34th ACM International Conference on Supercomputing, Barcelona, Spain, 29 June–2 July 2020. [Google Scholar]
  7. Edeling, W.N.; Cinnella, P.; Dwight, R.P.; Bijl, H. Bayesian estimates of parameter variability in the k-ε turbulence model. J. Comput. Phys. 2014, 258, 73–94. [Google Scholar] [CrossRef] [Green Version]
  8. Ling, J.; Templeton, J. Evaluation of machine learning algorithms for prediction of regions of high Reynolds averaged Navier Stokes uncertainty. Phys. Fluids 2015, 27, 085103. [Google Scholar] [CrossRef]
  9. Parish, E.J.; Duraisamy, K. A paradigm for data-driven predictive modeling using field inversion and machine learning. J. Comput. Phys. 2016, 305, 758–774. [Google Scholar] [CrossRef] [Green Version]
  10. Singh, A.P.; Medida, S.; Duraisamy, K. Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils. AIAA J. 2017, 55, 2215–2227. [Google Scholar] [CrossRef]
  11. He, C.X.; Liu, Y.; Gan, L. A data assimilation model for turbulent flows using continuous adjoint formulation. Phys. Fluids 2018, 30, 105108. [Google Scholar] [CrossRef] [Green Version]
  12. Yang, M.; Xiao, Z. Improving the k-ω-γ-Ar transition model by the field inversion and machine learning framework. Phys. Fluids 2020, 32, 064101. [Google Scholar]
  13. Ling, J.; Kurzawski, A.; Templeton, J. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. J. Fluid Mech. 2016, 807, 155–166. [Google Scholar] [CrossRef]
  14. Wang, J.X.; Wu, J.L.; Xiao, H. Physics informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data. Phys. Rev. Fluids 2017, 2, 034603. [Google Scholar] [CrossRef] [Green Version]
  15. Wu, J.L.; Wang, J.X.; Xiao, H.; Ling, L. A priori assessment of prediction confidence for data-driven turbulence modeling. Flow, Turbul. Combust. 2017, 99, 25–46. [Google Scholar] [CrossRef]
  16. Wang, J.X.; Huang, J.; Duan, L.; Xiao, H. Prediction of Reynolds stresses in high-Mach-number turbulent boundary layers using physics-informed machine learning. Theor. Comput. Fluid Dyn. 2019, 33, 1–19. [Google Scholar] [CrossRef] [Green Version]
  17. Wu, J.L.; Xiao, H.; Paterson, E. Physics-Informed Machine Learning Approach for Augmenting Turbulence Models: A Comprehensive Framework. Phys. Rev. Fluids 2018, 3, 074602. [Google Scholar] [CrossRef] [Green Version]
  18. Yin, Y.; Yang, P.; Zhang, Y.; Chen, H.; Fu, S. Feature selection and processing of turbulence modeling based on an artificial neural network. Phys. Fluids 2020, 32, 105117. [Google Scholar] [CrossRef]
  19. Yang, X.I.A.; Zafar, S.; Wang, J.X.; Xiao, H. Predictive large-eddy-simulation wall modeling via physics-informed neural networks. Phys. Rev. Fluids 2019, 4, 034602. [Google Scholar] [CrossRef]
  20. Weatheritt, J.; Sandberg, R.D. The development of algebraic stress models using a novel evolutionary algorithm. Int. J. Heat Fluid Flow 2017, 68, 298–318. [Google Scholar] [CrossRef]
  21. Jiang, C.; Mi, J.; Laima, S.; Li, H. A Novel Algebraic Stress Model with Machine-Learning-Assisted Parameterization Energies. Energies 2020, 13, 258. [Google Scholar] [CrossRef] [Green Version]
  22. Xie, C.; Wang, J.; Li, H.; Wan, M.; Chen, S. Artificial neural network mixed model for large eddy simulation of compressible isotropic turbulence. Phys. Fluids 2019, 31, 085112. [Google Scholar]
  23. Schmelzer, M.; Dwight, R.P.; Cinnella, P. Discovery of algebraic Reynolds stress models using sparse symbolic regression. Flow Turbul. Combust. 2020, 104, 579–603. [Google Scholar] [CrossRef] [Green Version]
  24. Fang, R.; Sondak, D.; Protopapas, P.; Succi, S. Neural network models for the anisotropic Reynolds stress tensor in turbulent channel flow. J. Turbul. 2020, 21, 9–10. [Google Scholar] [CrossRef]
  25. Zhu, L.; Zhang, W.; Kou, J.; Liu, Y. Machine learning methods for turbulence modeling in subsonic flows around airfoils. Phys. Fluids 2019, 31, 015105. [Google Scholar] [CrossRef]
  26. King, R.N.; Hamlington, P.E.; Dahm, W.J.A. Autonomic closure for turbulence simulations. Phys. Rev. E 2016, 93, 031301. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Gamabara, M.; Hattori, Y. Searching for turbulence models by artificial neural network. Phys. Rev. Fluids 2017, 2, 054604. [Google Scholar] [CrossRef]
  28. Zhou, G.; He, G.; Wang, S.; Jin, G. Subgrid-scale model for large-eddy simulation of isotropic turbulent flows using an artificial neural network. Comput. Fluids 2019, 195, 104319. [Google Scholar] [CrossRef] [Green Version]
  29. Yuan, Z.; Xie, C.; Wang, J. Deconvolutional artificial neural network models for large eddy simulation of turbulence. Phys. Fluids 2020, 32, 115106. [Google Scholar] [CrossRef]
  30. Maulik, R.; San, O.; Rasheed, A.; Vedula, P. Subgrid modelling for two-dimensional turbulence using neural networks. J. Fluid Mech. 2019, 858, 122–144. [Google Scholar] [CrossRef] [Green Version]
  31. Nathan, K.J. Deep learning in fluid dynamics. J. Fluid Mech. 2017, 814, 1–4. [Google Scholar]
  32. Bhushan, S.; Burgreen, G.W.; Martinez, D.; Brewer, W. Machine Learning for Turbulence Modeling and Predictions. In Proceedings of the ASME 2020 Fluids Engineering Division Summer Meeting FEDSM2020, Orlando, FL, USA, 12–16 July 2020. [Google Scholar]
  33. Bhushan, S.; Burgreen, G.W.; Bowman, J.; Dettwiller, I.; Brewer, W. Predictions of Steady and Unsteady Flows using Machine Learned Surrogate Models. In Proceedings of the 2020 IEEE/ACM Workshop on Machine Learning in High Performance Computing Environments (MLHPC) and Workshop on Artificial Intelligence and Machine Learning for Scientific Applications (AI4S), Supercomputing 2020, Online Conference, 12 November 2020. [Google Scholar]
  34. LeCun, Y.; Bengio, Y.; Hinton, G. Deep learning. Nature 2015, 521, 436–444. [Google Scholar] [CrossRef] [PubMed]
  35. Lee, M.; Kim, H.; Joe, H.; Kim, H.-G. Multi-channel PINN: Investigating scalable and transferable neural networks for drug discovery. J. Cheminform. 2019, 11, 46. [Google Scholar] [CrossRef]
  36. Kingma, D.; Ba, J. Adam: A method for stochastic optimization. arXiv 2015, arXiv:1412.6980. [Google Scholar]
  37. Maziar, R.; Paris Perdikaris, P.; Karniadakis, G.E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 2019, 378, 686–707. [Google Scholar]
  38. Warsi, Z.U.A. Fluid Dynamics: Theoretical and Computational Approaches; CRC Press: Boca Raton, FL, USA, 1998. [Google Scholar]
  39. Iwamoto, K.; Kasagi, N.; Suzuki, Y. Direct Numerical Simulation of Turbulent channel Flow at Reτ = 2320. In Proceedings of the 6th Symposium Smart Control of Turbulence, Tokyo, Japan, 6–9 March 2005. [Google Scholar]
  40. Alamo, J.C.; Jimenez, J.; Zandonade, P.; Moser, R.D. Scaling of the Energy Spectra of Turbulent Channels. J. Fluid Mech. 2004, 500, 135–144. [Google Scholar]
  41. Moser, R.D.; Kim, J.; Mansour, N.N. Direct Numerical Simulation of Turbulent Channel Flow Up to Reτ = 590. Phys. Fluids 1999, 11, 943–945. [Google Scholar] [CrossRef]
  42. Lee, M.; Moser, R.D. Direct numerical simulation of turbulent channel flow up to Reτ = 5200. J. Fluid Mech. 2015, 744, 395–415. [Google Scholar] [CrossRef]
  43. Abe, H.; Kawamura, H.; Matsuo, Y. Surface heat-flux fluctuations in a turbulent channel flow up to Reτ = 1020 with Pr = 0.025 and 0.71. Int. J. Heat Fluid Flow 2004, 25, 404–419. [Google Scholar] [CrossRef]
  44. Hoyas, S.; Jimenez, J. Scaling of the Velocity Fluctuations in Turbulent Channels up to Reτ = 2003. Phys. Fluids 2006, 18, 011702. [Google Scholar] [CrossRef] [Green Version]
  45. Bernardini, M.; Pirozzoli, S.; Orlandi, P. Velocity statistics in turbulent channel flow up to Reτ = 4000. J. Fluid Mech. 2014, 742, 171–191. [Google Scholar] [CrossRef] [Green Version]
  46. Schlatter, P.; Orlu, R. Assessment of direct numerical simulation data of turbulent boundary layers. J. Fluid Mech. 2010, 659, 116–126. [Google Scholar] [CrossRef]
  47. Jimenez, J.; Hoyas, S.; Simens, M.P.; Mizuno, Y. Turbulent boundary layers and channels at moderate Reynolds numbers. J. Fluid Mech. 2010, 657, 335–360. [Google Scholar] [CrossRef] [Green Version]
  48. Sillero, J.A.; Jimenez, J.; Moser, R.D. One-point statistics for turbulent wall-bounded flows at Reynolds numbers up to Re~2000. Phys. Fluids 2013, 25, 105102. [Google Scholar] [CrossRef] [Green Version]
  49. Schlatter, P.; Li, Q.; Brethouwer, G.; Johansson, A.V.; Henningson, D.S. Simulations of spatially evolving turbulent boundary layers up to Reθ = 4300. Int. J. Heat Fluid Flow 2010, 31, 251–261. [Google Scholar] [CrossRef]
  50. Eitel-Amor, G.; Örlü, R.; Schlatter, P. Simulation and validation of a spatially evolving turbulent boundary layer up to Reθ = 8300. Int. J. Heat Fluid Flow 2014, 47, 57–69. [Google Scholar] [CrossRef]
  51. Martinez, D.; Brewer, W.; Strelzoff, A.; Wilson, A.; Wade, D. Rotorcraft virtual sensors via deep regression. J. Parallel Distrib. Comput. 2020, 135, 114–126. [Google Scholar] [CrossRef]
  52. Rolnick, D.; Tegmark, M. The power of deeper networks for expressing natural functions. arXiv 2018, arXiv:1705.05502. [Google Scholar]
  53. Scotti, A.; Piomelli, U. Numerical Simulation of Pulsating Turbulent Channel Flow. Phys. Fluids 2001, 13, 1367–1384. [Google Scholar] [CrossRef] [Green Version]
  54. Bhushan, S.; Walters, D.K. Development of a Parallel Pseudo-Spectral Solver Using the Influence Matrix Method and Application to Boundary Layer Transition. Eng. Appl. Comput. Fluid Mech. 2014, 8, 158–177. [Google Scholar] [CrossRef] [Green Version]
  55. Bhushan, S.; Muthu, S. Performance and Error Assessment of Parallel Pseudo-Spectra Methods for Direct Numerical Simulations. Eng. Appl. Comput. Fluid Dyn. 2019, 13, 763–781. [Google Scholar]
  56. Jamal, T.; Bhushan, S.; Walters, D.K. Numerical Simulation of Non-Stationary Turbulent Flows using Double Exponential Dynamic Time Filtering Technique. In Proceedings of the ASME 2020 Fluids Engineering Division Summer Meeting FEDSM 2020, Orlando, FL, USA, 12–16 July 2020. [Google Scholar]
  57. Muthu, S.; Bhushan, S. Temporal Direct Numerical Simulation for Flat-Plate Boundary Layer Bypass Transition. J. Turbul. 2020, 21, 311–354. [Google Scholar] [CrossRef]
Figure 1. Block diagram summarizing machine learning approach inspired from LeCun et al. [34]. The above example shows a neural network with one input, two hidden and output layers. Each of the circles represent a feature in a layer, let’s say i features in the input layer, k features in the hidden layer and m features in the output layer. The variable “x” are the input parameters, and “z” and “y” are the inputs and outputs, respectively, of both hidden and output layers. The features on the extreme right are the true values (T) used for training the network. The forward arrows represent the calculation of the features on the next layer based on linear combination of features on previous layer using the (positive) weight matrix ( w i j ). Broken arrows show backpropagation of the network prediction error to readjust the weights. The network prediction error (E) is computed using a predefined cost function comparing the network output with the true values. The derivatives of the errors are computed to adjust the weight matrix, where α is learning rate.
Figure 1. Block diagram summarizing machine learning approach inspired from LeCun et al. [34]. The above example shows a neural network with one input, two hidden and output layers. Each of the circles represent a feature in a layer, let’s say i features in the input layer, k features in the hidden layer and m features in the output layer. The variable “x” are the input parameters, and “z” and “y” are the inputs and outputs, respectively, of both hidden and output layers. The features on the extreme right are the true values (T) used for training the network. The forward arrows represent the calculation of the features on the next layer based on linear combination of features on previous layer using the (positive) weight matrix ( w i j ). Broken arrows show backpropagation of the network prediction error to readjust the weights. The network prediction error (E) is computed using a predefined cost function comparing the network output with the true values. The derivatives of the errors are computed to adjust the weight matrix, where α is learning rate.
Energies 14 01465 g001
Figure 2. Correlation of shear stress ( τ u v ) with mean velocity gradient for channel and flat-plate boundary later DNS/LES dataset. The dataset shows significant overlap for du/dy between 5 and 20 (within the box region).
Figure 2. Correlation of shear stress ( τ u v ) with mean velocity gradient for channel and flat-plate boundary later DNS/LES dataset. The dataset shows significant overlap for du/dy between 5 and 20 (within the box region).
Energies 14 01465 g002
Figure 3. (a) Apriori predictions using data-driven machine learned (DDML) using different sets of input parameters in Equation (20) for channel flow datasets. Variation of shear stress ( τ u v ) for channel flow dataset with respect to (b) R e d and (c)   R e l . The latter shows better collapse in the DNS dataset in the buffer-layer region (within the box region).
Figure 3. (a) Apriori predictions using data-driven machine learned (DDML) using different sets of input parameters in Equation (20) for channel flow datasets. Variation of shear stress ( τ u v ) for channel flow dataset with respect to (b) R e d and (c)   R e l . The latter shows better collapse in the DNS dataset in the buffer-layer region (within the box region).
Energies 14 01465 g003
Figure 4. Apriori DDML and physics-informed machine learning (PIML) stress predictions trained using channel and flat-plate datasets for channel flows at Reτ = (a) 395, (b) 590, (c) 1000, and (d) 2000 using different weighting functions and PIML are compared with DNS results.
Figure 4. Apriori DDML and physics-informed machine learning (PIML) stress predictions trained using channel and flat-plate datasets for channel flows at Reτ = (a) 395, (b) 590, (c) 1000, and (d) 2000 using different weighting functions and PIML are compared with DNS results.
Energies 14 01465 g004
Figure 5. (a) Pressure pulse, d P 0 d x α cos ( ω t ) (broken black line), and variation of wall shear stress (solid red line: top wall, solid black line: bottom wall) over five pressure oscillation cycles, and (b) instantaneous streamwise and wall-normal velocity profiles at quarter cycle obtained from DNS for the medium frequency case. As marked in subfigure (a), t = 0 and 0.5 corresponds to the peak and trough of the pressure gradient pulse, respectively. Results are shown for the medium frequency case.
Figure 5. (a) Pressure pulse, d P 0 d x α cos ( ω t ) (broken black line), and variation of wall shear stress (solid red line: top wall, solid black line: bottom wall) over five pressure oscillation cycles, and (b) instantaneous streamwise and wall-normal velocity profiles at quarter cycle obtained from DNS for the medium frequency case. As marked in subfigure (a), t = 0 and 0.5 corresponds to the peak and trough of the pressure gradient pulse, respectively. Results are shown for the medium frequency case.
Energies 14 01465 g005
Figure 6. (a) Schematic diagram demonstrating the post-processing of the 3D DNS data to obtain 1D unsteady mean flow. (b) Variation of the mean velocity over three oscillation period.
Figure 6. (a) Schematic diagram demonstrating the post-processing of the 3D DNS data to obtain 1D unsteady mean flow. (b) Variation of the mean velocity over three oscillation period.
Energies 14 01465 g006
Figure 7. Clustering of the shear stress ( τ u v ) DNS data using sklearn. The different colored regions represent unique turbulence feature.
Figure 7. Clustering of the shear stress ( τ u v ) DNS data using sklearn. The different colored regions represent unique turbulence feature.
Energies 14 01465 g007
Figure 8. (a) τ u v magnitude predicted by DNS. τ u v regression map (top) and error in the ML model (bottom) trained using (b) 0.25%, (c) 1%, (d) 10% and (e) 100% of datapoints.
Figure 8. (a) τ u v magnitude predicted by DNS. τ u v regression map (top) and error in the ML model (bottom) trained using (b) 0.25%, (c) 1%, (d) 10% and (e) 100% of datapoints.
Energies 14 01465 g008
Figure 9. Schematic diagram demonstrating coupling between machine learning turbulence model with a CFD solver.
Figure 9. Schematic diagram demonstrating coupling between machine learning turbulence model with a CFD solver.
Energies 14 01465 g009
Figure 10. (a) Convergence of DDML (input feature set ( R e l , y + , u y ) and ML-3 weighting) turbulence model simulation on coarse, medium, and fine grids for simulations started from converged one-equation RANS results. (b) Convergence of PIML turbulence model simulation on coarse grid for simulation started from ill-conditioned initial condition u = U c ( 1 d 8 ) . The abscissa title “iterations” refers to pseudo timestep level during the solution of the governing equations.
Figure 10. (a) Convergence of DDML (input feature set ( R e l , y + , u y ) and ML-3 weighting) turbulence model simulation on coarse, medium, and fine grids for simulations started from converged one-equation RANS results. (b) Convergence of PIML turbulence model simulation on coarse grid for simulation started from ill-conditioned initial condition u = U c ( 1 d 8 ) . The abscissa title “iterations” refers to pseudo timestep level during the solution of the governing equations.
Energies 14 01465 g010
Figure 11. Predictions of mean velocity (top row), turbulent (middle row) and viscous (bottom row) shear stresses for Reτ = 590 obtained using DDML (left column) and PIML (right column) using input feature set ( R e l , y + , U y ) . The DDML used ML-3 weighting. Simulations were started from an initial velocity profile obtained from one-equation RANS model. Results are compared with DNS and one-equation RANS (1-Equation Model) results.
Figure 11. Predictions of mean velocity (top row), turbulent (middle row) and viscous (bottom row) shear stresses for Reτ = 590 obtained using DDML (left column) and PIML (right column) using input feature set ( R e l , y + , U y ) . The DDML used ML-3 weighting. Simulations were started from an initial velocity profile obtained from one-equation RANS model. Results are compared with DNS and one-equation RANS (1-Equation Model) results.
Energies 14 01465 g011
Figure 12. Predictions of mean velocity (left panel) and turbulent shear stress (right column) for Reτ = 590 obtained using (a) DDML and (b) PIML turbulence models. The simulations were started from ill-posed initial flow condition.
Figure 12. Predictions of mean velocity (left panel) and turbulent shear stress (right column) for Reτ = 590 obtained using (a) DDML and (b) PIML turbulence models. The simulations were started from ill-posed initial flow condition.
Energies 14 01465 g012aEnergies 14 01465 g012b
Figure 13. Predictions of mean velocity (top row), turbulent shear stress (middle row) and turbulent kinetic energy (bottom row) for oscillating channel flow case over three pressure oscillation cycles. Results obtained using DNS (leftmost column), one-equation URANS (second column), apriori DDML predictions (third column), and aposteriori DDML (rightmost column). The URANS and DDML simulations were started from channel flow velocity profile corresponding to Reτ = 350.
Figure 13. Predictions of mean velocity (top row), turbulent shear stress (middle row) and turbulent kinetic energy (bottom row) for oscillating channel flow case over three pressure oscillation cycles. Results obtained using DNS (leftmost column), one-equation URANS (second column), apriori DDML predictions (third column), and aposteriori DDML (rightmost column). The URANS and DDML simulations were started from channel flow velocity profile corresponding to Reτ = 350.
Energies 14 01465 g013
Figure 14. Predictions of mean velocity (top row), turbulent shear stress (middle row), and turbulent kinetic energy (bottom row) for oscillating channel flow case over three pressure oscillation cycles. Results obtained using DNS (leftmost column), one-equation RANS (middle column), and aposteriori DDML (rightmost column). The URANS and DDML simulations are started from an ill-posed initial flow condition.
Figure 14. Predictions of mean velocity (top row), turbulent shear stress (middle row), and turbulent kinetic energy (bottom row) for oscillating channel flow case over three pressure oscillation cycles. Results obtained using DNS (leftmost column), one-equation RANS (middle column), and aposteriori DDML (rightmost column). The URANS and DDML simulations are started from an ill-posed initial flow condition.
Energies 14 01465 g014
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bhushan, S.; Burgreen, G.W.; Brewer, W.; Dettwiller, I.D. Development and Validation of a Machine Learned Turbulence Model. Energies 2021, 14, 1465. https://0-doi-org.brum.beds.ac.uk/10.3390/en14051465

AMA Style

Bhushan S, Burgreen GW, Brewer W, Dettwiller ID. Development and Validation of a Machine Learned Turbulence Model. Energies. 2021; 14(5):1465. https://0-doi-org.brum.beds.ac.uk/10.3390/en14051465

Chicago/Turabian Style

Bhushan, Shanti, Greg W. Burgreen, Wesley Brewer, and Ian D. Dettwiller. 2021. "Development and Validation of a Machine Learned Turbulence Model" Energies 14, no. 5: 1465. https://0-doi-org.brum.beds.ac.uk/10.3390/en14051465

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop