{'fraction_non_alphanumeric': 0.029096661313903604, 'fraction_numerical': 0.007789468638400882, 'mean_word_length': 4.919137466307277, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An iterative data-driven turbulence modeling framework based on Reynolds stress representation

Data-driven turbulence modeling studies have reached such a stage that the fundamental framework is basically settled, but several essential issues remain that strongly affect the performance, including accuracy, smoothness, and generalization capacity. Two problems are studied in the current research: (1) the processing of the Reynolds stress tensor and (2) the coupling method between the machine learning turbulence model and flow solver. The first determines the form of predicting targets and the resulting physical completeness and interpretability. The second determines the training process and intrinsic relevance between the mean flow features and Reynolds stress. For the Reynolds stress processing issue, we perform the theoretical derivation to extend the relevant tensor arguments of Reynolds stress in addition to the strain rate and rotation rate. Then, the tensor representation theorem is employed to give the complete irreducible invariants and integrity basis. In addition, an adaptive regularization term is employed to enhance the representation performance. For the coupling issue, an iterative coupling data-driven turbulence modeling framework with consistent convergence is proposed. The training data preparation, predicting target selection, and computation platform are illustrated. The framework is then applied to a canonical separated flow for verification. The mean flow results obtained by coupling computation of the trained machine learning model and flow solver have high consistency with the direct numerical simulation true values, which proves the validity of the current approach.


I. Introduction
Modern engineering design requires high accuracy of flow separation prediction for the 2 computational fluid dynamics (CFD) computation. For the complex turbulent flow separation problem, the traditional turbulence simulation methods either produce unsatisfactory flow prediction or require large computational costs, which cannot meet the requirements of accuracy and efficiency. With the rapid development of data science and machine learning (ML) techniques, the influences of flow structure and physical features that are ignored or difficult to consider in traditional turbulence modeling can be extracted and mapped to the turbulence quantities. The obtained augmented turbulence is referred to as data-driven turbulence modeling. Duraisamy et al. [1] demonstrated the data-driven turbulence modeling as a mathematical expression. First, the baseline turbulence model can be expressed as M: (2) The expression shows four roles that the data can play. (1) Extend the set of independent variables w [2,3]. (2) Modify certain terms in the governing equations P () [4][5][6][7][8].
Different choices of correction terms reflect the researchers' view of where the main discrepancy is located and correspond to different upper limits of augmentation. Regardless of the direction, the final obtained model can be regarded as a new constitutive relation that can predict Reynolds stress and mean flow quantities closer to the true values.

A. Classification of predicting targets
The complexity of the Reynolds stress leads to different strategies for selecting the predicting targets, each having its advantages and limitations. Four commonly employed choices are separately reviewed as follows from the aspects of spatial invariance, physical interpretability, and smoothness. 3 (1) The first choice is the multiplying factor of the turbulence model equation term, marked as β(x) [4][5][6][7][8][26][27][28][29][30]. As a scalar field, predicting β can ensure spatial invariance. By multiplying a spatially non-uniform factor field β(x), the turbulent quantities in certain areas are altered and the overall performance of the turbulence model can be enhanced. The statistical inference method is generally employed to solve the inverse problem. After acquiring the β field, a machine learning model is trained to predict the β field using mean flow features η. The entire process is named field inversion machine learning (FIML) [6][7][8]. However, the inferred β field can yield a quite close result, but the Reynolds stress might be far from the truth because there might exist many different turbulence fields all producing the same mean flow field. The inference cannot guarantee the correct turbulence field.
(2) The second choice is the eddy viscosity νt. Also, as a scalar, predicting νt satisfies the spatial invariance. A similar problem as β is how to acquire the "correct" νt. There are two methods. The first is employing statistical inference. By changing the inferred variable from the turbulence model term to the turbulence model result, the optimal νt field can be directly acquired [31]. However, the inferred νt faces the same problem as FIML. The second is to compute the optimal viscosity using the pointwise least-square approximation [3,21]. However, the result might lose clear physical implications and smoothness in the complex flow region where the Reynolds stress anisotropy is remarkable, which deteriorates the model performance.
(3) The third choice involves the Reynolds stress eigenvalues (k, λ1, λ2) and eigenvectors (v1, v2, v3) [2,3,32]. Selecting such targets means discarding all the assumptions and modeling the entire stress as a second-order symmetric tensor. These features can be computed by eigendecomposition of the true value from high fidelity databases, or by inference from observed mean flow quantities [33].
The eigendecomposition method needs to deal with the spatial invariance problem. The invariance of three scalars (k, λ1, λ2) can be guaranteed but the vectors (v1, v2, v3) are naturally spatially variant. One solution is to introduce the baseline eigenvectors and change the targets to the discrepancy between two sets of vectors [3,32]. Such a treatment introduces spatial rotation invariance, but reflection invariance is still missing. In addition, the rotation angle faces discontinuity because of the switching of the eigenvalue ranking and needs further numerical treatment [32].
(4) The fourth choice is the Reynolds stress representation based on the tensor function representation theorem. This method comes from the nonlinear eddy viscosity model (marked as NEVM below) in traditional turbulence modeling. Researchers supposed the Reynolds stress as a tensor function of strain rate S and rotation rate Ω, which is: In addition, prior physics of turbulence modeling require the tensor function to be isotropic under the extended Galilean transformation [34], which means that the symmetric transformation group of the Reynolds stress function is the entire full orthogonal group (rotation and reflection), which is: where Q is a temporal-constant orthogonal matrix.
Pope [35] deduced 10 tensor bases (referred to as the integrity basis) and 5 invariants formed by S and Ω using the Cayley-Hamilton theorem. Any symmetric isotropic tensor can be obtained by the linear combination of the 10 tensor bases and the coefficients are functions of the 5 invariants.
One clear advantage of the Reynolds stress representation is a combination of the accuracy and realizability. Taking multiple tensor bases evades the poor performance of t m  . Meanwhile, the coefficients are all scalars, naturally guaranteeing spatial invariance. Therefore, current research follows this direction.

B. The coupling method: frozen and iterative
After the predicting targets are selected, we need to select a coupling method between the ML model and the CFD solver, which can be divided into two categories: frozen substitution and iterative substitution.
In frozen substitution, the ML model establishes the mapping from the mean flow features computed by the baseline model to the Reynolds stress true value. Therefore, when used for prediction, the baseline Reynolds-averaged Navier-Stokes (RANS) simulation is performed to acquire the baseline mean features. The model is executed once, and the predicted value corresponds to the true stress. The stress is then substituted into the RANS equations and frozen until convergence. The flow chart is shown in Fig. 1 (a).
In the iterative substitution, the model establishes the mapping from the true mean flow features obtained from the high-fidelity database, e.g. the direct numerical simulation (DNS) results, to the true Reynolds stress. When used for prediction, using the baseline input features leads to incorrect stress at the beginning. Therefore, iteration is needed, and the model is executed in each CFD iteration step. After the iteration converges, the mean flow and the Reynolds stress both converge to correct results. The flow chart is shown in Fig. 1 [21], and Duraisamy [37]. Given that there is no unified statement of the issue, the problem is referred to as "truth compatibility" in the current research.
Second, the coupling affects convergence and accuracy. Inappropriate Reynolds stress treatments may cause ill-conditioned problems [38], which means that the computed mean flow given the true stress is different from the true mean flow. Relevant studies reached a consensus that the iterative framework has better accuracy and decreases the propagation error. However, the iterative framework may encounter convergence difficulty because the model needs to take the intermediate flow quantities as input and ensure that the iteration ends in the final correct results, which requires the ML model to have strong dynamic robustness.
Third, importantly, the coupling method determines the physical relevance between the input and the stress. In the frozen framework, the mapping from baseline mean flow to true stress lacks sufficient physical implications. The excavated rule by the ML is more like identifying the baseline error region and modifying it, but not a constitutive relation. In contrast, in the iterative framework, the mapping is established from true mean flow to true stress. The 6 physical relevance is more reasonable, and the excavated rule is closer to the nature of turbulence.
In summary, the two coupling methods each have their advantages and limitations. In the current research, since we have selected the Reynolds representation coefficients as targets, the consistency between the mean flow and the stress is more necessary and should be guaranteed first. Therefore, we selected the iterative framework.

C. Motivation
In the present work, we construct an iterative data-driven turbulence modeling framework based on Reynolds stress representation. Two main processes have been developed. First, we review the selection of tensor arguments that the Reynolds stress depends on and further discuss the physical implications. We then reformulate the tensor invariants and integrity basis using the tensor function representation theorem. The results under two-dimensional flow and threedimensional flow are given out separately.
Second, a novel iterative framework is proposed. The framework is designed to manifest "consistent convergence" including the truth compatibility and the dynamic robustness mentioned above. In addition, the framework separates the ML process and the Reynolds stress representation, while nearly all earlier studies combined them. The treatment greatly enhances the physical interpretability and smoothness of the coefficients and the final prediction performance.
The rest of this paper is organized as follows. Section II introduces the methodology from three aspects: tensor representation analysis, framework construction, and representation coefficient computation. Section III presents the numerical results of ML model training and prediction of the canonical periodic hill flow. Section IV discusses the effect of the remaining part after the Reynolds stress representation and the ML model calling frequency. Section V summarizes the paper.

II. Methodology
The three-dimensional compressible RANS equations for a Newtonian fluid without body force, heat transfer, and heat generation are where λ and μ are the bulk viscosity and molecular viscosity, respectively. T represents the total stress tensor, including the pressure, the molecular viscous stress, and the Reynolds stress τ that must be closed.
Regardless of traditional modeling or data-driven modeling, τ can always be expressed as follows: ( ) To build a new constitutive relation, the first step is to select proper arguments. Then the tensor invariants and integrity basis can be acquired, which leads to the first part of the methodology.

Extension of the tensor arguments
Tensor representation analysis has been employed in turbulence modeling for decades. As mentioned above, since Lumley [39] and Pope [35] proposed a complete form of the NEVM, with the development of subsequent studies [40][41][42], the model using (S, Ω) as tensor arguments is currently being perfected. In the current research, we start from the basic nonlinear eddy viscosity model and explore other potential tensor arguments apart from (S, Ω). More specifically, we evaluate the hypotheses of the NEVM and introduce additional quantities when these hypotheses do not hold.
The original form of the NEVM deduced by Lumley [39] is as follows: where q is the turbulent velocity scale, and ξ = x/(q 3 /ε) is the nondimensional spatial coordinate normalized by q and the dissipation rate ε. Three assumptions are used during the deduction: 8 (1) incompressible turbulence, (2) weak historical effect.
The first assumption is the incompressible hypothesis. In the incompressible flow, the pressure can be obtained from the mean velocity field through the Poisson equation, and the velocity field contains all the mean field information. However, most flow problems in actual engineering are compressible. The pressure becomes an independent state variable, which should be added to the arguments set.
The second assumption is the weak historical effect and homogeneous hypothesis. The consideration of including the historical effect and inhomogeneous turbulence is also a key direction in traditional modeling research. A commonly used approach is to consider the temporal and spatial historical effects in S to produce an "effective" strain rate S [43], which can be expressed as a convolution form: where Λm is the turbulent time scale. We perform series expansion at local (x, t) on the equation We analyze the right-hand side of the equation in sequence. (S 2 +Ω 2 ) can be expressed by the integrity basis of (S, Ω). νΔS represents the viscous diffusion and is not an active source term dominating the historical effect. The remaining two terms are the pressure gradient related term and the Reynolds stress gradient related term. If we want to represent the effect of S , the two gradients should be included.
Based on the analysis above, this paper adds two additional tensor arguments into the original (S, Ω), which are the pressure gradient vector vp and the turbulent kinetic energy 9 (marked as TKE below) gradient vector vk, defined as follows: (11) The final Reynolds stress isotropic tensor function is ( ) turb turb ;, , , , pk TL = τ fSΩ vv (12) Generally, the TKE is solved by its own transport equation and the expression above can be normalized to where b is the nondimensional Reynolds deviatoric tensor / 2 / 3 k =− b τ I , and the superscript ()  means normalization using turbulence scales, as equation (14) shows.
It is worth mentioning that the final tensor argument set is basically the same as in earlier studies [3,32], except including the density into the pressure term and the alternative normalization. The main purpose of this part is to systematically deduce the additional tensor arguments rather than determine the arguments ad hoc.

Complete irreducible tensor invariants and integrity basis
The implication of the tensor function representation is briefly introduced in Section I A.
Here, we restate the representation of b in a mathematical manner. Given a set of tensor arguments, the isotropic tensor function representation theorem indicates that any composed tensor function can be expressed as a linear combination of several tensor bases: (I1 ~ Ia) are the complete and irreducible tensor invariants formed by the argument set.
Completeness means that all the other invariants can be represented by these invariants.
Irreducibility means that they are independent of each other. Ti, i=1~w are the complete and irreducible tensor bases. Ti are all second-order symmetric tensors and collectively referred to 10 as the integrity basis. gi are the representation coefficients corresponding to the tensor bases, which are all functions of the invariants (I1 ~ Ia).
The representation process is to acquire the invariants and integrity basis. Theoretically this set of invariants and tensor bases is applicable to all symmetric tensors. If identifying b as the target, a corresponding set of gi is settled simultaneously.
There are two methods to compute the invariants and integrity basis. The traditional method uses the Cayley-Hamilton theorem. This theorem indicates that a high-degree tensor polynomial can be expressed by low-degree polynomials. In the actual computation process, the general form of the tensor polynomial composed of this set of tensor arguments is given first. Then, the C-H theorem is repeatedly applied to the polynomial to simplify the expression.
Finally, a set of low-degree tensor bases and invariants appearing in the reduction process are acquired. The deduction of Pope [35] followed this method. However, there are deficiencies.
First, the general form of tensor polynomials is easy to acquire only if the number of arguments is small, for example, only S and Ω. If the number increases, the general form will be complex and multiple, increasing the difficulty of reduction. Second, it is also difficult to prove whether the final results are irreducible. For example, the 10 tensor bases given by Pope were proven not to be the minimal representation in later research [40].
We employ the other method proposed by Zheng [44] which directly constructs the basis rather than simplifying it from a complex situation. This method can directly and systematically deduce the results for any number of tensor arguments. We suppose a symmetric tensor function H is composed of L symmetric tensors Ai, m antisymmetric tensors Wp, and N vectors vm: It is difficult to directly find invariants and tensor bases and to verify their completeness and irreducibility. Therefore For invariants, because of the definition of invariance, they can be computed using the components of all arguments no matter in which coordinate system and remain unchanged.
Therefore, in turn, if we construct a set of invariants that can represent all the components in a certain coordinate, these invariants are complete and can represent all the scalar-valued functions. The requirements above for the invariants (I1 ~ Ia) can be described mathematically as: where χt represents each component of tensor arguments. Because a symmetric tensor contains 6 independent components, an antisymmetric tensor contains 3, and a vector contains 3, there are in total (6L+3M+3N) components.
For tensor bases, after acquiring the invariants above, to further represent a symmetric tensor-valued function, we should select a set of tensor bases that can express all 6 complete orthogonal bases expanding the entire symmetric tensor space, which is described as:  (18) where the coefficients ηw are isotropic functions of (I1 ~ Ia) and E1 -E6 are:  (19) e1 -e3 are three independent orthogonal basis vectors constructing the entire space.
One advantage of the method is that the coordinate can be properly selected to minimize the number of components needing representation. For example, when dealing with a symmetric tensor, the coordinate can be the same as its principal axes, and there only exist 3 independent components. When multiple tensor arguments exist, all the possibilities between the principal axes need to be considered and the final complete form is acquired.
In the current research, we employ the method above and further develop the conclusion.
The original work only gave the general expression in three-dimensional space of the situation in which the numbers of Ai, Wp are arbitrary but all the vm are collinear. We deduce the situation with arbitrary numbers of Ai, Wp and vm. In addition, the expression in two-dimensional space is also acquired. The two expressions are listed in the Appendix. Here we only give the twodimensional situation of the argument set )ˆ(, , , pk S Ω vv . There are a total of 11 invariants and 7 tensor bases, as follows: nvariantsˆˆˆˆˆˆt rˆˆˆt The form of vectors vp / vk is worthy of mention. In earlier studies, the vectors are transformed to corresponding antisymmetric tensors using where ε is the permutation tensor. In the current research, the vectors are directly introduced. To verify which treatment is rational, we take kk = −  A ε v in a two-dimensional situation as an example. The components of b, S, and Ak are as follows: 11 where 13 It can be found that this term makes no contribution to the component to the deviatoric stress b. Therefore, it is physically unreasonable. As a comparison, if vk is used, the term ( ) kk  vv is: (24) where 12 This term makes an effective contribution to b. 13 In summary, in this part, we extend the tensor argument set from (S, Ω) to (S, Ω, vp, vk) and deduce the corresponding invariants and integrity basis in three-dimensional and twodimensional situations. The difference between the current treatment and earlier studies is analyzed, proving the validity.

B. Framework construction
In this part, we illustrate the construction of the data-driven turbulence modeling framework according to the sequence of preprocessing, training, predicting, and solving.

Training data preparation
The first step of preprocessing is to determine what data should be used for the construction.  The following issue after the flow case selection is to determine the status in which the training data should be. In Section I B, we propose the truth compatibility requirement, which is simply described as "all the features are constructed from the true value". The subsequent problem is that the true values of some quantities cannot be acquired from the high-fidelity database. To overcome this problem, we extend the implication of truth compatibility by substituting the "true value" with the "end-state value".
The end-state value is defined as the quantity for which the computation converges. To better illustrate the end-state of different variables, the CFD solving process is explained first, as shown in Fig. 3.  At the end of this section, the truth compatibility problem is further discussed and the 15 differences between the current research and the former researches [17,21,36] are clarified.
Computing the frozen turbulence model equations using true mean flow quantities is a commonly employed method in iterative frameworks to acquire the compatible turbulent quantities. The first difference is the treatment of TKE. In the work of Liu et al. [21], the k pDNS is the only TKE appearing in the framework, which is reasonable because the predicting target in the work is νt and the true TKE is not necessary. Schmelzer et al. [17] employed a more When predicting the entire Reynolds stress, the TKE is predicted using a sole machine learning model corresponding to the true TKE from DNS results. The machine learning model of TKE will be further illustrated in Section B 3. As a summary, the current method retains the original simple and robust form of the turbulence model equations and predicts the true TKE for the reconstruction of the Reynolds stress simultaneously.

Input feature selection
Once the training data are acquired, the next step is to construct the input feature set. In this paper we follow the feature selection criteria proposed in our earlier research [32]. The input features are constructed from two perspectives: tensor analysis and flow characteristic.
The tensor analysis perspective is based on the Reynolds representation illustrated above.
The acquired 11 complete irreducible tensor invariants of (S, Ω, vp, vk) are shown in equation (20). The number and some expressions of invariants in the current research are different from those in earlier research. This is because a new treatment of vectors is employed, and the new one is more effective as proven in Section II A. Therefore, the final input feature set with 9 features is established and listed in Table 1. It is worth mentioning that the invariants employ a different normalization method and are marked with the superscript ()  , shown as follows: where the symbol "||⸱||" represents the tensor norm. Such a "a/(|a|+b)" treatment can constrain the value range to (-1, 1) without significantly changing the original distribution. The thought is also employed in some flow characteristic perspective features.

Predicting target selection
We have illustrated the concept of the Reynolds stress representation in Section II A. The representation form is given in equation (15) and the expression of tensor bases is given in equation (21). Therefore, the only unknown variables are the representation coefficients gi. In traditional modeling, the stress expression is substituted into the algebraic stress equation and simplified. The acquired coefficients are complex polynomials of the invariants [40]. In datadriven modeling, with the help of the DNS database, the determination of gi becomes a numerical optimization problem.
The predicting targets in former ML frameworks based on the Reynolds stress representation generally follow the idea of a tensor basis neural network (marked as TBNN below) [13]. This concept is characterized by embedding the combination of gi and Ti into the ML model, as shown in Fig. 5 (a). The output of the ML model is the directly predicted deviatoric stress b pred , and the loss function is defined as The gi term is not explicitly shown, but only performs as the latent variables. No more preprocessing is needed besides separating the DNS Reynolds stress τ true into the magnitude and b true . 18 Several deficiencies exist in the framework above. First, embedding giTi into the ML model makes gi inaccessible. The training process only minimizes the discrepancy of the final combination, while the distribution of each coefficient is ignored. This mixes the error of model training and the representation. Second, the estimation of the TKE is generally missing in former studies. The reason for this might be that the TBNN using dimensionless features as input cannot be directly used to map a dimensional quantity.
To overcome these shortcomings, we developed an improved framework and corresponding prediction targets, as shown in Fig. 5 (b). The framework is separated into two parts, the ML part and the Reynolds stress representation part. The purpose of the representation part is to acquire the true coefficient set true i g and the TKE discrepancy Δln k. The true i g is determined by the pointwise least-square approximation, as follows: The TKE discrepancy is defined as the logarithm of the ratio between the true value and the pDNS value, as follows: After representation, the ML model is trained to model true i g and Δln k.
The current framework can remedy the problems of the former framework in the following aspects. First, the true TKE value is predicted, which further improves the modeling accuracy of the Reynolds stress. The treatment of Δln k makes the target dimensionless and suitable for the ML model. In addition, by introducing the pDNS value computed by the baseline turbulence model, prior knowledge is included in the model, which can also decrease the modeling difficulty and benefit the convergence. Second, the separation of the ML and the Reynolds stress representation enables us to evaluate the physical interpretation and smoothness of representation coefficients before the training. As will be stated in Section II C, the direct computation of

Model training and prediction
The model training and prediction flow chart is shown in Fig. 6, which can be summarized in the following steps: 1) Interpolate the DNS mean flow results onto the RANS grid and discretize to acquire q DNS , interpolate the DNS Reynolds stress onto the RANS grid to acquire τ DNS .
2) Iterate the turbulence equations with the mean flow quantities frozen as the DNS results to acquire (k,ω) pDNS .

