Turbulence closure modeling with data-driven techniques: physical compatibility and consistency considerations

A recent thrust in turbulence closure modeling research is to incorporate machine learning (ML) elements, such as neural networks, for the purpose of enhancing the predictive capability to a broader class of flows. Such a turbulence closure framework entails solving a system of equations comprised of ML functionals coupled with traditional (physics-based - PB) elements. While combining closure elements from fundamentally different ideologies can lead to unprecedented progress, there are many critical challenges that must be overcome. This study examines three such challenges: (i) Physical compatibility (or lack thereof) between ML and PB constituents of the modeling system of equations; (ii) Internal (self) consistency of the ML training process; and (iii) Formulation of an optimal objective (or loss) function for training. These issues are critically important for generalization of the ML-enhanced methods to predictive computations of complex engineering flows. Training and implementation strategies in current practice that may lead to significant incompatibilities and inconsistencies are identified. Using the simple test case of turbulent channel flow, key deficiencies are highlighted and proposals for mitigating them are investigated. Compatibility constraints are evaluated and it is demonstrated that an iterative training procedure can help ensure certain degree of consistency. In summary, this work develops foundational tenets to guide development of ML-enhanced turbulence closure models.


I. INTRODUCTION
Many studies in recent literature aim to enhance the capabilities of turbulence models by incorporating data-driven techniques. These methods seek to better utilize the increasingly available high fidelity numerical data in turbulent flows along with recent developments in machine learning (ML) to improve turbulence model predictions in unseen complex engineering flows. In these approaches, functionals obtained using ML are used within existing turbulence modeling framework to incorporate capabilities that fall beyond the scope of traditional approaches. In principle datadriven methods may be used to improve model performance at any level of turbulence closure including large-eddy simulations (LES), scale-resolving simulations (SRS) and Reynolds-averaged Navier-Stokes (RANS) methods.
In their seminal work, Sarghini et al. [1] adopted neural networks to develop improved subgrid scale models for LES of turbulence. Gamahara and Hattori [2], Maulik and San [3], Maulik et al. [4] and Beck et al. [5] proposed ML-enhanced LES closures using neural networks. Different ML techniques are also adopted to enhance RANS turbulence models. Random forest technique [6] has been used to improve RANS turbulence models by Wang et al. [7] and Wu et al. [8]. Field inversion and neural networks techniques are employed to introduce correction factors in RANS modeled transport equations by Singh and Duraisamy [9], Singh et al. [10], Zhang et al. [11] and Parish and Duraisamy [12]. Galilean invariant Reynolds stress anisotropy models are trained using tensor basis neural networks by Ling et al. [13]. Explicit algebraic Reynolds stress models have been developed using gene expression programming (GEP) by Weatheritt and Sandberg [14] and applied to different test flows by Weatheritt and Sandberg [15], Weatheritt et al. [16], Akolekar et al. [17] and Zhao et al. [18]. GEP is also applied in unsteady RANS and PANS (partially averaged Navier-Stokes) simulations by Lav et al. [19]. A review of the important studies in this area can be found in Kutz [20] and Duraraisamy et al. [21].
Incorporating data-driven techniques into the turbulence closure modeling process can have a transformative influence on the field. While preliminary studies show basis for optimism, more research is needed to understand the physical underpinnings of data-driven methods in order to ensure their generalizability to unseen test flows. In order to maximize their impact, data-driven approaches must leverage the physical understanding and closure modeling knowledge already incumbent in traditional models. Two-equation RANS models are still the most widely used tools in practical flow calculations and any improvement of their capabilities can have significant impact on engineering applications. Towards this end, the goal of this work is to investigate data-driven approaches for RANS enhancement. For the sake of clarity and contrast, through the remainder of this paper, we denote traditional closures with the prefix PB (to indicate physics-based) and the novel data-driven approach with ML (for machine learning).
PB-RANS computations of a turbulent flow involves the solution of a dynamically interacting system of equations. The two-equation RANS model is often called the lowest-order complete closure model [22] as it solves independent model equations for length and velocity scales to compute eddy viscosity. There are three main closure elements in a two-equation RANS model: algebraic (linear or non-linear) Reynolds stress constitutive relation; a modeled transport (partialdifferential) equation for kinetic energy to provide the turbulence velocity scale; and, a modeled transport (partial-differential) equation for dissipation or turbulence frequency to specify the turbulence length scale. The closure models and coefficients are typically developed in canonical flows that highlight individual turbulence processes. For reliable predictive computations of complex flows, the individual models must be independently accurate, and even more importantly, the dynamical interplay between the various equations must be compatible and consistent with overall flow physics. In PB-methods, the required compatibility between the various equations is accomplished (to the extent possible) by performing a dynamical system analysis. The fixed-point behavior at various asymptotic limits is examined for consistency with known physics. Such a systematic strategy assures some degree of generalizability to unseen flows.
Data-driven approaches for two-equation RANS proposed in literature employ machine learning (ML) methods for certain closures and retain traditional models for other aspects. For example Ling et al. [13] and Weatheritt and Sandberg [14,15] use ML for obtaining improved Reynolds stress constitutive relations while the modeled transport equations for turbulence length and velocity scales are retained without changes. On the other hand, Zhang et al. [11] and Parish and Duraisamy [12] use ML to optimize transport equation coefficients for best performance in flows of their interest without changing Reynolds stress constitutive relation. Thus the data-driven closure framework represents a mix of ML and traditional (PB) models. For generalizability of these methods to untested flows, it is important that the data-driven framework satisfy key requirements.
In this study, we examine three such prerequisites: 1. Physical Compatibility: In the PB-RANS model the various coefficients are carefully orchestrated to yield reasonable and holistic behavior in a set of canonical cases. Any MLbased modification of a subset of these coefficients can have deleterious effect on the overall computed outcome. Therefore, the importance of ensuring compatibility between ML functionals and PB elements is investigated.
2. Training Consistency: ML training requires input features that are currently obtained from baseline RANS models which employ PB-closure coefficients. Then, the training produces an improved functional for some of the same closure coefficients. In current methods, there is no process to ensure consistency between the a priori PB-closure coefficients that produce the input features and the a posteriori ML values. The consequences of such inconsistency is examined and means of imposing consistency are proposed.

