A Novel Algebraic Stress Model with Machine-Learning-Assisted Parameterization

: Reynolds-stress closure modeling is critical to Reynolds-averaged Navier-Stokes (RANS) analysis, and it remains a challenging issue in reducing both structural and parametric inaccuracies. This study ﬁrst proposes a novel algebraic stress model named as tensorial quadratic eddy-viscosity model (TQEVM), in which nonlinear terms improve previous model-form failure due to neglection of nonlocal e ﬀ ects. Then a data-driven regression model based on a fully-connected deep neural network is designed to determine the TQEVM coe ﬃ cients. The well-trained data-driven model using high-ﬁdelity direct numerical simulation (DNS) data successfully learned the underlying input-output relationships, further obtaining spatial-dependent optimal values of these coe ﬃ cients. Finally, detailed validations are made in wall-bounded ﬂows where nonlocal e ﬀ ects are expected to be signiﬁcant. Comparative results indicate that TQEVM provides improvements both for the stress-strain misalignment and stress anisotropy, which are clear advantages over linear and quadratic eddy-viscosity models. TQEVM extends to the scope of resolution to the wall distance y + ≈ 9 as well as provides a realizable solution. RANS simulations with TQEVM are also carried out and the obtained mean-ﬂow quantities of interest agree well with DNS. This work, therefore, results in a high-ﬁdelity representation of Reynolds stresses and contributes to further understanding of machine-learning-assisted turbulence modeling and regression analysis.


Introduction
There are many kinds of turbulence phenomena in energy engineering, e.g., the wake meandering.To understand these complicated flow behavior, turbulence simulation and modeling is an important and useful tool, which has been extensively investigated and comprehensively reviewed [1][2][3][4][5].Although more accurate scale-resolving simulations, e.g., large-eddy simulation (LES) and direct numerical simulation (DNS), are increasingly popular, Reynolds-averaged Navier-Stokes (RANS) analysis remains a widely used option in computational fluid dynamics due to its efficiency [3] and is expanding into new stages [5,6].However, RANS analysis relies heavily on accurate closure model for Reynolds stresses, including structural (model-form) and parametric (model-parameter) adequate representations of real physics.
Various Reynolds-stress closures have been proposed.Most popular are still the linear eddy-viscosity models (EVMs), mainly including the Spalart-Allmaras model [7], SST k-ω model [8] and Realizable kmodel [9].A number of improvements have been made on these models aiming at Energies 2020, 13, 258 2 of 21 improving sensitivity to curvature and/or rotation effects, such as introducing corrections into the eddy viscosity equation [10], the dissipation rate equation [11], the specific dissipation rate equation [12], and introducing corrections into eddy-viscosity coefficient [12,13].However, none of the essential changes occur in the baseline framework of closure modeling.Obviously, the variant closure coefficients cannot improve the inherent limitations due to the underlying model-form inaccuracy rooted in Boussinesq hypothesis.
Linear stress-strain relationship origins from the Boussinesq hypothesis that turbulence react locally to the changes of mean-flow behavior in equilibrium flows.However, that is not really the case.In most turbulent flows, there exists a large misalignment in stress-strain response which has been repeatedly confirmed with high-fidelity DNS investigation [14][15][16].In view of pressure-strain correlations in the exact Reynolds-stress transport equations, the stress-strain misalignment is theoretically generated from the nonlocal response of the pressure to the entire flow field through Green's functions [17].As a result, the in-phase type EVMs certainly experience failure, especially in flows with strong spatial inhomogeneity (e.g., rotation/curvature [10,18], three-dimensionality [19], secondary-flow [20], etc.).These flow effects make changes in the mean shear (i.e., spatial inhomogeneity) and turbulence structures [21][22][23], further enhancing nonlocality on the stress-strain misalignment.Hamlington and Dahm [24] also observed that EVMs over-predict shear stresses and give inappropriate isotropic normal stresses even in simple parallel shear flows, and further explained the failure as a neglection of nonlocal effects.Edeling, et al. [25] argued that the failure of EVMs is attributed to the lack of physical interpretation of stress-strain misalignment.Although recent progress in representation of shear-stress-strain misalignment (see, e.g., [26,27]) has been made by three-equation EVMs (as an aside, misrepresentation of real physics [19,28] and increasing computational cost), these variants still fail in predictions of normal stresses, which are closely related to turbulent structures in channel flows [29] and secondary motions in non-circular ducts [30].Model-form inaccuracy cannot be reduced by corrections of model-parameter inaccuracy.
Nonlinear eddy-viscosity models (NLEVMs), referring to the second class of RANS closure in Figure 1b, is a generalization of the Boussinesq hypothesis, which was first proposed by Lumley [31] and investigated in a systematic way by Pope [32].A step forward over EVMs is that NLEVMs have the potential to improve normal stresses through added nonlinear terms (e.g., a k-ω 2 model of Wilcox and Rubesin [33]).However, no further improvements in shear stresses (same as EVMs) are obtained in most flows.According to RANS equations, the predictions of normal stresses also depend on the axial-flow motion which is controlled by the dominant shear stresses.Therefore, the key to improvements at a given NLEVM model-form returns to model-parameter optimization as EVMs do.At this point, this problem is not easy to settle down through conventional regression analysis due to following reasons (more details in Section 3): (i) lack of prior knowledge about how to select features as arguments, (ii) strong nonlinearity, resulting in difficulty in design of functional forms, (iii) high-dimensionality and computational cost, and (iv) spatial independence of the obtained coefficients.For quadratic eddy-viscosity models (QEVMs) commonly used, for instance, different values of closure coefficients [34][35][36][37][38] were introduced, however, no further improvements are obtained in most flows.
Recently, machine learning algorithms may take a breakthrough for turbulence modeling.Weatheritt and Sandberg [39] employed gene expression programming to calibrate the closure coefficients of QEVMs.Zhu, et al. [40] approximated the eddy-viscosity coefficient only existed in EVMs with a single-layer neural network and tested by flows around simple airfoils.The former obtained algebraic models with explicit analytical forms, while the latter is a data-driven model.These machine learning-based regression analyses are successful in reducing model-parameter inaccuracy.Nevertheless, normal stresses cannot be predicted correctly using the model trained by Zhu, et al. [40] because of the model-form inaccuracy of EVMs.Similarly, the model-form inaccuracy in QEVMs proposed by Weatheritt and Sandberg [39] leads to the fact that not all dominant shear stresses can be predicted correctly with the only sharing closure coefficient, thus resulting in failure in flows where Energies 2020, 13, 258 3 of 21 there are different misalignments between shear components and corresponding mean strain (e.g., secondary flow [20] and three-dimensional boundary layer [41]).Therefore, an adequate model-form representation of Reynolds stress is critical before using machine-learning-assisted regression analysis.Ling, et al. [42] adopts a more general form of NLEVMs based on 10 independent tensor bases from representation theory and then learns the closure coefficients using 5 invariants as inputs of the neural networks.Testing results show that embedded invariance properties can improve significantly the performance of the neural networks.Besides, Wang, et al. [43] applied machine learning to directly predict the Reynolds stresses discrepancies (corresponding to nonlinear terms in NLEVMs) compared to the truth, forming implicit algebraic NLEVMs.The main drawbacks of these implicit models are poor interpretability of black-box models, numerical convergence and transfer (or generalization) to another flow configuration [44].Edeling, Iaccarino, and Cinnella [25] also argued that the complex input-output mapping of black-box models are unable to increase confidence.A more comprehensive overview of data-driven turbulence modeling has been presented in Refs.[5].Recently, machine learning algorithms may take a breakthrough for turbulence modeling.Weatheritt and Sandberg [39] employed gene expression programming to calibrate the closure coefficients of QEVMs.Zhu, et al. [40] approximated the eddy-viscosity coefficient only existed in EVMs with a single-layer neural network and tested by flows around simple airfoils.The former obtained algebraic models with explicit analytical forms, while the latter is a data-driven model.These machine learning-based regression analyses are successful in reducing model-parameter inaccuracy.Nevertheless, normal stresses cannot be predicted correctly using the model trained by Zhu, et al. [40] because of the model-form inaccuracy of EVMs.Similarly, the model-form inaccuracy in QEVMs proposed by Weatheritt and Sandberg [39] leads to the fact that not all dominant shear stresses can be predicted correctly with the only sharing closure coefficient, thus resulting in failure in flows where there are different misalignments between shear components and corresponding mean strain (e.g., secondary flow [20] and three-dimensional boundary layer [41]).Therefore, an adequate model-form representation of Reynolds stress is critical before using machine-learningassisted regression analysis.Ling, et al. [42] adopts a more general form of NLEVMs based on 10 independent tensor bases from representation theory and then learns the closure coefficients using 5 invariants as inputs of the neural networks.Testing results show that embedded invariance properties can improve significantly the performance of the neural networks.Besides, Wang, et al. [43] applied machine learning to directly predict the Reynolds stresses discrepancies (corresponding to nonlinear terms in NLEVMs) compared to the truth, forming implicit algebraic NLEVMs.The main drawbacks of these implicit models are poor interpretability of black-box models, numerical convergence and transfer (or generalization) to another flow configuration [44].Edeling, Iaccarino, and Cinnella [25] also argued that the complex input-output mapping of black-box models are unable to increase confidence.A more comprehensive overview of data-driven turbulence modeling has been presented in Refs.[5].
NLEVMs also can be derived by Reynolds-stress transport models (RSTMs), referring to the third class of RANS closure in Figure 1b, and examples are models mentioned above (e.g., k-ω 2 model of Wilcox and Rubesin [33], and that of Hamlington and Dahm [24]).RSTMs are rarely used in NLEVMs also can be derived by Reynolds-stress transport models (RSTMs), referring to the third class of RANS closure in Figure 1b, and examples are models mentioned above (e.g., k-ω 2 model of Wilcox and Rubesin [33], and that of Hamlington and Dahm [24]).RSTMs are rarely used in engineering due to the high computational cost and complexity, and not further reviewed here (refer to [2][3][4]).Therefore, it is concluded that NLEVMs are increasingly popular for RANS-based turbulence modeling, especially in the age of data.
This work seeks to reduce both structural and parametric inaccuracies with a novel proposed NLEVM and machine-learning-assisted parameterization.The remainder of this article is organized as follows: In Section 2, we develop a novel algebraic stress-strain relationship as RANS closures, in which nonlinear terms can reduce model-form inaccuracy due to the stress-strain misalignment.Then a machine learning strategy is applied to the symbolic regression of the coefficients in the closure proposal in Section 3. Section 4 covers a prior validation in testing datasets and a posterior numerical simulation to evaluate model performance and followed by conclusions in Section 5.