4) Train the ML model:
5) During the application, the computation restarts from the baseline RANS results. In each iteration step, the ML model is executed using q and (k,ω) to predict the Reynolds stress τ pre .
6) The final mean flow q|final is acquired after the computation converges. It is found that the addition of the pressure gradient vp correlation terms does not significantly improve the representation effect and leads to divergence during the substitution computation. In addition, the high-degree tensor polynomials such as T5 -T7 are small-valued and exhibit strong numerical oscillation, which also has no effect on the representation process.
Therefore, 3 tensor bases are selected and remarked as T1 -T3: The components in the shear stress direction are shown in Fig. 7. a b 21 c Fig. 7 The components in the shear stress direction of the selected tensor bases. a 1 = TS .

Adaptive regularization method
The representation coefficients can be computed by solving the optimization problem given by Eqn. (27). However, it is found that such acquired coefficients lack smoothness and physical interpretability if the original expression is directly used. Taking the α = 1.0 case as an example, the coefficient g1 acquired by directly solving Eqn. (27) is shown in Fig. 8 (a), and the corresponding Reynolds shear stress is shown in Fig. 9 (a). The true value is shown in Fig.   9 (e). Although the representation effect is quite promising, the coefficient distribution is far from satisfactory. The presented value range is truncated for the following comparison, while the actual value range is far greater. Such severely varying coefficients will result in substantial prediction error during the ML model training.
In addition to the magnitude, the distribution of directly solved coefficients also faces nonphysical and unsmooth problems. The relation between g1 and eddy viscosity νt can be deduced by taking the first term in Eqn.
Note that whether Cμ is included in the definition of νt varies in different programs. In the current program, Cμ is included.
Therefore, a negative g1 corresponds to a positive eddy viscosity and positive dissipation. 22 However, the value of the directly solved g1 is positive in a considerable part of the flow field.
The CFD solving process diverges when substituting the coefficients into the RANS equations, which confirms the deterioration of robustness because of negative dissipation. The unsmoothness is marked by the dashed box in Fig. 8 (a), where the coefficient jump occurs.
In consideration of the coefficient distribution, it is found that the non-physical and unsmooth areas are mainly located at the main flow near the upper surface. The flow field in these areas has almost no mean flow characteristics such as the velocity gradient; therefore, there is almost no turbulence production. However, because of the spatial transport effect of the turbulence, the Reynolds stress can still be conveyed from other areas. If the coefficients are still computed using the original expression (Eqn. (27)) in these areas, the representation process can be analogous to "dividing by zero" and lead to the value jump. Although the Reynolds stress in these areas cannot be accurately represented, the actual magnitude is relatively small, and the influence on the mean flow is limited, which reminds us that the coefficients in these areas can be specifically processed to ensure smoothness.
To overcome the problems above, two methods are proposed in the current research. First, the value range of g1 is constrained to (-∞, 0) to ensure positive dissipation, which can be solved using the constrained least square method.
Second, we introduce the regularization term R. Regularization is a commonly employed method to limit the drastic changes of coefficients. By adding an additional term into the target function, the optimization not only minimizes the original target but also considers the effect of R. Taking the neural network training as an example, the modified loss function L with the commonly used L2 regularization term is generally defined as: where wi represents the model weights and λ represents the regularization parameter controlling the regularization intensity. The training process minimizes the original loss L and the magnitude of wi simultaneously, which can avoid overfitting due to large model coefficients.
In the current research, the purpose is to ensure smoothness and computational robustness.
As proven above, a small coefficient g1 does not fulfill the requirement. Therefore, physicsinformed regularization is introduced, ensuring that the computed representation coefficients do not deviate much from the baseline turbulence model result, which is represented by base i g .
The regularization term for each representation coefficient is defined as: Such a definition is compatible with different baseline models. If linear eddy viscosity models are employed, only base 1 g exists. If a nonlinear eddy viscosity model is employed, additional prior knowledge can also be included. The optimization problem with regularization is: To verify the regularization effect, large and small λ values are selected, and the g1 distributions are shown in Fig. 8 (b) and (c). The corresponding Reynolds shear stress distributions are shown in Fig. 9 (b) and (c). The comparison shows that employing regularization surely constrains the coefficient near the baseline value but also decreases the representation effect. λ = 0.1 leads to better smoothness, but the discrepancy between the true stress and represented stress is larger. In contrast, λ = 0.001 cannot eliminate the unsmooth area, but the represented stress is closer to the true value.
where λmin and λmax are the minimum and maximum values of the representation parameter, respectively, β is the multiplying factor function of the tensor basis norm ||Ti||, and the expression is: where θi is the threshold for different Ti, which is predetermined manually based on the distribution of ||Ti||. When ||Ti|| > θi, β approaches 0, and λ approaches λmin; conversely, when 24 ||Ti|| < θi, β approaches 1 and λ approaches λmax. In the current research, the θi for each Ti is [0.1, 0.05, 0.01]. λmin and λmax are 0.001 and 0.1, respectively. It is worthy mentioning that the selection of θi seems to be ad hoc. The entire strain rate and rotation rate magnitude will no doubt vary a lot for different flow conditions. But after normalized by the turbulent quantities, the nondimensional ||Ti|| actually doesn't vary so much. The final prediction results of the stress field and mean flow field at different geometries can prove the validity.
The computed coefficient distribution and corresponding Reynolds shear stress are shown in Fig. 8 (d) and Fig. 9