Loss function formulation:
The success of the ML training approach hinges on the formulation of an appropriate loss function. The challenges in identifying the optimal flow statistics contributing to the loss function are examined.
It bears repeating that the goal of the study is not to propose a specific ML closure approach, but it is to establish fundamental guidelines for ML-RANS model development.
The paper is organized as follows. The RANS closure framework is discussed in Sec. II. Key challenges in applying ML techniques to two-equation turbulence model are identified in Sec. III.
A closed loop training framework is proposed in this section. Section IV formulates proof-ofconcept studies and demonstrates importance of defining appropriate loss function for ML. The results and inferences are presented in Sec. V. Finally, the conclusions of this study are summarized in Sec. VI.

II. RANS CLOSURE FRAMEWORK
The Navier-Stokes equations for a viscous and incompressible flow can be written as, where V i is the instantaneous velocity, p is the instantaneous pressure, ρ is the density and ν is the kinematic viscosity. Upon applying the Reynolds averaging operator [23] on the Navier-Stokes equations, RANS equations for incompressible flows are obtained: where U i is the mean velocity and P is the mean pressure. The Reynolds stress tensor ( u i u j ) in this equation is the subject of closure modeling. This symmetric, second order tensor can be decomposed into isotropic and anisotropic (b i j ) parts, where k is the turbulent kinetic energy, and δ i j is the Kronecker delta. In Reynolds stress closure modeling (RSCM) approach, modeled transport equations are solved for all independent components of Reynolds stress tensor ( u i u j ) [24][25][26][27]. At the current time, ML methods have not been employed much for enhancing RSCM closure approach.