A Novel Stress Closure Proposal
It is increasingly important to develop an adequate model-form representation of real physics, which helps to capture missed features in traditional models.As mentioned above, the failure of EVMs and NLEVMs is mainly attributed to stress-strain misalignment due to nonlocal effects [24,25].Spalart [2] also argued that modeled physics should be enhanced by nonlocal models which are categorized into two types: (i) models incorporating a scalar involving entire information, e.g., the wall distance (although it is not an invariant) in the Spalart-Allmaras model [7], and (ii) models involving high-order derivatives of the mean velocity.Here, to consider nonlocal effects with high accuracy, a stress closure model is proposed, which should fall into the second type model.
The key idea here is to consider nonlocal effects on the anisotropy evolution of the Reynolds stress tensor by constructing velocity fluctuations with nonlocal velocity gradients.Assume that the velocity increment between every two locations contributes to velocity fluctuations with a two-point correlation weight coefficient.
At an instantaneous time, the increment of flow velocity between two locations x and x can be expressed in terms of a Taylor series expansion with respect to the current location x, where From a physical point of view, turbulence develops as instability of shear layers, and the turbulent fluctuations are closely related to the velocity gradients.Thus, the velocity increment between two locations contributes to velocity fluctuations.The motion of fluid at any location in a turbulent flow is always influenced by the motion at other locations through the pressure field, i.e., nonlocal effects (see Appendix A).Then the contribution of velocities at other locations in the whole flow field to the velocity at a reference location can be accounted for through weight coefficients.An expression of velocity fluctuations u i at current location x is given as where r ≡ x − x, the space Λ is taken to be theoretically infinite, and the dimensionless weight coefficient α represents the influence intensity of location x on reference location x, varying inversely with the relative distance r m .Substituting Equation (1) into Equation (2), we obtain: where The ensemble average of coefficients α in Equation (4) represents averaging influence intensity between two points, and thus can be defined as a function of two-point correlation of flow quantities (e.g., pressure fluctuation p ), i.e., where R(x, t; r) = p (x, t)p (x + r, t).
The coefficients L m in Equation (4) represent equivalent influence distances or equivalent nonlocal length scales.The ensemble averages of these coefficients are approximately equal to themselves, i.e., L m ≈ L m , because the summation over the whole domain nearly experiences all states in all realizations.Other coefficients in Equation (3) also hold this property.
Obviously, nonlocal effects are involved in Equation (3) along with closure coefficients.The linear expansion in Equation ( 3) is suitable for flows in which the variations of the velocity gradients are not too large.When very strong spatial inhomogeneity exists (e.g., a three-dimensional turbulent Energies 2020, 13, 258 5 of 21 boundary layer), a second-order expansion may be needed.Terms up to quadratic velocity gradients in Equation (3) are sufficient for the vast majority of turbulent flows.Thus, Equation (3) reduces to: Applying an ensemble average to products of any two fluctuating components in Equation ( 7), the corresponding form for Reynolds stress is derived and further expressed with Reynolds decomposition of the flow variables into ensemble mean and fluctuating components as: where The first term in Equation ( 9) accounts for direct effects of mean velocity gradients on the anisotropy evolution, which corresponds to the rapid part Π r ij relating to nonlocal effects due to variations in the mean shear (see Appendix A), the second term represents turbulence-turbulence interactions and forces the turbulence to become isotropic, which corresponds to the slow part, Π s ij .Similar to Π s ij , the second term can be split into a purely local contribution (i.e., Ek ij / ) and the remaining nonlocal part.The former plays a dominant role in the second term in Equation ( 9), while the latter is relatively small and can be incorporated in the first term.Equation ( 9) can be rewritten as where E is a coefficient and is the dissipation rate of k.All of the terms inside the square brackets on the right side of Equation ( 13) represent nonlocal effects due to spatial variations in the mean shear (including local contributions of the mean shear) and the adjustment of turbulence itself.Equation ( 13) can be expressed with a decomposition of the velocity gradients into symmetric and antisymmetric parts, as follows: where the closure coefficients are defined as: Energies 2020, 13, 258 6 of 21 where S ij ≡ (∂u i /∂x j + ∂u j /∂x i )/2 and Ω ij ≡ (∂u i /∂x j − ∂u j /∂x i )/2, the closure coefficients B mn , C mn , and D mn involve nonlocal length scales L m , k, , and/or mean strain and rotation rate tensors.Now we discuss the last term in Equation ( 14).Since the dissipation is in reality anisotropic (but relatively weak), particularly close to solid boundaries, a form of ij similar to the formula proposed by Hanjalić and Launder [45] is chosen, in which the weakly anisotropic part is related to stress anisotropy modeled by linear EVMs.Thus, Furthermore, to guarantee that the trace of u i u j is 2k, it gets E = 1 (for incompressible flows, S kk = 0).Thus, the general formulation of Equation ( 14) simplifies to A similar analysis can be applied to Equations ( 15)-( 17).Here, the local contribution from turbulence itself to the anisotropy can be neglected due to its small contribution.Then, we obtain: where B pqkl , C pqkl , D pqkl , B mkl , C mkl , D mkl , B klm , C klm , and D klm are closure coefficients.Obviously, all coefficients in Equations ( 19)-( 22) satisfy Combining Equations ( 19)- (22) and Equation ( 23), we have: Energies 2020, 13, 258 7 of 21 Substituting Equations ( 19)-( 22) into Equation ( 8), we obtain a model that satisfies symmetry and contraction requirements.Theoretically, the corresponding coefficients can be determined through calibrations with either experimental or numerical data with the incorporation of physical consistency constraints.Rather than a "shared" coefficient as in previous models, each term of the proposed model has its unique coefficient, indicating greater ability to describe the Reynolds-stress anisotropy.
If only linear expansion is adopted in Equation ( 3) the nonlinear formulation for Reynolds stresses is expressed as Equation (19).On dimensional grounds, all of the coefficients (except for ν T ) can be expressed in terms of the length scale k 3/2 / .Thus, a new nonlinear kmodel can be obtained as: where b mn , c mn , and d mn are dimensionless closure coefficients, and the eddy-viscosity coefficient ν t here is represented by a two-equation formulation, where C µ is a dimensionless coefficient.In the standard kmodel, C µ = 0.09.According to Equation ( 13), the linear terms in Equation ( 25) describe the local relaxation of the turbulence toward isotropy while the nonlinear terms with the closure coefficients reflect nonlocal effects on anisotropy due to spatial variations in the mean shear (including local contributions of the mean shear) and turbulence itself.The mechanism of nonlocal effects incorporated in Equation ( 25) is clearly understood during its derivation.The variations in the mean shear are directly represented by the nonlinear terms as the products of mean strain and rotation rate, while the differences of turbulence structure scales in three directions lead to differences of nonlocal length scales in Equation ( 4), which can be realized by the coefficients with different values in Equation (25).
In particular, a set of scalar coefficients can be used to replace the tensor coefficients when they are identical.Furthermore, when cross terms (m n) are neglected, the proposed model in Equation (25) degenerates into a generalization of previous QEVMs (see, e.g., [34][35][36][37][38]), where c 1 , c 2 , and c 3 are dimensionless closure coefficients.It is noted that the nonlinear kmodel in Equation ( 25) is a first-order version of the proposed model, which has a similar structure to previous QEVMs in Equation (27), but with a set of tensor coefficients and extra cross terms (m n).Therefore, the tensorial quadratic eddy-viscosity model (denoted as TQEVM, hereafter) in Equation ( 25) has a greater ability to describe the anisotropy.By contrast, nonlocal effects on the anisotropy are inadequately considered in previous QEVMs and therefore TQEVM reduces model-form inaccuracy.

