Analysis on numerical stability and convergence of RANS turbulence models from the perspective of coupling modes

Reynolds-averaged Navier-Stokes simulations are still the main method to study complex flows in engineering. However, traditional turbulence models cannot accurately predict flow fields with separations. In such situation, machine learning methods provide an effective way to build new data-driven turbulence closure models. Nevertheless, a bottleneck that the data-driven turbulence models encounter is how to ensure the stability and convergence of the RANS equations in posterior iterations. This paper studies the effects of different coupling modes on the convergence and stability between the RANS equations and turbulence models. Numerical results demonstrate that the frozen coupling mode, commonly used in machine learning turbulence models, may lead to divergence and instability in posterior iterations; while the mutual coupling mode can maintain good convergence and stability in the process of iterations. This research can provide a new perspective to the coupling mode for machine learning turbulence models with RANS equations in posterior iterations.


Introduction
Turbulence is considered as the last major unsolved problem in classical physics [1] . Numerical simulation methods play an increasingly important role in the research of turbulent problems. However, direct numerical simulation (DNS) and large eddy simulation (LES) are expensive in computational storage and cost for high Reynolds number flows. Therefore, Reynolds-Averaged Navier-Stokes (RANS) simulation is still the dominant tool for industrial problems in the near future.
Traditional RANS models, especially based on the eddy viscosity models (SA model [2] , k-ε/k-ω model [3]~ [5] , Menter's SST model [6] ), have been widely used to solve turbulence problems in engineering. These models can generally predict reliable aerodynamic coefficients for attached flows around complex configurations. However, the key difficulty that traditional RANS models encounter is that they cannot accurately simulate turbulent flows with separations.
Consequently, in recent years, researchers have developed many data-driven turbulence closures, including tensor based neural network (TBNN) models [8]~ [9] , machine learning augmented SA models [10]~ [13] , Reynolds stress discrepancy models [14]~ [16] and physics informed neural network (PINN) models [17]~ [18] . The models are trained by machine learning techniques based on DNS or LES data, which can significantly improve the predicted accuracy for separated flows compared with traditional RANS models. Although the recently developed models have made remarkable achievements in turbulence simulation, a serious problem is how to ensure the numerical stability and convergence when a machine learning turbulence model is coupled with the RANS equations in posterior iterations. The TBNN model, constructed by Ling et al. [8] , was able to reach convergence only based on the converged solution of traditional RANS models, and the authors only pointed out that future studies would investigate the feasibility of iterative convergence which the TBNN model was coupled at every iteration. Thompson et al. [19] found that in the channel flows directly substituting Reynolds stresses from DNS databases into RANS equations led to mean velocity fields with large numerical errors, and the errors would increase with the Reynolds numbers. Poroseva et al. [20] also confirmed that unphysical solutions would occur when the fixed Reynolds stresses were directly inserted to the RANS equations. Wu et al. [21] demonstrated that treating Reynolds stresses as an explicit source term in the RANS equations may result in the ill-conditioned problem, and they also considered that the implicit treatment of the Reynolds stresses (decompose Reynolds stresses into the eddy viscosity part and the nonlinear part) could effectively alleviate the convergence problem. Cruz et al. [22] proposed to replace the Reynolds stress tensor with its divergence, the Reynolds force vector, as the machine learning target, which can more accurately reconstruct the average flow fields from DNS data. More recently, Brener et al. [23] elaborately explored different methods to solve the ill-conditioned problem of the RANS equations according to the plane channel, the square duct and the periodic hill geometries, and results indicated that the RANS equations were still ill conditioned in the whole range of cases even when the Reynolds stress tensor was treated implicitly. They also recommended that the Reynolds force vector with implicit treatment using optimized eddy viscosity could mitigate the error propagation to the mean velocity field in RANS equations.
Guo et al. [24] considered that both the model error and the propagation error existed in RANS simulations, and the authors also mathematically derived the propagation error and claimed that it was an inevitable factor affecting the convergence of RANS simulations.
Mean velocity fields will produce large numerical errors when the Reynolds stresses extracted from DNS data with small errors are fixedly injected to the RANS equations, which is a key issue that the machine learning turbulence models are encountering. Many researchers proposed different strategies to circumvent the problem. Most existing work adopted the idea that inserting the fixed Reynolds stresses directly into the RANS equations, which did not consider the coupling effect between the Reynolds stresses and the RANS equations. The commonly used approach is to make the RANS solver firstly converge to the steady state with traditional turbulence models, and re-converge the simulation using the high-fidelity Reynolds stresses. However, this cannot ensure the iterative convergence if the simulation is started from uniform flow-fields. From a new perspective of the coupling mode between turbulence models and the RANS equations, this paper compares the convergence properties with different types of coupling modes, and we find out that the frozen coupling mode may result in numerical divergence and instability for the RANS simulation at high Reynolds numbers.
The remainder of this paper is organized as follows. Section 2 introduces the CFD governing equations and the methodology of RANS coupling modes. Section 3 presents the numerical examples to verify the proposed method, and finally the conclusions are drawn in Section 4.

Methodology
The steady-state RANS equations for incompressible flows can be written as: where u and p are the mean flow velocity and pressure;  is molecular viscosity; and ij R is the Reynolds stress tensor, which needs to be modeled. The term ij R can be decomposed into a linear part and a nonlinear part based on the eddy viscosity hypothesis, which is defined as: where t  denotes the turbulent eddy viscosity; is the velocity strain rate tensor; and ij R  is the nonlinear part of the Reynolds stresses. The momentum equation can be rewritten as: The main purpose to construct turbulence models is to resolve the RANS equations by modeling the unresolved term t  or ij R , which are injected into RANS equations in order to solve the mean velocity field. There are mainly two parts to We discover that the mutual coupling mode can ensure the numerical convergence and the stability in the process of iteration from uniform flows, while the frozen coupling mode may induce divergence and instability in the simulation of separated flows at high Reynolds numbers. This paper mainly investigates the influence of two coupling modes on numerical convergence and stability when the RANS solver is initiated from uniform flows.
According to the Boussinesq hypothesis of linear eddy viscosity, the nonlinear part of the Reynolds stresses is ignored, and Equation (3) can be simplified as first-order closure models: where M  denotes the modeling error introduced by the approximation of turbulence models. In the process of solving RANS equations by numerical methods, the propagation error P  will also occur, which can be expressed as follows: 6 The modeling error M  has no influence on the convergence and the stability since it is a statistic error if a turbulence model is determined. And the propagation error P  is produced during solving the RANS equations. P  approaches to zero if the solver converges; otherwise, the error will continue to increase with the iteration of the solver and eventually lead to numerical instability or divergence.
This paper mainly studies the influence on numerical stability and convergence of the propagation error P  with different coupling modes, and the model error M  is not taken into consideration currently. For the mutual coupling mode, the propagation error is denoted by MC P  , and Equation (5) can be rewritten as: For the frozen coupling mode, we can firstly solve the RANS equations and obtain the turbulent eddy viscosity in the convergent state. And then, re-compute the RANS solver by fixedly inserting the converged turbulent eddy viscosity. The propagation error is denoted by FC P  , that is * ( ) The RANS calculations are conducted with our in-house code which is based on the cell-centered finite volume method. In this paper, Roe scheme [29] is used to evaluate convective fluxes, and interface values are reconstructed by the second-order least-squares approach. The traditional SA model is used to solve the turbulent eddy viscosity. For the semi-discretized form of the governing equations, the implicit Symmetric Gauss-Seidel scheme [30] is adopted to march in time integral, and the local time stepping and residual smoothing techniques are also employed to accelerate convergence. Our in-house RANS solver has been validated completely by simulating many complex configurations in engineering at high Reynolds numbers, and the details can refer to our previous work [31]~ [33] .

Curved backward facing step
The first test case is the curved backward facing step, and the definition of the geometry is the same as Ref. [34]. The step height is H and the Mach number is set as 0.2, and the Reynolds number is 13,700 based on H and inflow velocity. The computational mesh near the step is displayed in Fig. 2, which is locally refined. The total number of elements is 24,849, and the mesh growth rate is set as 1.

Flow around a S809 airfoil at high Reynolds number
The second test case is the flow around a S809 airfoil at high Reynolds number.
The computational mesh is shown in Fig. 5 We further investigate whether the frozen coupling mode can ensure the numerical stability when the RANS solver is based on the converged flow. Firstly, the mutual coupling mode is adopted to iterate 30,000 steps and the turbulent eddy viscosity is obtained in the convergent state. And then, the RANS solver is re-computed by the frozen coupling mode to a large number of iterative steps.
Comparisons of convergence histories and lift coefficients for different coupling modes are shown in Fig. 11 and Fig. 12. With the increase of iterative steps, the mutual coupling mode can consistently maintain numerical stability and convergence.
However, for the frozen coupling mode, the propagation error of the RANS equations will be gradually enlarged with the process of iterations, and both the residual and lift coefficients will eventually present numerical oscillations and non-convergence. Fig.   13 depicts the comparisons of streamlines near the trailing edge at the end of iterative step, which clearly shows that the vortex in the separated region is unsteady and constantly shaking. The results demonstrate that the frozen coupling mode may induce numerical instability even though the RANS solver is based on converged solutions.
Step As for the frozen coupling mode, the turbulent eddy viscosity is fixed, and there is no information feedback during the whole process of numerical calculation, which cannot eliminate propagation errors effectively. As the iterative step continues to increase, propagation errors may gradually accumulate and are enlarged, which will eventually lead to numerical non-convergence and instability. On the contrary, the mutual coupling mode can dynamically adjust the balance between the unclosed term

Conclusions
Reynolds stresses from DNS data are injected to the RANS solver as closure terms, which can cause significant errors for the mean velocity flow-fields. This paper analyzes this problem from a new perspective of the coupling modes. Based on the RANS equations with the traditional SA model, we compare the influences of the frozen coupling mode and the mutual coupling mode on numerical convergence and stability. As a result, the widely used frozen coupling mode may induce divergence and instability for separated flows at high Reynolds numbers even based on converged mean flow-fields; while the mutual coupling mode can dynamically adjust the balance between unclosed terms and mean flows, which can effectively eliminate propagation errors during iterative steps and ensure numerical convergence and stability. The curved back facing step and the high Reynolds flow around the airfoil are used to verify the proposed idea. In conclusion, this paper demonstrates the problem from a novel perspective that the RANS equations cannot obtain accurate mean flow-fields by injecting high fidelity Reynolds stresses. The proposed viewpoint can provide a new breakthrough to solve convergence and stability problems currently encountered in the machine learning modeling.