A. Two-equation RANS
In the two-equation RANS approach, which is the subject of this study, a constitutive relationship for Reynolds stress is postulated in terms of the strain and rotation rates of the mean flow field. Modeled transport equations are solved for turbulence velocity and length scales to yield eddy viscosity.
Reynolds stress constitutive relationship: Using representation theory, a general form of the constitutive relationship for the normalized anisotropic tensor can be written as [28], in terms of ten basis tensors T (n) i j and their scalar invariant functions λ 1 , ..., λ 5 . Here ε is the turbulent dissipation. The basis tensors and scalar invariants are known functions of the normalized mean strain (s i j ) and rotation rates (r i j ) [28], where The scalar function G n of each basis tensor T (n) i j must be modeled. Through the reminder of this study the G's are referred to as constitutive closure coefficients (CCC). Different Reynolds stress constitutive relations of varying degree of complexity have been proposed in literature. The simplest constitutive relation is the Boussinesq model [22] given by, and turbulent viscosity (ν t ) can be written as, In this model, the CCC are: G 1 = −C µ = −0.09 and G n = 0 for n > 1 [29]. More complex nonlinear eddy viscosity models [30,31] and algebraic Reynolds stress models (ARSM) [28,[32][33][34] have also been proposed. Various non-linear and ARSM models determine the G's using different approaches to match equilibrium anisotropies in various flows.
The goal of ML-enhancement is to learn ML functionals for CCC using high fidelity data in flows of choice.
Modeled transport equations: The required turbulence velocity and length scales are obtained by solving the modeled transport equations for turbulent kinetic energy (k) and dissipation (ε) or specific dissipation rate (ω = ε β * k ). The standard k − ω modeled transport equations are: Here α, β , β * , σ and σ * are the transport closure coefficients (TCC). In traditional modeling, the values of the TCC are determined to satisfy known asymptotic or equilibrium behavior in canonical flows. Each calibration flow is chosen to highlight a key turbulence process.
Decaying Isotropic Turbulence (DIT): Decaying homogeneous isotropic turbulence is the simplest non-trivial turbulent flow wherein production and transport terms vanish and there is no spatial variation of the flow statistics. This case is used to determine the ratio β β * from the decay rate of turbulent kinetic energy [22]. The modeled transport equations for k and ω (Eq. (9)) reduces to: leading to the following asymptotic power-law decay of kinetic energy and turbulence frequency: In the above equation k 0 , ε 0 and ω 0 are values for k, ε and ω at the reference time t 0 = n k 0 ε 0 , respectively. It is known from a variety of experiments and DNS that the kinetic energy power-law decay exponent n is nearly a constant in the range 1.15 < n < 1.45. In standard k − ω model, the ratio β β * is determined by selecting n = 1.25, Equilibrium behavior of homogeneous turbulence: In a homogeneous flow, production term is non-zero and all transport terms are still zero. In energetic homogeneous turbulent flows three key dimensionless quantities -turbulence frequency (ω), production-to-dissipation ratio (P/ε) and mean-to-turbulence frequency ratio (Sk/ε)-evolve to their respective equilibrium states. The equilibrium values of the these quantities can be related to the unclosed model coefficients (TCC) by performing a fixed point analysis of Eq. (9): Invoking the definition of the turbulence dissipation (ε = β * kω), Eq. (13) can be simplified to Employing Eq. (8), the production term (P) can be written as, whers S is defined as S ≡ 2S i j S i j . Equation (15) can be rewritten as, Using Eqs. (14) and (16), the following relationship can be obtained, The fixed-point solutions relate α and β to the equilibrium values of Sk/ε and P/ε. From a suite of experiments and numerical simulations of homogeneous turbulence, it is known that the range of mean-to-turbulence frequency ratio is 4.0 < Sk/ε < 6.5 [29]. The production-to-dissipation ratio range is 1.5 < P/ε < 2.0 [29]. In standard k − ω model, Sk/ε = 4.13 and P/ε = 1.54 are selected to determine the values for α and β coefficients using Eqs. (14) and (17).
Body forces and other factors such as reference frame rotation and streamline curvature can change the above equilibrium behavior of turbulence. For turbulence in a rotating frame, Speziale et al. [35], demonstrate that dissipation and eddy viscosity (CCC values) must vanish asymptotically in the limit of infinite rotation. Thus there are many implied or explicit relationships between TCC and CCC. The exact nature of the relationship between TCC and CCC will depend on the order of the constitutive relation (e.g., linear, quadratic, cubic) and the complexity of flow (e.g., rotating reference frame) and type of flow.
Equilibrium behavior of log-layer: Turbulent transport model coefficients (σ and σ * ) are developed from the analysis of equilibrium boundary layer [22]: In the above κ is the Karman constant. Literature suggests that the value of this constant is in the range 0.384 < κ < 0.41 [36].
In the equilibrium log-layer, the Reynolds shear stress and kinetic energy are related as follows: [22], It has been shown β * = 0.09 leads to a log-layer solution consistent with experimental measurements [22].
Based on the above analyses, the TCC for the standard k-ω model (assuming G 1 = −0.09) are specified as: In the above, the Karman constant is taken to be 0.41 [29]. It is evident from the above discussion that the transport equation closure coefficients (TCC) and the the Reynolds stress constitutive equation coefficients (CCC) are strongly interconnected. Any change in a subset of coefficients without corresponding modification of others can lead to erroneous model behavior. It is essential that any closure procedure must make allowance for these relationships in the model development process.
While the two-equation models with advanced Reynolds stress constitutive relations perform adequate in some applications of engineering interest, the models are generally found wanting in flows which include complexities not accounted for in the model derivation procedure, e.g., non-equilibrium turbulence states, largescale unsteadiness, underlying instabilities, and spatiallydeveloping features.