Model Training and Coefficients Regression
Although the proposed TQEVM has clear advantages over EVMs and QEVMs, its accuracy still depends on adequate representations of closure coefficients.When classical regression tools are applied to Equation (25), the coefficients are kept constant, e.g., spatial independence.Traditionally, different coefficients are optimized with different flow configurations.For instance, a typical value 0.09 of C µ in Equation ( 26) is determined using high-fidelity data in narrow regions close to the channel centerline.In this flow, EVMs and QEVMs return −u 1 u 2 /k = 2C µ S 12 k/ , and therefore only C µ needs to be determined.The residual coefficients could be determined in other flow configurations.
Energies 2020, 13, 258 8 of 21 However, the method fails when the number of coefficients in a closure model exceeds the number of nonzero independent Reynolds stresses, which just exists in the case with our proposed closure model.In addition, these coefficients are not globally optimal.Considering the superiority of machine learning algorithms, a deep neural network (DNN) is designed as an advanced regression tool to determine the closure coefficients in Equation (25).
The channel flow is a good alternative for training and testing the DNN-based regression model because several high-fidelity DNS databases are available.In this work, five databases at various Reynolds numbers (based on the friction velocity u τ and half-channel height h) are used to determine the closure coefficients and validate the proposed model: (i) Re τ = 590 from Moser et al. [46], (ii) Re τ = 650 from Iwamoto et al. [47], (iii) Re τ = 950, 2000, 4200 from Hoyas and Jimenez [48], (iv) Re τ = 1000, 2020, 4100 from Bernardini et al. [49], and (v) Re τ = 1000, 1990, 5200 from Lee and Moser [50].The data at Re τ = (1000, 1990, 2020, 4100) is used to train the DNN-based data-driven model for coefficients regression.Then the spatial-dependent optimal coefficients are applied to the proposed model to evaluate its performance by comparisons with DNS results at Re τ = [590, 650, 950, 2000, 4200, 5200].
The fully-connected DNN architecture used herein is shown in Figure 2a, which is a cascade of several hidden layers of nonlinear processing units.The dimensional shear parameters Sk/ = (2S ik S ik ) 1/2 k/ as mean-flow invariants at various Re τ are selected as initial inputs x 0 and the final outputs y nL = (C µ , b mn , c mn , d mn ) are the closure coefficients.Each successive layer uses the outputs from the previous layer as inputs, thus resulting in a hierarchical representation of learned features.Finally, a data-driven model is established to represent the mapping relationship between the closure coefficients and mean-flow quantities, viz, where F σ is the equivalent function representing the overall n L -layer DNN-based data-driven model and θ denotes all learned parameters, including the weights w i and the basis b i in each layer.
Energies 2020, 13, x 9 of 22 0 0 ( ; ) Where Fσ is the equivalent function representing the overall -layer DNN-based data-driven model and denotes all learned parameters, including the weights w and the basis b in each layer.Figure 2b shows a complete optimization procedure to train this DNN-based data-driven model, which is controlled by the loss function: where and are the real stress vector and reconstructed stress vector according to Equation ( 25), Figure 2b shows a complete optimization procedure to train this DNN-based data-driven model, which is controlled by the loss function: where τ and τ are the real stress vector and reconstructed stress vector according to Equation ( 25), respectively, and λ is a regularization coefficient to prevent overfitting.N is the number of training samples.The rectified linear unit (ReLU) is used as the activator.The Adam (adaptive-momentestimation) algorithm [51] is adopted to update parameters θ.Other hyper-parameters are set to λ= 10 −3 and ϕ =10 −5 (the step size).The number of hidden layers is set to 9, and the number of nonlinear units in each layer is 12, 18, 21, 27, 32, 35, 30, 28, and 27, respectively.After 5000 training epochs, iterative calculations are stopped when the tolerance error converges to δ =10 −4 .
To improve the smoothness of input-output pairs, the explicit analytic forms of F σ : x 0 → y n L is further given in Figure 3, which enhances the robustness of RANS-based numerical simulations.In traditional pointwise-based trainings (e.g., Wang et al. [43]), numerical convergence may be a problem due to the spatial non-smoothness of Reynolds stresses.In addition, the derivatives of Reynolds stresses are used in RANS equations, therefore, the smoothness is more important than the pointwise accuracy.The explicit re-parameterization adopted herein is very important otherwise alternative methods are taken to improve the smoothness.To reduce the sensitivity to re-parameterization coefficients, we select proper function forms based on three rules below: (i) That Reynolds stresses can become infinite for any possible / should be avoided [52].For example, a quadratic function is unsuitable for .As an approximate alternative, the negative exponential function of / is selected.(ii) To avoid numerical instabilities due to the high sensitivity to strain-dependent coefficients, the nested transcendental functions should be averted [53].For an exponential function, ⁄ ≈ 1 for a limited small .(iii) All coefficients should reduce to non-zero values at the high level of / to preclude numerical instability [53].According to the rules above, we select the negative exponential function, of which the advantages have been confirmed by Weatheritt and Sandberg [54].Finally, these coefficients in the Reynolds shear stress reads where the values of , , and are 30.8,250, 1.0, −0.03 and 0.22, 0, 0.41, 0.02 for ⁄ < 5.0 and ⁄ ≥ 5.0, respectively.The coefficients in the normal Reynolds stress obey To reduce the sensitivity to re-parameterization coefficients, we select proper function forms based on three rules below: (i) That Reynolds stresses can become infinite for any possible Sk/ should be avoided [52].For example, a quadratic function is unsuitable for C µ .As an approximate alternative, the negative exponential function of Sk/ is selected.(ii) To avoid numerical instabilities due to the high sensitivity to strain-dependent coefficients, the nested transcendental functions should be averted [53].For an exponential function, e a /e a+δ ≈ 1 for a limited small δ. (iii) All coefficients should reduce to non-zero values at the high level of Sk/ to preclude numerical instability [53].According to the rules above, we select the negative exponential function, of which the advantages have been confirmed by Weatheritt and Sandberg [54].Finally, these coefficients in the Reynolds shear stress reads where the values of α, β, γ and C 0 are 30.8,250, 1.0, −0.03 and 0.22, 0, 0.41, 0.02 for Sk/ < 5.0 and Sk/ ≥ 5.0, respectively.The coefficients in the normal Reynolds stress obey where ξ mn stands for the learned closure coefficients b The learned C µ agrees well with 0.09 in narrow regions far from the wall, which is the optimal value of C µ in traditional QEVMs as well as the slope of shear stress-strain corresponding to the regions Sk/ ≤ 2.2 ∼ 3.2 in Figure 4.It can be, to some extent, used as an indicator of the success of the DNN-based regression model in Figure 2. Besides, Figure 3 shows that the closure coefficients related to shear stresses perform strong linear correlations (Figure 3a), and the closure coefficients related to normal stresses behavior in a very similar pattern (Figure 3b).The fact that linear coefficient C µ reduces with increasing Sk/ (corresponding to the reduction of the wall distance) indicates an increase in the nonlinearity of shear stress when close to the wall.
Energies 2020, 13, x 11 of 22 related to normal stresses behavior in a very similar pattern (Figure 3b).The fact that linear coefficient reduces with increasing ⁄ (corresponding to the reduction of the wall distance) indicates an increase in the nonlinearity of shear stress when close to the wall.