A. ML model training
The ML model is trained after the training data preparation and representation coefficient computation according to the methods above. The random forest (RF) is selected as the ML model in the current research. The RF model is a representative of integrated learning that is realized by combining multiple simple learning machines. During the training of RF, the entire training set is separated into several subsets and a decision tree is trained on each subset.
Random feature selection is introduced to ensure the "good but different" characteristics of the decision trees and avoid overfitting. The final prediction result is acquired by averaging all the decision tree results. The entire flow chart is shown in Fig. 11. The ML training is realized using the sci-kit learn package [48] in Python.

C. Prediction performance
The Reynolds stress and mean flow prediction are evaluated after the iterative computation converges. To verify the generalization capability, we mainly focus on the performance of the three testing sets. From the aspect of the hill slope, the α = 0.5 and 1.5 cases are extrapolation cases, and the α = 1.0 case is an interpolation case. The comparison between the predicted and theoretical values of Reynolds shear stress for the testing cases is shown in Fig. 14. The theoretical values are computed by the Reynolds stress representation using the true Reynolds stress and mean flow field. Therefore, the theoretical value is only the part that can be represented by three tensor bases but is close enough to the true stress, as compared in Fig. 9.
More details about the treatment of the discrepancy between the two distributions are discussed in Section IV.
The theoretical values indicate that the Reynolds stress of the testing cases are quite different in distribution and magnitude, which requires the ML model to have a strong Further evaluation of the mean flow field. The velocity contours are shown in Fig. 15, and the profiles are shown in Fig. 16. The comparison of lower surface friction is shown in Fig. 17.
In these comparisons, the reference values are selected as the true value directly from the DNS results since the final purpose of the data-driven turbulence modeling is to correct the mean flow field.
The velocity contour comparison indicates that the mean flow field has better smoothness than the Reynolds stress distribution. This is reasonable because the Reynolds stress can be