III. DATA DRIVEN FRAMEWORKS
To improve the performance of two-equation RANS models in complex flows, several recent studies have considered replacing some traditional model elements by trained functionals obtained from ML. Proposed modifications include improving the Reynolds stress constitutive relations as described in [13][14][15][16][17][18][19], and optimizing coefficients in the modeled transport equations as given in [11,12,37,38]. In both instances the resulting system of equations are composed of traditional closure elements and data-derived functionals obtains from ML. In what follows we describe the current training methodology which is the open loop approach. Then, we identify shortcomings and propose potential improvements.

A. Open loop framework
The ML approaches in many current studies employ open loop training framework. Here we briefly describe the approach adopted by Ling et al. [13]. A schematic of the approach comprising of training and predictive computation stages is shown in Fig 17) and (18). This incompatibility can possibly lead to unphysical behavior of the model, even if the flow variables in the loss function behave reasonably well. To improve compatibility, one or both of the Eqs. (17) and (18) can be used to modify TCC values to be in accordance with ML values of CCC. However, over-constraining the coefficients may lead to other drawbacks. This will be discussed in the later sections. It is desirable to incorporate some form of dynamical systems analysis into the ML training process to ensure consistency and compatibility between the various coefficients. However, due to the implicit nature of the learned functionals, analytical approaches are not straightforward.
Instead, we propose embodying some degree of physical compatibility and training consistency into the modeling process by adopting a closed loop training framework. A schematic of one such framework is given in Fig. 2. Constructing the objective function in terms of Reynolds stress tensor ( u i u j ) [40] is the other option. This will certainly lead to an adequate computation of the mean velocity field. However, this can lead to another important inconsistency. As mentioned earlier, the DNS and RANS kinetic energy can be quite different. Thus the ML-functional for turbulent kinetic energy (k) and the value of k obtained from the modeled transport equation will not be the same leading to a disparity in the computed results. It is evident that some degree of disparity in the computed results is inevitable when ML-functionals and modeled transport equations are used in combination to simulate turbulent flows.
The objective function must be constructed based on the desired features of the computed flow.
Here, we designate the following features as the required elements of ML-RANS computation: 1. Accurate log-law velocity profile, 2. Accurate Reynolds shear stress ( u 1 u 2 ),
Towards this end, the objective or mean square error (MSE) loss function for ML training is defined as, where predicted outputs of the ML algorithm are denoted by b αα , u 1 u 2 and the true DNS values are shown by b αα and u 1 u 2 . Here, N represents the number of data points. In this definitions, shear stress component is normalized by true DNS value of the friction velocity (u τ ) to ensure a consistent velocity. Therefore, for the choice of flow considered in this study, i.e., channel flow, the ML objective function is defined based on Reynolds shear stress and normal anisotropy components.
Models based on fundamental physical principles can be expected to yield reasonable results for quantities not invoked in the coefficient calibration process. One of the limitations of databased methods is that the accuracy of quantities not involved in the loss (objective) function is unclear. In the computations we will examine the ability of the ML-RANS models to predict important flow quantities not used in the ML loss function: production-dissipation ratio (P/ε) and mean flow to turbulence frequency (Sk/ε).
In summary, while turbulent channel flow is a simple benchmark problem, it is ideally suited for examining important concepts on how ML can be used to enhance RANS.