Reynolds Stress
Figure 4 shows that Reynolds shear stress offered by all models except for QEVM-SZL falls almost exactly on the DNS results within the range of approximately Sk/ ≤ 2.2 ∼ 3.2.However, once 2.2 ∼ 3.2 < Sk/ ≤ 18 ∼ 19, the Reynolds shear stress by EVM and QEVMs seriously deviates from DNS results, while TQEVM still shows good agreement.
EVM and previous QEVMs give the same prediction for the Reynolds shear stress as −u 1 u 2 /k = C µ (Sk/ ) where C µ = 0.09 (except for QEVM-SZL).However, the TQEVM provides −u 1 u 2 /k = C µ (Sk/ ) + (b 21 − d 21 )(Sk/ ) 2 /4, in which the nonlinear terms improve predictions in the near-wall region where nonlocal effects are significant due to spatial inhomogeneity.Figure 5 shows the dimensionless shear parameter Sk/ , and dimensionless quadratic mean velocity gradient ∂ 2 u 1 /∂x 2 2 as functions of dimensionless wall distance y + .Within Sk/ ≤ 2.2 ∼ 3.2, which corresponds to a narrow region far from the wall, the variation of ∂u 1 /∂x 2 (u 1 is the axial mean velocity and x 2 is the direction normal to the wall) is small, resulting in weak nonlocal effects.Therefore, the in-phase stress-strain relationships in EVM and QEVMS achieve good performance in this region.However, in the range of 2.2 ∼ 3.2 < Sk/ ≤ 18 ∼ 19, which corresponds to a near-wall region, the dominant component ∂u 1 /∂x 2 of the mean velocity gradients vary dramatically, resulting in the large stress-strain misalignment due to nonlocal effects induced by strong spatial variations in the mean shear.As a result, the nonlinear terms in TQEVM play a significant role and improve predictions of shear stress, by contrast, EVM and QEVMs fail.( is the axial mean velocity and is the direction normal to the wall) is small, resulting in weak nonlocal effects.Therefore, the in-phase stress-strain relationships in EVM and QEVMS achieve good performance in this region.However, in the range of 2.2~3.2< ≤ ⁄ 18~19, which corresponds to a near-wall region, the dominant component ∂ ∂ ⁄ of the mean velocity gradients vary dramatically, resulting in the large stress-strain misalignment due to nonlocal effects induced by strong spatial variations in the mean shear.As a result, the nonlinear terms in TQEVM play a significant role and improve predictions of shear stress, by contrast, EVM and QEVMs fail.The normal Reynolds stress determined by various turbulence models are shown in Figures 6-8.TQEVM agrees well with DNS results in a large range of .EVM always gives the same value of 2k/3 for the normal Reynolds stress components.However, the normal stresses are in reality anisotropic even very close to the centerline of the channel.Therefore, EVM always fails.QEVM-S and QEVM-SZL offer nearly identical results with EVM, while QEVM-RB provides a slight improvement.All QEVMs under-predict the axial component of the normal Reynolds stress and over-predict the other two components.When close to the wall, the difference between normal stresses increases rapidly and the predictions by QEVMs seriously deviate from DNS results.Therefore, previous QEVMs only achieve slightly better performance than linear EVM.The normal Reynolds stress determined by various turbulence models are shown in Figures 6-8.TQEVM agrees well with DNS results in a large range of y + .EVM always gives the same value of 2k/3 for the normal Reynolds stress components.However, the normal stresses are in reality anisotropic even very close to the centerline of the channel.Therefore, EVM always fails.QEVM-S and QEVM-SZL offer nearly identical results with EVM, while QEVM-RB provides a slight improvement.All QEVMs under-predict the axial component of the normal Reynolds stress and over-predict the other two components.When close to the wall, the difference between normal stresses increases rapidly and the predictions by QEVMs seriously deviate from DNS results.Therefore, previous QEVMs only achieve slightly better performance than linear EVM.
2k/3 for the normal Reynolds stress components.However, the normal stresses are in reality anisotropic even very close to the centerline of the channel.Therefore, EVM always fails.QEVM-S and QEVM-SZL offer nearly identical results with EVM, while QEVM-RB provides a slight improvement.All QEVMs under-predict the axial component of the normal Reynolds stress and over-predict the other two components.When close to the wall, the difference between normal stresses increases rapidly and the predictions by QEVMs seriously deviate from DNS results.Therefore, previous QEVMs only achieve slightly better performance than linear EVM.

