{'fraction_non_alphanumeric': 0.03043197446905513, 'fraction_numerical': 0.014893051175867178, 'mean_word_length': 4.723418134377038, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, '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}}
Data-driven turbulence modeling in separated flows considering physical mechanism analysis

Accurate simulation of turbulent flow with separation is an important but challenging problem. In this paper, a data-driven Reynolds-averaged turbulence modeling approach, field inversion and machine learning is implemented to modify the Spalart-Allmaras model separately on three cases, namely, the S809 airfoil, a periodic hill and the GLC305 airfoil with ice shape 944. Field inversion based on a discrete adjoint method is used to quantify the model-form uncertainty with limited experimental data. An artificial neural network is trained to predict the model corrections with local flow features to extract generalized modeling knowledge. Physical knowledge of the nonequilibrium turbulence in the separating shear layer is considered when setting the prior model uncertainty. The results show that the model corrections from the field inversion demonstrate strong consistency with the underlying physical mechanism of nonequilibrium turbulence. The quantity of interest from the observation data can be reproduced with relatively high accuracy by the augmented model. In addition, the validation in similar flow conditions shows a certain extent of generalization ability.


Introduction
The emergence and development of turbulent flow separation is crucial to the performance of a wing or turbomachinery. For example, trailing edge stall, which causes loss of lift coefficient, is induced by the breakup of the separation bubble appearing from the trailing edge [1]. In contrast, separation on an airfoil with horn ice usually begins at the leading edge at a small angle [2], which is characterized by the abrupt emergence of a separation bubble behind the ice shape and vast loss of the lift coefficient.
The numerical simulation of turbulent flow separation is essential for aerodynamic design applications. Although high-fidelity methods such as direct numerical simulation (DNS) and large eddy simulations (LES) can produce reliable results, the 3 computational cost is too high to apply to daily engineering design applications.
Efficient and reliable Reynolds-averaged Navier-Stokes (RANS) models will continue to be applied in daily engineering design in the foreseeable future [3]. The consensus is that the turbulence model is the most important factor to accurately evaluate the behavior of flow separation. Nevertheless, widely used RANS models such as the one- Studies show that the results of SA and SST models are far from satisfied in predicting the clean airfoil stall performance caused by boundary layer separation [4] or the iced airfoil stall performance induced by separating shear layer [2].
Previous studies have revealed that the nonequilibrium effect in the separating shear layer is one of the reasons for the failure of the present RANS models. Most turbulence models are calibrated by equilibrium turbulent shear flow, where the production-to-dissipation rate of turbulence kinetic energy is approximately 1.4 to 1.8 [5]. However, experimental and numerical studies have shown that the production of turbulence kinetic energy could be significantly larger than dissipation in a separating shear layer [6] or in the region of corner separation, as well as the wake and boundary layer of a compressor flow [7]. There have been some modifications for the turbulence models to include the nonequilibrium effect in a separated flow. In our previous studies, Li et al developed a separating shear layer fixed (SPF) − 2 ̅̅̅ − model by applying 4 a multiplier to the destruction term of the equation to improve the accuracy of stall prediction for an iced airfoil [8][9] [10].
However, ad hoc modifications for turbulence modeling usually have poor generality on different flow problems. For example, Rumsey [11] found that although the enhancement of turbulent mixing in a separated shear layer improves the prediction accuracy for a periodic hill, the same corrections obtained worse results for backstep flows. This indicates that the nonequilibrium turbulence behavior might be different in different flow cases, which is difficult to precisely describe by manual functional corrections. On the other hand, the previous correction of Li et al [8] merely increased the production-to-dissipation rate in the shear layer region, while the LES results by Fang et al [7] showed that the production can be either increased or decreased in such regions. Consequently, a more rigorous process of quantifying and reconstructing model corrections, especially for the nonequilibrium effect of the separated flow, is needed.
In recent years, data-driven approaches have been increasingly applied to turbulence modeling. This benefits from the increasing data from experiments or highfidelity simulations and the application of statistical and machine-learning methods in turbulence modeling. From the aspect of data-driven turbulence modeling, the simplifications for Reynolds stress closure bring uncertainty to the RANS simulation.
The uncertainty propagates forward through RANS computations and ultimately leads 5 to uncertainty in the predictions. Data-driven methods can be used for quantifying or reducing the uncertainties of the RANS models [12] and can be further used for the augmentation of RANS models.
There are typically two types of studies focused on the uncertain model parameter θ and its prior probability distribution p(θ) . From the forward perspective, the probability distribution of QoI can be determined by sampling from p(θ) and computing QoIs for each sample [12] [13]. From the backward perspective, when given observable data D from DNS or experiments, the posterior probability distribution of model parameters p(θ|D) can be inferenced by the Bayes formula:  [17]. Dow and Wang used the adjoint method 6 to infer the viscosity field in a plain channel flow to fit the velocity field obtained from DNS [17]. Singh and Duraisamy studied the flow around an airfoil at a high angle of attack, where a strong adverse pressure gradient and flow separation existed [18]. They modified the SA model by multiplying the production term by a factor (  [20]. Apart from steady flow problems, data assimilation has also been used for unsteady problems. Meldi and Poux proposed the sequential data assimilation method using Kalman filtering [21], and He et al. improved the method using the continuous adjoint method [22]. Field inversion is usually converted to an optimization problem based on the adjoint method [23]. However, there may be innumerable local optimal solutions for the optimization problem. Many of them are physically unfeasible and may be obtained by gradient-based optimization algorithms. It is a critical problem to consider the physical mechanism in the field inversion process and locate the model uncertainties with physical realizability. Setting the prior based on physical intuitions has been proven to be an effective approach for this problem [16] [20]. Machine learning (ML) is another frequently used approach for data-driven turbulence modeling. Many studies have used ML models to assist RANS modeling by predicting the source terms in model equations [18], eddy viscosity [24] [29]. They used the wall isentropic Mach number distribution as the experimental data for correction, and they limited the range of correction to improve the robustness of the correction factor.
Previous studies have established FIML as an effective framework for the correction of RANS models. However, how to consider the physical mechanism in field inversion is a problem that has drawn little attention. For example, nonequilibrium turbulence is an important physical mechanism of a separated shear layer. Constraining the correction 8 area and limiting the correction range in the field inversion based on the nonequilibrium turbulence phenomenon might be an effective method to obtain a physically feasible modification of turbulence modeling.
In this paper, FIML is used for data-driven turbulence modeling considering the nonequilibrium effect. The SA model is modified by multiplying the production term of ̃ with a spatially distributed correction factor ( ), which is closely related to the equilibrium state of the turbulence. The corrections are applied to several problems, namely, the S809 airfoil, a periodic hill, and the GLC305 airfoil with ice shape 944, to Finally, the augmented model is examined with flow conditions similar to those of the datasets, and a certain generalization capability is demonstrated.

Framework of field inversion and machine learning
The framework used in this paper is similar to that in references [28] [29], as shown 9 in Fig. 1. First, field inversion is implemented in different flow cases to obtain corresponding model corrections of the SA model. Then, the dataset from field inversion is used to train an ANN. Finally, the ANN is integrated into the original SA model to form an ML augmented model. Fig. 1 Framework of FIML used in this paper

Field inversion
The SA model [30] contains a transport equation for eddy viscosity ̃, which is written as: Here, is the molecular viscosity. = 2/3, 2 is a model parameter, and ̃ is the solution variable, which is used to compute the eddy viscosity : is the production term, and is the destruction term: where is the strain rate tensor of the mean velocity field, d is the minimum distance to the wall, = 0.41, and 1 , 1 , 2 , 3 , 1 , 3 and 4 are model constants.
Although the SA model has been validated in many flow cases, there are several flow problems that the model has poor accuracy, such as the massive flow separation problem. Many efforts have been made to improve the SA model for flow separation problems by implementing modifications to the model terms, such as SAE [31] and SA-QCR [32]. In this paper, a data-driven augmentation approach is applied to the SA model by enhancing its model parameters. The production term of ̃ is multiplied by a spatially distributed factor , making the transport equation: Statistically speaking, rather than physical analysis, ( ) might be a random field. The prior probability distribution p( ) can be determined manually with prior knowledge of turbulence modeling before any observation data from experiments or 11 DNS is given. Once given observation data , the posterior probability distribution of can be obtained by Bayesian inference. We assume that and are Gaussian and independent of each other on different grid points while having the same variance.
Furthermore, supposing a prior of 1.0 for in each grid point, the posterior of can be written as: where 2 and 2 are the variance for and , respectively, and ℎ ( ) is the predicted QoI at the i-th observed point with a given . In this paper, MAP inference is used to find the optimal correction , which is equivalent to an optimization problem with an object function . The first part of represents the prediction error with a given , and the second part of represents the deviation from the prior of , which can also be interpreted as regularization. The weight between these two terms is adjusted by tuning and .

Discrete adjoint method
The adjoint method is an effective approach to compute the derivative of an object function with respect to a large set of variables when involving partial differentiation equations (PDEs). This method is often combined with gradient-based optimization algorithms to solve optimization problems [33] [34]. Due to its low computational cost to the number of variables, the method is particularly suitable for the optimization problem in field inversion, where the parameter on each grid point might be a design 12 variable. The discrete adjoint method [35][36] is used in this paper. The algorithm is briefly introduced as follows.
The aerodynamic object function and the residual of RANS equations can be expressed as functions of the coordinate and the flow variables at each grid point: ( , )， ( , ). Given the constraint of RANS equations ( , ) = , the derivative of with respect to the design variables can be expressed as follows: If we take the differential of equation ( , ) = , we obtain: Thus: We introduce an adjoint variable , which satisfies: Then, we obtain: The discrete adjoint computation involves first obtaining and , solving the discrete adjoint equation to obtain , then obtaining and , and finally obtaining according to Eq. (11). The RANS solver and the discrete adjoint solver in this paper are based on the open-source CFD code ADFlow [37] with secondary development.
The automatic differentiation code Tapanade [38] is used to generate the 13 forward/reverse mode differentiation code for sensitivity evaluations.
In this paper, an augmented discrete adjoint method [39] is applied to handle the constraint of the aerodynamic optimization problem, such as the fixed lift coefficient of airfoil computation or specified mass flow rate of the internal flow.

Machine learning
Field inversion provides an optimal correction field for each flow condition of the training datasets. This case-specific message can be generalized by establishing the functional relationship between the correction and some local physics features. The choice of the features is crucial for the effect of the augmentation model and is usually dependent on physics intuition. The features used in this paper rely on [28] and [29]: where =̃/ is the eddy viscosity ratio, and are respectively the production and destruction terms of the transport equation of ̃, and Ω are respectively the magnitudes of the strain rate and vorticity, is the magnitude of the Reynolds stress, and |∇̃| is the magnitude of the gradient of ̃.
The parameter ′ is proposed by [14] as a modification of : 3 ) was initially proposed for detached eddy simulations [39] and was used as a feature for ML in reference [28].
Various ML algorithms have been used to assist RANS modeling. Parish et al. and 14 Duraisamy et al. used the Gaussian process as the ML algorithm for FIML [23][27].
One advantage of Gaussian process regression is that it can predict both the mean and the variance of the output, so the uncertainty of the correction can be determined simultaneously. Wang et al. used a random forest to reconstruct the Reynolds stress model discrepancy [40]. Schmelzer used sparse symbolic regression to learn Reynolds stress expression [41], and Zhao et al used the gene-expression programming (GEP) method to train a symbolic turbulence model [42]. Singh et al. [28] and Ferrero et al. [29] used an ANN for FIML. Among those ML algorithms, the evaluation cost of ANNs is independent of the size of the training data. Because the ML model is called at each grid point and each timestep, the computational cost is an important factor for integration into an RANS solver. Consequently, an ANN is used as the ML algorithm in this paper. The structure of the ANN is shown in Fig. 2. Note that the specific structure of the ANN is chosen according to not only the error for regression but also the error for reproducing the experimental data when integrated into the SA model, and there might be no certain correlation between these two characteristics. This is discussed in the following section. 15 Fig. 2 Structure of the ANN used in this paper

S809 airfoil stall case
The S809 airfoil is designed for a wind-turbine blade with 21% relative thickness.
The experimental data were obtained in the low-turbulence wind tunnel of the Delft University of Technology Low Speed Laboratory [43]. The flow around the S809 airfoil near the stall shows a strong adverse pressure gradient (APG), which causes separation at the tailing edge. Previous studies indicated that the SA model tends to overestimate the lift coefficient of the S809 airfoil near the stall [28]. The nonequilibrium turbulence in the shear layer might be one reason for the discrepancy [28]. In this paper, FIML is used to increase the accuracy of the SA model. Note that this case has been applied by many researchers [28] [44]. It is also used as a validation case of the present FIML framework.
The computation grid is shown in Fig. 3. It has 457 points along the 16 circumferential direction and 97 points in the normal direction. The height of the first grid layer is less than 1 × 10 −6 c (c is the chord length), which ensures Δ + < 1.
where , is the variance of the prior of at the -th grid point, which creates a nonuniform variance field. For this case, , is set to be: This set of , comes from the physical knowledge that the uncertainty is mostly attributed to the strong APG near the boundary layer. A uniform variance field with , = 20 is also used for field inversion for comparison. The convergence of and 17 the object function at α = 14° at the two variance distributions are shown in Fig. 4.
The error of the predicted is reduced sufficiently regardless of whether it has uniform variance or nonuniform variance. Fig. 5 shows the correction fields for the two variance distributions. The production of turbulence viscosity is suppressed mainly in the region with a strong APG near the boundary layer. On the other hand, the difference between the two correction fields indicates that there could be multiple solutions for the correction field, all of which can reduce the error of QoIs sufficiently. The introduction of physical knowledge by setting the prior can help exclude solutions that are inconsistent with physical expectations. Consequently, the nonuniform variance prior is used in the following to consider the physical mechanism of field inversion.  Fig. 8. The open-source tool PyTorch [45] is used to train the ANN. The loss function is the mean square error (MSE), the training algorithm is Adam, and the activation function is sigmoid. Z-score normalization is used for the input features, and batch normalization is implemented between the hidden 20 layers. (1)̃ (3) / (4) Fig. 8 Features from the field inversion at α = 14° The ANN has three hidden layers. Four ANNs with different sizes of hidden layers, (8,4,2), (8,8,8), (16,16,8) and (64,32,16), are tested with different training epochs.
The root-mean-square errors on the training set and the validation set and the prediction errors of for the enhanced model are shown in Table 1. The reduction in the loss function during training does not necessarily lead to the reduction in the prediction error for the enhanced SA model. Specifically, when the size of the hidden layer increases, the prediction error can be even larger while the loss function is still decreasing, which is a kind of overfitting. In fact, when the ANN is (1) α = 8° (2) α = 14° Reynolds numbers (1 × 10 6 , 2 × 10 6 , and 3 × 10 6 ) are tested with a wide range of angles of attack. and computed by the original and enhanced SA model are 23 compared with the experimental data in Fig. 10. The enhanced SA model shows considerable improvement not only in the two cases of FI but also in a wide range of Reynolds numbers and angles of attack. This indicates that the ML process has extracted generalizable modeling knowledge from the specific FI cases. Moreover, although the drag coefficient is not provided for the FI, the enhanced model could still make better predictions on . This shows the powerful generalization ability of the FIML framework. The S809 airfoil case demonstrates that the present FIML framework is reliable.
(1) = 1 × 10 6 (2) = 2 × 10 6 (3) = 3 × 10 6  24 The periodic hill is a typical flow separation problem that is widely used for evaluating the performance of turbulence models. The large flow separation region makes this phenomenon challenging to simulate accurately. Extensive studies have focused on the correction of RANS models for periodic hill flows, using either manual turbulence modeling [11] or data-driven approaches [46] [47]. It is known from physical analysis that the model uncertainty in this problem mostly comes from the nonequilibrium turbulence effect in the separation shear layer [11], where the turbulence shear stress can be underestimated. To fix the failure of turbulence models in separation regions, Rumsey [11] modified Menter's k − ω SST model. As shown in Fig. 11, they increased the destruction term of ω by a multiplier in the separating shear layer. The fixed region is characterized by a nonequilibrium turbulence region, where the production-to-dissipation rate of turbulence kinetic energy is significantly larger than 1. With the same physics intuition, FIML is implemented on the periodic hill case to correct the original SA model in this paper.

Periodic hill case
(1) / (2) Fig. 11 Nonequilibrium turbulence effect and the corresponding modification in the separation shear layer of a periodic hill from Rumsey [11] 25 The geometries of the periodic hills used by Xiao et al. [48]  The constraint-augmented adjoint method is used for the sensitivity computation.
The bulk flow velocity is used as the aerodynamic constraint, and the volume force in the field can be adjusted automatically. layer where the nonequilibrium turbulence effect is significant according to Rumsey [11]. The rectangular region is shown as the red dashed box in Fig. 13, where λ = 8 × 10 −6 is set inside the region. The third case sets = 8 × 10 −6 in the whole computation domain, which assumes a uniform variance. The setting for the three cases are listed in Table 2. The convergence histories of the squared error of velocity in 27 different cases are shown in Fig. 14. The reduction in error is more difficult with a nonzero than with = 0.   Fig. 16.
Although the friction coefficient on the lower surface is not explicitly included in the objective function, the FI results are also closer to the DNS results at most locations. 28 (1) Overall view of the velocity profiles  Fig. 18.
The correction field of the first case is very irregular, and no general tendency can be found from the field. A local optimized solution is obtained by the adjoint solver.
However, the solution might be not physical. The second case constrains the corrections inside the prescribed region of the separation shear layer, which is marked as the dashed box in Fig. 18(2). The error of QoIs is still sufficiently reduced, with the production of  The correction field is also compared with the separation shear layer fixed region of the SPF k − 2 ̅̅̅ − model [8]. The SPF k − 2 ̅̅̅ − model makes use of the ratio of local turbulence production to dissipation to trigger an increase in nonequilibrium turbulence at the separation shear layer, and such modification can be represented by the term f NE [8]. Fig. 18(4) shows the resulting contour of f NE from the SPF k − 2 ̅̅̅ − model. This shows that the modification of the nonequilibrium turbulence is triggered right in the separated shear region. For the FI, the corrected area is almost in 31 the same region, as shown in Fig. 18(2) and Fig. 18(3), which further confirms the rationality of the physical mechanism of the flow FI. In this paper, FIML is implemented on airfoil GLC305 with ice shape 944 to modify the SA model. The experimental data used for FI were obtained from the NASA John H. Glenn Research Center at Lewis Field [49].
The geometry of the iced airfoil and the computation grid are shown in Fig. 23.
There are 457 and 97 points along the streamwise and normal directions, respectively.
The freestream condition is = 0.12, = 3.5 × 10 6 . Previous studies have shown that the SA model tends to underestimate and at relatively large angles of 36 attack [8], which needs to be corrected by the FIML approach in this paper.
The second case uses only the lift coefficient as the observation data, which is referred to as "FI2" in the following section. As shown in Eq. (16), the prior of is  The flow fields before and after the FIs are shown in Fig. 25. The correction for the flow field is relatively small at an 4° angle of attack. However, at an angle of attack of 6°, the separation bubble predicted by the original SA model covers the whole airfoil, while the corrected bubbles take approximately half of the chord length. The flow reattachment after the shortened separation bubble helps to increase the lift coefficient , which is closer to the experimental results in Table 3. The distribution of the correction is also compared with the SPF k − 2 ̅̅̅ − model [8], which made accurate predictions for flow around iced airfoils. As shown in This demonstrates the physical realizability of the correction by the present FI.
Moreover, both a decrease and an increase in turbulence mixing are obtained by FI, which shows the capability of more accurate inference for the specific flow case. 40 (1) /

Conclusions
Flow separation poses challenges to the precise modeling of turbulence. Previous studies have indicated that the model-form uncertainty in current RANS models mostly comes from the nonequilibrium effects related to the separation region. In this paper, the data-driven RANS modeling framework FIML is applied to the SA model to quantify and reconstruct the corrections related to nonequilibrium effects. The FIML framework is applied on three flow cases, namely, the S809 airfoil, the periodic hill, and the GLC305 airfoil with ice shape 944.
(1) FI is used to quantify the model-form uncertainty, with the prior of model uncertainty set manually using physical knowledge. The given observation data can be reproduced with significantly high accuracy after FI. In addition, the accuracy of the prediction of other QoIs can also be promoted even though they did not explicitly appear in the object function. Moreover, comparative analysis demonstrates that the correction made by the FI on iced airfoil and periodic hills accounts for the 46 nonequilibrium effect in the separation shear layer, which is also a crucial consideration for the SPF k − 2 ̅̅̅ − model.
(2) The ML model is used to predict the model correction given local flow features.
With the ANN trained with elaborately chosen hyperparameters, the machine-learning augmented SA model can reproduce the results from FIs with relatively high accuracy, and a certain extent of generalizability can be guaranteed at similar flow conditions.
Nevertheless, the generalizability of the augmented turbulence model might be further promoted in future research.