B. Neural networks
Recently, various ML algorithms including neural networks and random forests have been used for modeling fluid dynamics in general and turbulence in particular. Neural networks have been shown to have superior performance in modeling non-linear and complicated relationships with high-dimensional data [13]. It has been shown that incorporating Galilean invariant features further enhances the generalizability of the neural networks [41]. In this study Tensorflow [42], which is widely used and comparatively well documented library, is employed for ML computations. A fully connected feed-forward neural network is trained using backpropagation with gradient descent method. The schematic of the selected neural network architecture is shown in Fig. 3. To control the overfitting of the neural network during the training, a L2-norm regularization term is imposed on the loss function (also known as Ridge-regression [43]) to constrain the magnitudes of the learning parameters (weights and biases of neurons). A grid search approach is adopted for hyperparameter optimization. The details of the optimized network architecture used in this work     [22].
Given that the channel flow is statistically two-dimensional, the constitutive equation we seek to train in this study requires only four basis tensors [29], Consideration is restricted to two scalar input features, λ 1 and k ων to determine the CCC from ML functional, i.e., G n = g n (λ 1 , k ων ). Datasets used in training, validation and predictive calculations are shown in Table II. For training purposes, DNS data [46] within the wall-normal distance range 0 < y/h < 0.8 is employed, where h is the channel half-width. The points in the regions near the channel center (y/h > 0. 8) are excluded to prevent the model coefficients from becoming unphysical as stresses tend to zero [47].
We examine the use of two compatibility constraints. The relationship given in Eq. (18) Case  Fig. 4. It is seen that by training the ML algorithm on multiple loops, the magnitude of the G 1 coefficient gradually reduces and finally converges to a value around 0.083. Some variation of G 1 on wall-normal distance is seen. Other coefficients converge to non-zero functions. Marked improvements in computing normal components of anisotrpy tensor are observed in ML-RANS simulations and turbulent shear stress is accurately reproduced.
The results of the baseline and ML-RANS simulations for other quantities of interest (QoI), i.e., turbulent kinetic energy (k), mean flow to turbulence frequency (Sk/ε) and production-dissipation ratio (P/ε) are nearly identical and no significant improvements are observed near the wall with ML-RANS. This is to be expected as the baseline (standard) k − ω model is tuned and calibrated to yield good agreement of these QoI.  In many applications of interest, the baseline Reynolds stress constitutive equation can be quite inadequate due to incorrect values of CCC. For such cases, it is expected that ML training with high fidelity data can lead to a ML functional for CCC that is significantly more accurate. Clearly, it is important to establish the capability of ML training procedure to recover from a poor baseline model. Rather than seek a flow in which the standard model is not correct, we simulate the scenario by intentionally modifying the standard Boussinesq model coefficient (G 1 ).
Computation results of the baseline and closed loop ML-RANS for this case are presented in Fig. 5. It can be seen that the baseline RANS simulation with modified CCC is inaccurate for most of the QoI in turbulent channel flow. By performing multiple loops of training, initially incorrect G 1 coefficient (= -0.045) used in baseline model recovers to a more reasonable value G 1 ∼ −0.074.
Other CCC converge to non-zero values similar to those in Case-1. The ML-RANS leads to improved computations of turbulent shear stress and normal anisotropy components. Other QoI such as mean velocity, Sk/ε and P/ε that have not been used in definition of the ML loss function are also significantly better than the baseline case. It should be noted that the TCC constraint for σ is imposed here. The results without this constraint are significantly worse. Thus the approach of (i) enforcing TCC constraint and (ii) closed loop training leads to recovery from an inaccurate baseline model.

Case-3: Modified TCC baseline model
In many applications, the coefficients in the modeled transport equations (TCC) can be inaccurate. To simulate this effect, in this case a key TCC value (β ) is modified from the standard value as shown in Table III. This case is of interest as the incorrect coefficient is not modified with ML training process. Nonetheless, such a scenario can occur in a practical application.
The computed results of various quantities by the baseline and the closed loop ML-RANS simulations are compared against DNS data in Fig. 6. To compensate for the altered TCC, ML process changes CCC away from their correct values. For instance, the value of G 1 drifts away from the 'correct value' of about −0.09 to about −0.42. Despite the incorrect CCC, the anisotrpoies and turbulent shear stress are captured reasonably well. Furthermore, the mean velocity and P/ε profiles are also adequately computed. It must be noted that these QoI are directly related to the ML loss function. This exhibits the strength of the closed loop training. Despite physically incompatible closure model coefficients, the training process provides reasonable prediction of QoI included in the optimization process. The incompatibility of the closure coefficients leads to poor predictions of other quantities such as Sk/ε and k. Thus, one of the important challenges in ML-RANS is to ensure reasonable behavior of QoI not related to the ML loss function. In the next subsection, we demonstrate that the behavior can be improved by imposing further TCC constraints. at Re τ = 1000 are used to perform predictive simulations of channel flow at Re τ = 5200 (Table   I)  The normal stress anisotropies are also well captured. The production-to-dissipation ratio (P/ε) is reasonably captured with all models due to the imposition of the σ constraint given in Eq. (18).
Evidently, the closed loop training leads to markedly improved results. The above findings lead to the question whether additional physical compatibility constraints can improve the predictive capability of the closed loop models.

Two compatibility constraints
In all of the closed loop ML-RANS results presented so far, only one CCC-TCC compatibility constraint (Eq. (18)) was enforced. As mentioned earlier, this constraint is critically important for two reasons: (i) obtaining the correct log-layer slope; and (ii) yielding P = ε in the log-layer. captured. The production-to-dissipation ratio is also in reasonable agreement with data. This is to be expected as a consequence of enforcement of the first compatibility constraint. The benefit of imposing the second compatibility constraint is evident from the profile of Sk/ε. There is a significant improvement in the predicted profile compared to the Case-3 one-constraint training.
Overall, the results presented in this section provides evidence of the importance of (i) closed loop training; and (ii) imposition of appropriate constraints. Although the demonstration has been provided only in the case of the simple channel flow, internal consistency (closed loop) and physical compatibility will be even more important in complex engineering flows.