Stress-Strain Misalignment
Here we further evaluate the validity of turbulence models using modeled stress-strain misalignment.Analogous to the cosine of the angle between vectors, a new indicator is introduced as a measure of the alignment trend between the Reynolds stress and mean strain tensors, viz, where ′ and are the Reynold stress and mean strain tensors, respectively.For symmetric tensors, : ∶ = , ‖ ‖= ( ) ⁄ and ‖ ‖ = ( ) ⁄ .Therefore, is a number between 0 and 1. = 1 represents in-phase relationships between and (an alignment behavior), while = 0 corresponds to a "perpendicular" state (maximum misalignment between and ).Accordingly, can be used to characterize the validity of turbulence models.For example, EVM always returns modeled = 1 , therefore, EVM is valid only in the regions of a turbulent flow where the corresponding real values are = 1.When is far from 1, the simple linear relationships between the stress and strain are no longer satisfied, thus nonlinear terms must be added to turbulence models to correct the predicted alignment trend.
Figure 9 shows alignment indicators ( ) obtained from DNS result and predicted by previous QEVMs and TQEVM in channel flow at various .The stress-strain alignment gradually improves as increases, however, it suddenly becomes worse after reaches a very large value, which corresponds to the region very far from the wall (around the centerline of the channel).Also, the minimum phase lag (corresponding to the maximum value of ) for Reynolds stress is no less than 4 ⁄ , and the value gradually increases as the Reynolds number increases.Therefore, there is a clear stress-strain misalignment in a flow.