IV. Discussion
The main framework, key methods, and prediction performance are illustrated in the two sections above. Two more issues are discussed in this section. The first concerns the remaining Reynolds stress. Although the discrepancy between the representation stress and the true value stress is small, it still exists. Can the prediction result be further improved if the remaining part is also included in the targets? The second concerns the ML model calling frequency. Because the data are transferred between two different programs, the efficiency still needs further improvement. What if the ML model is not executed in each iteration step, but instead executed after an interval of several steps? The effects are evaluated below.

A. Effect of the remaining part after the Reynolds stress representation
The remaining part of the Reynolds stress, τb, is defined as the discrepancy between the represented Reynolds stress, τL, and the true value of Reynolds stress τ: 32 As τ and τL are second-order symmetric tensors, τb is also a symmetric tensor. Therefore, it can be eigen-decomposed: (38) where Λ = diag(λ1, λ2, λ3) and R = [v1, v2, v3] is the rotation matrix constructed by three eigenvectors. Note that 1 2 3 0    ++= because the TKE is included in τL., therefore, only two eigenvalues are independent.
To better utilize the represented result, the targets of eigenvalues are defined as the discrepancy between τL and τ: The original eigenvectors cannot be directly used as predicting targets because they are not spatially invariant. Therefore, the same process as used in the previous research [3,32] is needed, which is to compute the Euler angles describing the rotation from τL to τ. The details can be found in the mentioned literature and are not listed here. A schematic of the process is shown as: where only 3 z  is effective in the two-dimensional case.
In summary, if τb is considered, there are 3 more predicting targets: (1) Not considering τb: 1 2 3 {Δ ln , , , } k g g g (reference) (2) Considering only the eigenvalues: {Δ ln , , , , Δ , Δ ,} z k g g g    The computation of the three scenarios is performed for the α = 1.0 case. To eliminate the effect of the ML model prediction error, the targets above are all frozen as the theoretical values 33 during the computation. The mean flow velocity contours are shown in Fig. 18. The comparison indicates that the introduction of τb does not improve the prediction and even worsens smoothness and convergence. Therefore, the method of discarding τb employed in the current research is rational. Considering only the eigenvalues. c Considering the eigenvalues and eigenvectors.

B. Effect of the ML model calling frequency
As described in Fig. 3 and  Table 2. It can be observed that the iteration cost increases considerably, but the time cost relative to unsteady simulations such as LES or DNS is still acceptable.
A natural idea is to execute the ML model after an interval of several steps. To verify the effect, we modify the program and test four intervals: executing the ML model every 1, 3, 5, and 10 steps. The time costs are also listed in Table 2. The results show that the computation will diverge if the interval is too large. The ratio of total time cost between different intervals is essentially the same as the ratio of the interval steps, which indicates that the data transfer and Python computation occupy the main proportion of the time cost. The 3 converged mean flow results are generally the same, but the per 5 step case shows vibration near the periodic hill top, as shown in Fig. 19. 34 In summary, executing the ML model after an interval can accelerate the computation.
However, the smoothness and the convergence will be affected if the interval step is too large.  A. Three-dimensional situation