Stress-Strain Misalignment
Here we further evaluate the validity of turbulence models using modeled stress-strain misalignment.Analogous to the cosine of the angle between vectors, a new indicator β RS is introduced as a measure of the alignment trend between the Reynolds stress and mean strain tensors, viz, where R and S are the Reynold stress and mean strain tensors, respectively.For symmetric tensors, R :: Therefore, β RS is a number between 0 and 1.
β RS = 1 represents in-phase relationships between R and S (an alignment behavior), while β RS = 0 corresponds to a "perpendicular" state (maximum misalignment between R and S).Accordingly, β RS can be used to characterize the validity of turbulence models.For example, EVM always returns modeled β RS = 1, therefore, EVM is valid only in the regions of a turbulent flow where the corresponding real values are β RS = 1 When β RS is far from 1, the simple linear relationships between the stress and strain are no longer satisfied, thus nonlinear terms must be added to turbulence models to correct the predicted alignment trend.Figure 9 shows alignment indicators β RS (y + ) obtained from DNS result and predicted by previous QEVMs and TQEVM in channel flow at various Re τ .The stress-strain alignment gradually improves as y + increases, however, it suddenly becomes worse after y + reaches a very large value, which corresponds to the region very far from the wall (around the centerline of the channel).Also, the minimum phase lag (corresponding to the maximum value of β RS ) for Reynolds stress is no less than π/4, and the value gradually increases as the Reynolds number increases.Therefore, there is a clear stress-strain misalignment in a flow.Previous QEVMs always provide in-phase relationships for the Reynolds shear stress as EVM (nonlinear terms in QEVMs have no effects on the shear stress) and a mild separation of the normal Reynolds stress.Consequently, these QEVMs cannot accurately predict .However, the results from TQEVM agree well with DNS, owing to nonlinear relationships both for the Reynolds shear stress and the normal Reynolds stress.

Realizability Constraints
Also, the realizability of turbulence models are examined using the barycentric map [55], which is a linear domain providing a non-distorted visual representation of anisotropy, as an advantage over commonly used anisotropy invariant map [56].The eigen decomposition of the Reynolds stress tensor reads where are eigenvectors and Λ = diag( , , ) is the diagonal matrix of eigenvalues ( ≥ ≥ , + + = 0).Here, turbulent kinetic energy, eigenvectors, and eigenvalues represent the magnitude, orientation, and shape of the Reynolds stress, respectively.In the anisotropy phase space of the barycentric map, any anisotropy state is located as a point such that the linear relationship holds: ( , ) = ) correspond respectively to one component (1C), two-component axisymmetric (2C-axis), and three-component isotropic (3C-iso) states.The corresponding weights ( , , ) are entirely Previous QEVMs always provide in-phase relationships for the Reynolds shear stress as EVM (nonlinear terms in QEVMs have no effects on the shear stress) and a mild separation of the normal Reynolds stress.Consequently, these QEVMs cannot accurately predict β RS .However, the results from TQEVM agree well with DNS, owing to nonlinear relationships both for the Reynolds shear stress and the normal Reynolds stress.

Realizability Constraints
Also, the realizability of turbulence models are examined using the barycentric map [55], which is a linear domain providing a non-distorted visual representation of anisotropy, as an advantage over commonly used anisotropy invariant map [56].The eigen decomposition of the Reynolds stress tensor reads where v ik are eigenvectors and Λ kl = diag(λ 1 , λ 2 , λ 3 ) is the diagonal matrix of eigenvalues (λ 1 ≥ λ 2 ≥ λ 3 , λ 1 + λ 2 + λ 3 = 0).Here, turbulent kinetic energy, eigenvectors, and eigenvalues represent the magnitude, orientation, and shape of the Reynolds stress, respectively.In the anisotropy phase space of the barycentric map, any anisotropy state is located as a point such that the linear relationship holds: (x, In Figure 10, all the QEVMs except for the proposed TQEVM give unrealizable solutions for Reynolds stresses with C 2C > 1 and C 3C < 0 when close to the wall.Moreover, QEVMs fail in most regions with large discrepancies of C 1C between prediction and truth, and they also provide a larger C 3C in low Reynolds numbers and a larger C 2C in high Reynolds numbers.When close to the channel centerline, all models return to a 3C-iso state with modeled C 3C = 1 as the strain vanishing, which conflicts with a mild anisotropic state in reality.As a result, a corresponding sudden descent of stress-strain misalignment in Figure 9 arises.It is noted that this circumstance is encountered in all algebraic models (not limited to models mentioned here), but fortunately, the influence can be ignored because of extremely low anisotropic degree and only a small region existing (shown in Figures 6-8).
Energies 2020, 13, x 16 of 22 determined by the eigenvalues as = − , = 2( − ), and = 3 + 1.Thus, the weights are used herein to characterize the degree of stress anisotropy and any realizable turbulence should satisfy 0 ≤ , , ≤ 1.Note that weights are very sensitive to a slight difference of Reynolds stress components, especially the values of and (amplitude modulation with a multiplier of 2 or 3).
In Figure 10, all the QEVMs except for the proposed TQEVM give unrealizable solutions for Reynolds stresses with > 1 and < 0 when close to the wall.Moreover, QEVMs fail in most regions with large discrepancies of between prediction and truth, and they also provide a larger in low Reynolds numbers and a larger in high Reynolds numbers.When close to the channel centerline, all models return to a 3C-iso state with modeled = 1 as the strain vanishing, which conflicts with a mild anisotropic state in reality.As a result, a corresponding sudden descent of stressstrain misalignment in Figure 9 arises.It is noted that this circumstance is encountered in all algebraic models (not limited to models mentioned here), but fortunately, the influence can be ignored because of extremely low anisotropic degree and only a small region existing (shown in Figures 6-8).Form the above discussion, it can be concluded that the proposed TQEVM has clear advantages over EVM and previous QEVMs.TQEVM offers a more approximate description both for stressstrain misalignment and stress anisotropy as well as extends the scope of the resolution.Form the above discussion, it can be concluded that the proposed TQEVM has clear advantages over EVM and previous QEVMs.TQEVM offers a more approximate description both for stress-strain misalignment and stress anisotropy as well as extends the scope of the resolution.

A Posterior Numerical Simulation
Four cases of fully-developed channel flows are selected for numerical comparisons: Case I with Re τ = 650 [46], Case II with Re τ = 1000 [49], Case III with Re τ = 5200 [50] and Case IV with Re τ = 8000 [57].Case I and III are investigated to show the Reynolds-number extrapolation for which DNS data are readily available to examine the prediction performance.Case IV is chosen to demonstrate the capability of the proposed TQEVM at higher Reynolds number where DNS is computationally expensive and high-fidelity data are extremely rare.Besides, Re τ = 8000 is such a special case that the logarithmic region of the streamwise Reynolds stress exists and does not overlap with that of the axial mean velocity [57].In the mean velocity profile, an incipient logarithmic region appears until Re τ = 4000 and an unambiguous logarithmic region at Re τ = 5200.
Two-dimensional RANS simulations with TQEVM are performed using structured meshes due to homogeneity in the spanwise direction.The incompressible flow solver simpleFOAM [58] based on OpenFOAM (an open-source package) is adopted.The height of the computation domain is 2h and the length is 8πh (h is the half-channel height).Figure 11 shows the predicted axial mean velocity profile, where all quantities are normalized by the half-channel height and axial mean velocity at the channel centerline.The mean velocity here is averaged over three locations randomly selected in the streamwise direction.The results show that the axial mean velocity predicted by TQEVM is basically consistent with DNS results.Two-dimensional RANS simulations with TQEVM are performed using structured meshes due to homogeneity in the spanwise direction.The incompressible flow solver simpleFOAM [58] based on OpenFOAM (an open-source package) is adopted.The height of the computation domain is 2h and the length is 8 ℎ (h is the half-channel height).Figure 11 shows the predicted axial mean velocity profile, where all quantities are normalized by the half-channel height and axial mean velocity at the channel centerline.The mean velocity here is averaged over three locations randomly selected in the streamwise direction.The results show that the axial mean velocity predicted by TQEVM is basically consistent with DNS results.Furthermore, skin friction coefficients predicted by TQEVM agree well with DNS as shown in Table 2.Note that the skin friction coefficients at = 650, 1000 and 5200 are from DNS of Schultz and Flack [59] and that of = 8000 is estimated according to classical Prandtl's smooth flow formula with a set of log-law coefficients suggested by Bernardini, Pirozzoli and Orlandi [49], viz,   Furthermore, skin friction coefficients predicted by TQEVM agree well with DNS as shown in Table 2.Note that the skin friction coefficients at Re τ = 650, 1000 and 5200 are from DNS of Schultz and Flack [59] and that of Re τ = 8000 is estimated according to classical Prandtl's smooth flow formula with a set of log-law coefficients suggested by Bernardini, Pirozzoli and Orlandi [49], viz, where

Conclusions
In the work, we present a machine learning strategy to assist the development of RANS closure modeling, aiming at reducing both structural (model-form) and parametric (model-parameter) inaccuracies.First, a novel algebraic Reynolds stress model named as TQEVM is developed, where nonlocal effects on the anisotropy evolution of Reynolds stresses are accounted for by nonlinear terms with corresponding coefficients, thus achieving a more appropriate description of Reynolds stresses.The proposed TQEVM also guarantees the symmetry and contraction requirements.Furthermore, a data-driven model by fusion of a fully-connected DNN architecture is designed to determine the coefficients of the proposed TQEVM.The well-trained data-driven model using high-fidelity DNS data at various Reynolds numbers successfully learned the underlying relationships between the closure coefficients and mean-flow quantities, resulting in spatial-dependent optimal values of these coefficients.The learned coefficient related to the linear terms C µ agrees well with the real slope of shear stress-strain in regions close to the channel centerline.Thus, machine-learning-assisted regression analysis used herein offers a new route to reducing model-parameter inaccuracy, further advancing the performance of a given model.
Finally, a complete validation including a prior and posterior testing is performed in detail to evaluate the performance of TQEVM.Comparative results show that the proposed TQEVM has clear advantages over traditional EVM and QEVMs in predictions both for the Reynolds shear stress and normal Reynolds stress by a wide margin.Besides, the proposed TQEVM extends the scope of resolution to y + ≈ 9 as well as provides a realizable solution.Nonlocal effects due to spatial inhomogeneity are significant in the near-wall region, resulting in a large stress-strain misalignment and strong stress anisotropy.Linear EVM and previous QEVMs misrepresent the real misalignment trend and under-predict the anisotropy magnitude.The failure is in large part due to insufficient consideration of nonlocal effects, which is a key source of error in previous RANS models.By contrast, the effects of nonlinear terms in the proposed TQEVM are dramatic, which contributes to improvements both for stress-strain misalignment and stress anisotropy.Accurate mean-flow quantities of interest are also obtained using numerical simulations with the proposed TQEVM.
It should be noted that TQEVM is expected to achieve improved predictions in other flows with strong spatial inhomogeneity, such as flows with secondary motion, three-dimensionality or curvature/rotation.TQEVM will be tested in these flows in the near future.Also, our model brings a new vision about how to select input features of a data-driven model for directly predicting Reynolds stresses: the second and higher-order derivatives of the mean velocity field should be introduced as supplements.
where D/Dt ≡ ∂/∂t + u k ∂/∂x k is the material derivative; the overbar denotes an ensemble average; u i and u i are the respective mean and fluctuating velocity components in the x i direction; p is the kinematic pressure; ν is the kinematic viscosity; and the viscous dissipation tensor ij , pressure-strain correlation tensor Π ij and turbulent-transport tensor C ijk are defined as ), C ijk ≡ u i u j u k + p u i δ jk + p u j δ ik .(A2) The pressure fluctuation p appearing in Π ij is governed by a Poisson equation: Thus, the pressure p can be obtained by an integral over the entire spatial domain of the flow with Green's functions.The integral results in the pressure at a given point being influenced by the whole flow field, i.e., nonlocal effects.It indicates that nonlocality is intrinsic for incompressible flows, i.e., the pressure field responds nonlocally to changes in the flow to maintain incompressibility.Furthermore, the pressure p in Equation (A3) is obviously influenced by three parts: mean-flow field, turbulence itself and solid-wall boundaries.These correspond to "rapid", "slow" and wall parts as p ≡ p (r) + p (s) + p (w) , governed by their respective Poisson equations: The effect of wall proximity from p (w) is significant only in the very near-wall region [60], which is neglected here.
Based on the "rapid" and "slow" parts, the rapid and slow terms in the pressure-strain correlation Π ij [17] are: The slow part Π s ij represents changes caused by changes in turbulence itself and reflects "slow" relaxation of the turbulence towards isotropy.The most common representation for Π s ij is the return-to-isotropy [61], which has a linear relation.Sarkar et al. [62] argued that the additional quadratic terms should be included in Π s ij , while it has also been pointed out that these terms are small [63].Therefore, the nonlocal effects in Π s ij have little influence on the anisotropy.where r m ≡ xm − x m .It is common to consider the mean velocity gradients sufficiently homogeneous (i.e., A kl (x) ≈ A kl (x)) to be brought outside the integral [17,[64][65][66].However, due to strong spatial variations of mean-flow quantities in the wall-normal direction, nonlocal effects in Π r ij depends on resulting from p (r) are considerable in wall-bounded flows.Other inhomogeneous flows (e.g., free shear flows) also have a certain degree of nonlocality.Thus, A kl (x) − A kl (x) in Equation (A7) is not equal to zero and generates nonlocal effects on flow behaviors.
Based on the analysis above, the pressure-strain correlation Π ij in Equation (A3) is not a localized process.The nonlocality of pressure has effects on the anisotropy evolution of the Reynolds-stress tensor, and Π r ij contributes much more than the other two parts to nonlocal effects on the anisotropy.Therefore, a more accurate anisotropy representation of Reynolds stresses should incorporate nonlocal effects due to spatial variations in mean shear.

Energies 2020, 13 , x 3 of 22 Figure 1 .
Figure 1.A schematic showing the hierarchy of turbulence simulation and modeling at the level of Navier-Stokes equations: (a) four classes of simulation methods based on modeling requirements and computational cost, and (b) three classes of turbulence models for Reynolds-averaged Navier-Stokes (RANS) closures which implied two sources of model inaccuracies due to structural (model-form) and parametric (model-parameter) inadequate representations of real physics.

Figure 1 .
Figure 1.A schematic showing the hierarchy of turbulence simulation and modeling at the level of Navier-Stokes equations: (a) four classes of simulation methods based on modeling requirements and computational cost, and (b) three classes of turbulence models for Reynolds-averaged Navier-Stokes (RANS) closures which implied two sources of model inaccuracies due to structural (model-form) and parametric (model-parameter) inadequate representations of real physics.

Figure 2 .
Figure 2. A schematic of the designed data-driven model to determine the closure coefficients: (a) the deep neural network (DNN) architecture to map the initial inputs to the finial outputs, and (b) the optimization procedure to train this DNN-based data-driven model in (a).

Figure 2 .
Figure 2. A schematic of the designed data-driven model to determine the closure coefficients: (a) the deep neural network (DNN) architecture to map the initial inputs to the finial outputs, and (b) the optimization procedure to train this DNN-based data-driven model in (a).

Energies 2020, 13 , x 10 of 22 Figure 3 .
Figure 3. Closure coefficients as a function of dimensionless shear parameter ⁄ in fullydeveloped turbulent channel flow: (a) coefficients appearing in the Reynolds shear stress , and , and (b) coefficients appearing in the normal Reynolds stress , , , , and .The wall distance reduces with increasing ⁄ when ⁄ < 18.

Figure 3 .
Figure 3. Closure coefficients as a function of dimensionless shear parameter Sk/ in fully-developed turbulent channel flow: (a) coefficients appearing in the Reynolds shear stress C µ , b 21 and d 21 , and (b) coefficients appearing in the normal Reynolds stress b 11 , b 22 , c 11 , c 22 , d 11 and d 22 .The wall distance reduces with increasing Sk/ when Sk/ < 18.

Energies 2020 ,
13, x 12 of 22 shows the dimensionless shear parameter ⁄ , and dimensionless quadratic mean velocity gradient ∂ ∂ ⁄ as functions of dimensionless wall distance .Within ⁄ ≤ 2.2~3.2, which corresponds to a narrow region far from the wall, the variation of ∂ ∂ ⁄

Figure 5 .
Figure 5. (a) Dimensionless shear parameter ⁄ and (b) dimensionless quadratic mean velocity gradient ⁄ as functions of dimensionless wall distance in channel flow at = 650, 1000, 2000 and 5200.and is the axial mean velocity, and is the direction normal to the wall.

Figure 5 .
Figure 5. (a) Dimensionless shear parameter Sk/ and (b) dimensionless quadratic mean velocity gradient d 2 u 1 /dx 2 2 as functions of dimensionless wall distance y + in channel flow at Re τ = 650, 1000, 2000 and 5200.u 1 and is the axial mean velocity, and x 2 is the direction normal to the wall.

Figure 11 .
Figure 11.Comparison of solved axial mean-velocity by TQEVM (dash line) and DNS (solid line) for four different cases of channel flow:= 650, = 1000, = 5200 and = 8000 (from left to right), where h is the half-height of the channel and is axial mean velocity at the channel centerline.

Figure 11 .
Figure 11.Comparison of solved axial mean-velocity U by TQEVM (dash line) and DNS (solid line) for four different cases of channel flow: Re τ = 650, Re τ = 1000, Re τ = 5200 and Re τ = 8000 (from left to right), where h is the half-height of the channel and U c is axial mean velocity at the channel centerline.
[50] cases of fully-developed channel flows are selected for numerical comparisons: Case I with = 650[46], Case II with = 1000[49], Case III with = 5200[50]and Case IV with = 8000 [57].Case I and III are investigated to show the Reynolds-number extrapolation for which DNS data are readily available to examine the prediction performance.Case IV is chosen to demonstrate the capability of the proposed TQEVM at higher Reynolds number where DNS is computationally expensive and high-fidelity data are extremely rare.Besides, = 8000 is such a special case that the logarithmic region of the streamwise Reynolds stress exists and does not overlap with that of the axial mean velocity [57].In the mean velocity profile, an incipient logarithmic region appears until = 4000 and an unambiguous logarithmic region at = 5200.

Table 2 .
Comparison of skin friction coefficients by TQEVM and DNS for four given cases.

Table 2 .
Comparison of skin friction coefficients by TQEVM and DNS for four given cases.
Case I: Re τ =650Case II: Re τ =1000 Case III: Re τ =5200 Case IV: Re τ =8000 The rapid part Π r ij depends on the mean velocity gradients ∂u k /∂ xl , whose variations directly affects Reynolds stress.The spatial variations of the mean velocity gradient A kl (x) ≡ ∂u k /∂ xl can be expressed in terms of a local Taylor expansion at A kl (x) ≡ ∂u k /∂x l :A kl (x) − A kl (x) = r m