Hostname: page-component-6898cfcb55-kwdml Total loading time: 0 Render date: 2026-11-11T14:35:15.868Z Has data issue: false hasContentIssue false
  • English
  • Français

Reynolds-averaged Navier–Stokes equations with explicit data-driven Reynolds stress closure can be ill-conditioned

Published online by Cambridge University Press:  29 April 2019

Jinlong Wu,
Heng XiaoOpen the ORCID record for Heng Xiao[Opens in a new window],
Rui Sun and
Qiqi Wang
Jinlong Wu
Affiliation:
Kevin T. Crofton Department of Aerospace and Ocean Engineering, Virginia Tech, Blacksburg, VA 24060, USA
Heng Xiao*
Affiliation:
Kevin T. Crofton Department of Aerospace and Ocean Engineering, Virginia Tech, Blacksburg, VA 24060, USA
Rui Sun
Affiliation:
Kevin T. Crofton Department of Aerospace and Ocean Engineering, Virginia Tech, Blacksburg, VA 24060, USA
Qiqi Wang
Affiliation:
Department of Aeronautics and Astronautics, MIT, Cambridge, MA 02139, USA
*
†Email address for correspondence: hengxiao@vt.edu
Rights & Permissions [Opens in a new window]

Abstract

Reynolds-averaged Navier–Stokes (RANS) simulations with turbulence closure models continue to play important roles in industrial flow simulations. However, the commonly used linear eddy-viscosity models are intrinsically unable to handle flows with non-equilibrium turbulence (e.g. flows with massive separation). Reynolds stress models, on the other hand, are plagued by their lack of robustness. Recent studies in plane channel flows found that even substituting Reynolds stresses with errors below 0.5 % from direct numerical simulation databases into RANS equations leads to velocities with large errors (up to 35 %). While such an observation may have only marginal relevance to traditional Reynolds stress models, it is disturbing for the recently emerging data-driven models that treat the Reynolds stress as an explicit source term in the RANS equations, as it suggests that the RANS equations with such models can be ill-conditioned. So far, a rigorous analysis of the condition of such models is still lacking. As such, in this work we propose a metric based on local condition number function for a priori evaluation of the conditioning of the RANS equations. We further show that the ill-conditioning cannot be explained by the global matrix condition number of the discretized RANS equations. Comprehensive numerical tests are performed on turbulent channel flows at various Reynolds numbers and additionally on two complex flows, i.e. flow over periodic hills, and flow in a square duct. Results suggest that the proposed metric can adequately explain observations in previous studies, i.e. deteriorated model conditioning with increasing Reynolds number and better conditioning of the implicit treatment of the Reynolds stress compared to the explicit treatment. This metric can play critical roles in the future development of data-driven turbulence models by enforcing the conditioning as a requirement on these models.

JFM classification

Mathematical Foundations: Computational methodsTurbulent Flows: Turbulence modellingTurbulent Flows: Turbulence simulation

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 869 , 25 June 2019, pp. 553 - 586
Copyright
© 2019 Cambridge University Press 

1 Introduction

Reynolds-averaged Navier–Stokes (RANS) simulations play an important role in industrial simulations of turbulent flows. The state-of-the-art eddy-viscosity models, e.g. the $k$ – $\unicode[STIX]{x1D700}$ , $k$ – $\unicode[STIX]{x1D714}$ and Spalart–Allmaras (SA) models (Launder & Sharma Reference Launder and Sharma1974; Wilcox Reference Wilcox1988; Spalart & Allmaras Reference Spalart and Allmaras1994; Menter Reference Menter1994), are based on the assumption that the turbulence production and dissipation are in equilibrium, and thus they perform poorly in flows with non-equilibrium turbulence (Speziale & Xu Reference Speziale and Xu1996; Hamlington & Dahm Reference Hamlington and Dahm2008; Hamlington & Ihme Reference Hamlington and Ihme2014), e.g. flows with massive separations or abrupt mean flow changes. On the other hand, Reynolds stress models (RSMs), also referred to as second-moment closures, take into account the transport of Reynolds stresses and thus have better performance than eddy-viscosity models for flows with non-equilibrium turbulence (Pope Reference Pope2000). NASA’s CFD Vision 2030 Study white paper identified advanced turbulence modelling based on Reynolds stress models as a priority for aeronautic technological advancement in the coming decades (Slotnick et al. Reference Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Mavriplis2014). However, so far Reynolds stress models have not been widely used in engineering applications despite their theoretical superiority. A key shortcoming of Reynolds stress models is their lack of stability and robustness, i.e. the difficulty in achieving convergence (Basara & Jakirlic Reference Basara and Jakirlic2003; Maduta & Jakirlic Reference Maduta and Jakirlic2017) and the high sensitivity to the modelling of unclosed terms (particularly the pressure–rate-of-strain tensor) in the Reynolds stress transport equations.

Among the leading causes of numerical instability of Reynolds stress models is the intricate coupling as dictated by the Reynolds stress transport equations (Pope Reference Pope2000, chap. 7). First, the shear stress $\unicode[STIX]{x1D70F}_{xy}$ is responsible for the production of streamwise fluctuations  $\unicode[STIX]{x1D70F}_{xx}$ , where $x$ and $y$ denote streamwise and wall-normal coordinates, respectively. Then, the turbulent kinetic energy in $\unicode[STIX]{x1D70F}_{xx}$ is redistributed by the pressure–rate-of-strain to the other two normal components of the Reynolds stress tensor, including $\unicode[STIX]{x1D70F}_{yy}$ , which in turn generates the turbulent shear stress $\unicode[STIX]{x1D70F}_{xy}$ through interactions with the mean strain field. Another possible cause for the numerical instability is the ill-conditioning of the RANS momentum equations themselves. Here model conditioning refers to the sensitivity of the solved quantities (e.g. mean velocity and pressure fields) to the modelled terms (Reynolds stresses). In computational fluid dynamics codes that solve RANS equations and turbulence transport equations in a segregated manner, ill-conditioned systems lead to numerical instability due to the sensitivity of solved velocities to residuals in the Reynolds stress from iteration to iteration. While the chain coupling mechanism described above is well known in the turbulence modelling literature, the conditioning of RANS equations has rarely been mentioned. This is probably due to the fact that the two causes are closely intertwined in traditional models based on Reynolds stress transport equations and the former dominated.

1.1 Unique challenges in data-driven Reynolds stress closures

Recently, data-driven turbulence modelling has emerged as a promising field of research with a number of approaches having been proposed in the past few years (e.g. Ling, Kurzawski & Templeton Reference Ling, Kurzawski and Templeton2016; Singh & Duraisamy Reference Singh and Duraisamy2016; Weatheritt & Sandberg Reference Weatheritt and Sandberg2016; Wang, Wu & Xiao Reference Wang, Wu and Xiao2017). Of particular relevance to the present discussions are the data-driven Reynolds stress closures (Ling et al. Reference Ling, Kurzawski and Templeton2016; Wang et al. Reference Wang, Wu and Xiao2017; Geneva & Zabaras Reference Geneva and Zabaras2019) in which the Reynolds stresses are obtained directly from machine learning models trained on high-fidelity simulation databases without solving partial differential equations (PDEs) or using explicit algebraic models. In such models, the chain coupling mechanism related to the Reynolds stress transport equations as described above is not a relevant cause of numerical instability, and thus the possible ill-conditioning of RANS equations is placed under the spotlight. Moreover, these data-driven Reynolds stress models do not have explicit expressions for the Reynolds stress (Ling et al. Reference Ling, Kurzawski and Templeton2016; Wang et al. Reference Wang, Wu and Xiao2017), which makes it difficult to treat the Reynolds stresses implicitly in the RANS equations to improve model conditioning. For example, Wang et al. (Reference Wang, Wu and Xiao2017) used random forest regression to predict Reynolds stresses based on direct numerical simulation (DNS) databases and reported that the solved mean velocity field does not improve over the original RANS simulations, even though the predicted Reynolds stresses showed significant improvements. Such a paradox clearly demonstrates the gap between a priori and a posteriori performances in turbulence models based on data-driven Reynolds stress closures. That is, an ill-conditioned model can amplify small a priori errors in the modelled terms to large a posteriori errors in the solved quantities, which defeats all efforts in the improvement of the closure models. Similar gaps have been observed in data-driven subgrid-scale (SGS) stress models for large-eddy simulations (LES). For example, Gamahara & Hattori (Reference Gamahara and Hattori2017) reported that their neural network model predicted better SGS stresses than the Smagorinsky model for turbulent channel flows, but the mean velocity predictions were much less satisfactory. Given the large number of efforts in developing data-driven Reynolds stresses model (e.g. Ling et al. Reference Ling, Kurzawski and Templeton2016; Wang et al. Reference Wang, Wu and Xiao2017; Zhu et al. Reference Zhu, Zhang, Kou and Liu2019) and SGS models (e.g. King, Hamlington & Dahm Reference King, Hamlington and Dahm2016; Maulik & San Reference Maulik and San2017; Vollant, Balarac & Corre Reference Vollant, Balarac and Corre2017; Maulik et al. Reference Maulik, San, Rasheed and Vedula2019) for turbulent flows, the conditioning in such models must be closely examined along with their a priori predictive performances. Note, however, that some data-driven approaches, e.g. those based on correcting turbulence transport equations (Parish & Duraisamy Reference Parish and Duraisamy2016; Singh, Medida & Duraisamy Reference Singh, Medida and Duraisamy2017) or discovering analytical turbulent constitutive relations by using symbolic regression (Weatheritt & Sandberg Reference Weatheritt and Sandberg2016), are less affected by the conditioning discussed here.

1.2 DNS data as the ideal scenario for data-driven Reynolds stress closures

Using Reynolds stress obtained from DNS data to solve the RANS equations for velocity can be considered the ideal scenario for data-driven Reynolds stress closure, at least for those without analytical expressions. Solving for mean velocities with a given Reynolds stress field is referred to as ‘propagation’ in this work, i.e. propagation of Reynolds stresses to mean velocities by solving the RANS equations. Such a methodology represents an absolute upper limit of performances for data-driven Reynolds stress models as pursued by Ling et al. (Reference Ling, Kurzawski and Templeton2016), Wang et al. (Reference Wang, Wu and Xiao2017) and Zhu et al. (Reference Zhu, Zhang, Kou and Liu2019). It allows us to concentrate on the conditioning of this class of machine-learning-based Reynolds stress models without considering the a priori performance of any specific model.

Somewhat surprisingly, even solving RANS equations with Reynolds stresses from DNS data can lead to large errors in the velocities, which has been demonstrated in the recent work of Thompson et al. (Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016). They propagated Reynolds stresses obtained from DNS to mean velocities for turbulent channel flows at a wide range of Reynolds numbers. These DNS were performed with extreme caution by reputable groups (Del Alamo & Jiménez Reference Del Alamo and Jiménez2003; Bernardini, Pirozzoli & Orlandi Reference Bernardini, Pirozzoli and Orlandi2014; Lee & Moser Reference Lee and Moser2015), and it was verified that the errors in the reported Reynolds stresses were indeed very small, typically less than 0.5 % as shown in table 1. Thompson et al. (Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016) showed that the solved mean velocity has an unsatisfactory agreement with the DNS data at high frictional Reynolds numbers (e.g.  $Re_{\unicode[STIX]{x1D70F}}=5200$ ). Poroseva et al. (Reference Poroseva, Colmenares, Juan and Murman2016) also confirmed such observations. To motivate our work, we reproduced the two studies of solving for mean velocities by using the DNS Reynolds stresses obtained from Lee & Moser (Reference Lee and Moser2015), and the results are summarized in table 1. Although the solved mean velocities shown here are accurate up to $Re_{\unicode[STIX]{x1D70F}}=1000$ , researchers have reported that large errors in the propagated velocities can be found at Reynolds number as low as $Re_{\unicode[STIX]{x1D70F}}=395$ depending on the DNS data used (S. Poroseva, 2017, personal communication).

Table 1.

Summary of the results of the channel flow test, showing percentage errors in the turbulent shear stresses and the propagated mean velocities. The large errors in the high-Reynolds-number case $Re_{\unicode[STIX]{x1D70F}}=5200$ are highlighted in bold. The true Reynolds shear stresses are obtained by analytically integrating the RANS equation with the DNS mean velocities (Thompson et al. Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016).

The percentage errors shown in table 1 are computed by comparing the DNS Reynolds stresses and the mean velocities propagated therefrom against their respective truths, which are obtained according to Thompson et al. (Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016). Specifically, the true velocities are taken as the DNS mean velocities, and the true Reynolds shear stresses are computed based on the analytical solution

(1.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\left(1-\frac{y}{h}\right)\unicode[STIX]{x1D70F}_{w}-\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\frac{\text{d}U}{\text{d}y}\end{eqnarray}$$

from the DNS mean velocities and wall shear stresses $\unicode[STIX]{x1D70F}_{w}$ , where $y$ denotes the wall distance ( $y=0$ indicates the wall), $h$ denotes the half channel width, and $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D708}$ represent density and kinematic viscosity, respectively.

It can be seen in table 1 that the errors in the Reynolds stresses are less than 0.5 %. The errors in the Reynolds stress can be attributed to different sources, e.g. statistical sampling errors (due to inadequate averaging of the DNS results), iterative errors (due to lack of convergence when solving the linear equations in DNS) and interpolation errors (when interpolating the DNS data to the RANS mesh). All these errors are inevitable in any DNS data. It is noted that errors in the Reynolds stresses vary non-monotonically with respect to the Reynolds number. Such a coincidental trend should not be overly or literally interpreted, since the data shown here are a result of complex interactions that lead to superposition and cancellation among various sources as discussed above. It suffices to point out that all the DNS were obtained with extreme caution, as is evident from the very small errors in the Reynolds stresses. Nevertheless, the errors in the mean velocity solved from the Reynolds stresses clearly increase monotonically with the Reynolds number, because they are dominated by the conditioning of RANS equations, as we will show later. Specifically, while the errors in the solved mean velocity are also less than 0.5 % at low Reynolds number ( $Re_{\unicode[STIX]{x1D70F}}=180$ ), such errors can be as high as 35 % at high Reynolds number ( $Re_{\unicode[STIX]{x1D70F}}=5200$ ). This observation suggests that the RANS equations can be ill-conditioned by directly substituting the modelled Reynolds stress into the equations, which is a common practice in current data-driven RANS modelling. Such results raise several critical questions on Reynolds-stress-based data-driven turbulence models, specifically the following:

  1. (i) How to explain the deteriorated conditioning (i.e. increased amplification of errors in Reynolds stresses to velocities) with increasing Reynolds numbers observed in these studies?

  2. (ii) Is there a quantitative metric that can characterize the conditioning of RANS equations with a given turbulence model?

  3. (iii) What are the implications of the above observations to data-driven turbulence models that treat Reynolds stresses as source terms in the RANS equations?

1.3 Summary and contribution of present work

Our present work aims to answer the above questions by proposing a quantitative metric and elucidating the relevance of the above-mentioned studies to data-driven turbulence modelling. Throughout this study, no turbulence models are used. Rather, Reynolds stresses obtained from DNS databases are used to represent the ideal performance for any data-driven Reynolds stress closures. A local condition number function based on the work of Chandrasekaran & Ipsen (Reference Chandrasekaran and Ipsen1995) is derived as a metric to assess the conditioning property for turbulence models. Numerical tests on turbulent channel flows at various Reynolds numbers and two more complex flows suggest that the proposed metric can adequately explain observations in previous studies, i.e. deteriorated model conditioning with increasing Reynolds number, and better conditioning of the implicit treatment of the Reynolds stress compared to the explicit treatment. As an application of the proposed approach, Wu, Xiao & Paterson (Reference Wu, Xiao and Paterson2018) improved the conditioning of the data-driven Reynolds stress model of Wang et al. (Reference Wang, Wu and Xiao2017) by training machine learning models for the linear and nonlinear parts of the Reynolds stress separately. The obtained Reynolds stress model achieved good conditioning and satisfactory predictive accuracy simultaneously.

The rest of this paper is organized as follows. Section 2 introduces the global condition number and shows that it fails to explain the deteriorated conditioning with increasing Reynolds numbers. A local condition number function is derived to achieve such an objective. In § 3, the local condition number is used to evaluate the conditioning of RANS equations with data-driven Reynolds stress closures. Section 4 discusses the conditioning of RANS equations with traditional Reynolds stress transport models and the monolithic coupling for both traditional and data-driven models. Finally, conclusions are presented in § 5.

2 Methodology

Consider the steady-state RANS equations for incompressible, constant-density fluids:

(2.1) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}=0, & \displaystyle\end{eqnarray}$$
(2.2) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the mean flow velocity, $\unicode[STIX]{x1D708}$ is molecular viscosity, $p$ is the pressure normalized by the constant density of the fluid, and $\unicode[STIX]{x1D749}$ is the Reynolds stress tensor, which needs to be modelled. For simplicity, we first consider a Reynolds-stress-based model where $\unicode[STIX]{x1D749}$ is obtained by solving a transport equation in a segregated manner with the RANS equations or by a data-driven function (see e.g. Ling et al. Reference Ling, Kurzawski and Templeton2016). The objective of this work is to investigate the sensitivity of the obtained mean velocity with respect to small perturbations on the Reynolds stress.

For simplicity of notation, we introduce nonlinear operator ${\mathcal{N}}$ to include the convection and diffusion terms with

(2.3) $$\begin{eqnarray}{\mathcal{N}}(\boldsymbol{u})=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}.\end{eqnarray}$$

The RANS momentum equation in (2.1) can be written as

(2.4) $$\begin{eqnarray}{\mathcal{N}}(\boldsymbol{u})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}-\unicode[STIX]{x1D735}p.\end{eqnarray}$$

In numerical solvers, the convection term is first linearized around the current velocity $\overline{\boldsymbol{U}}_{0}$ to obtain the linearized RANS equations

(2.5) $$\begin{eqnarray}{\mathcal{L}}(\boldsymbol{u})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}-\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where ${\mathcal{L}}$ is the linearized operator of ${\mathcal{N}}$ , i.e.

(2.6) $$\begin{eqnarray}{\mathcal{L}}(\boldsymbol{u})=\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}.\end{eqnarray}$$

The linearized equation (2.5) is then discretized on a computational fluid dynamics (CFD) mesh to obtain a linear system of the following form:

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D63C}\boldsymbol{U}=\boldsymbol{b}.\end{eqnarray}$$

Here we have denoted $\boldsymbol{b}=[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}-\unicode[STIX]{x1D735}p]$ as the imbalance between the two forces, pressure gradient and Reynolds stress divergence; and $\boldsymbol{U}=[\boldsymbol{u}]$ is the discretized velocity field to be solved for. Both $\boldsymbol{b}$ and $\boldsymbol{U}$ are $n\times 1$ vectors, where $n$ is the number of cells or grid points in the mesh. The matrix $\unicode[STIX]{x1D63C}$ with dimension $n\times n$ comes from the implicit discretization of the linearized convection term and the molecular diffusion term. In this work, we focus on the conditioning of linearized RANS equations. This is because most CFD codes deal with the linearized RANS equations in each iteration when solving the nonlinear RANS equations. Therefore, the conditioning metrics studied here are valid within each iteration, even though the flow of concern may deviate from the linearized RANS equations.

2.1 Matrix norm as a measure of model conditioning

We first show the derivation of the traditional matrix-norm-based condition number and explain why it fails to distinguish the sensitivities of solving for mean velocity at different Reynolds numbers as shown in table 1. Following the definition of matrix norm, the norm of the error in the velocity is bounded as follows:

(2.8) $$\begin{eqnarray}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{U}\Vert }{\Vert \boldsymbol{U}\Vert }\leqslant {\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{b}\Vert }{\Vert \boldsymbol{b}\Vert },\end{eqnarray}$$

where

(2.9) $$\begin{eqnarray}{\mathcal{K}}_{\unicode[STIX]{x1D63C}}\equiv \Vert \unicode[STIX]{x1D63C}\Vert \,\Vert \unicode[STIX]{x1D63C}^{-1}\Vert\end{eqnarray}$$

denotes the condition number of matrix $\unicode[STIX]{x1D63C}$ (see e.g. Strang Reference Strang2016).

As explained in the notation, the norms $\Vert \overline{\boldsymbol{U}}\Vert$ and $\Vert \boldsymbol{b}\Vert$ are taken of the discretized vectors $[\overline{\boldsymbol{U}}]$ and $[\boldsymbol{b}]$ , respectively, with the brackets inside the norm omitted for clarity. Such a brief notation does not cause confusion, because norms in this work are always taken for the discretized vectors or matrices with dimensions of $n\times 1$ or $n\times n$ , respectively, and never for the velocity or force vectors at any particular location.

Considering that the objective is to assess the effects of Reynolds stress perturbation $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}$ on the velocities, the inequality in (2.8) above is formulated as

(2.10) $$\begin{eqnarray}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{U}\Vert }{\Vert \boldsymbol{U}\Vert }\leqslant {\mathcal{K}}_{\unicode[STIX]{x1D70F}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert }{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert },\end{eqnarray}$$

where

(2.11) $$\begin{eqnarray}{\mathcal{K}}_{\unicode[STIX]{x1D70F}}={\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert }{\Vert \boldsymbol{b}\Vert }.\end{eqnarray}$$

Detailed derivations are omitted here for brevity and are presented in § A.1. It can be seen that the model condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ consists of the condition number of the matrix $\unicode[STIX]{x1D63C}$ and the ratio

(2.12) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6FC}}=\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert /\Vert \boldsymbol{b}\Vert .\end{eqnarray}$$

For plane channel flows, the convective term disappears, and thus $\boldsymbol{b}$ is the force due to the divergence of the viscous stress, i.e. $\boldsymbol{b}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D708}(\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{vis}$ . Consequently, the ratio $\overline{\unicode[STIX]{x1D6FC}}$ indicates the overall relative importance of the forces due to Reynolds stress and viscous stress.

The proposed condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ based on matrix condition number ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ is a natural first attempt in explaining the increasing sensitivity of the velocities to the Reynolds stress with increasing Reynolds numbers as shown in table 1. However, surprisingly, it turns out that the condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ is more or less the same across all Reynolds numbers from $Re_{\unicode[STIX]{x1D70F}}=180$ to 5200, which is shown in figure 1. This observation suggests that the matrix-based condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ cannot explain the ill-conditioning of the $Re_{\unicode[STIX]{x1D70F}}=5200$ case and the better conditioning of the lower-Reynolds-number cases as observed in table 1.

Figure 1.

Conditioning measure of Reynolds-stress-based turbulence models based on ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ and the ratio $\overline{\unicode[STIX]{x1D6FC}}$ as defined in (2.12). The Reynolds stress is computed from DNS data to study the ideal scenario of the RANS modelling.

The following two factors explain why the matrix-based condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ is almost the same at different Reynolds numbers:

  1. (i) The matrix condition number ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ is constant for all Reynolds numbers, because the matrix $\unicode[STIX]{x1D63C}$ itself is independent of the Reynolds number.

  2. (ii) The ratios $\overline{\unicode[STIX]{x1D6FC}}=\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert /\Vert \boldsymbol{b}\Vert$ are very similar at vastly different Reynolds numbers, because both norms (which involve integration or square sums) are dominated by the viscous wall regions of each flow. It is well known that the Reynolds number only determines the thickness of this region in outer coordinates, and the Reynolds-number effect is weak here. Each factor above will be detailed as follows.

First, it can be established through simple algebra that for plane channel flows the matrix $\unicode[STIX]{x1D63C}$ is independent of the Reynolds number but depends on the discretization scheme and the mesh used. Since the convection term disappears for plane channel flows, the matrix $\unicode[STIX]{x1D63C}$ results solely from the discretization of the diffusion operator $\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}(\cdot )$ . When discretized with central difference on a uniform mesh of $n$ cells, matrix  $\unicode[STIX]{x1D63C}$ can be written as follows (Strang Reference Strang2016):

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D708}\left[\begin{array}{@{}ccccc@{}}2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & \ddots & \\ & & \ddots & \ddots & -1\\ & & & -1 & 2\end{array}\right].\end{eqnarray}$$

The condition number for matrix $\unicode[STIX]{x1D63C}$ is ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}=4n^{2}/\unicode[STIX]{x03C0}^{2}$ . Therefore, ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ does not explicitly depend on the viscosity or the Reynolds number. This analysis is confirmed by the results shown in figure 1, which shows that ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ is strictly constant for all five flows at different Reynolds numbers. Moreover, ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ depends on the mesh size  $n$ , which is a critical shortcoming of the matrix-norm-based condition number as a measure of the conditioning property of a turbulence model. To exclude the influences of the mesh, we used the same mesh with 1040 cells in all the flows at different Reynolds numbers presented in figure 1.

Figure 2.

The balance among forces due to turbulent shear stress ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}$ ), viscous shear stress ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{vis}$ ) and pressure gradient ( $\unicode[STIX]{x1D735}p$ ) for two plane channel flows at frictional Reynolds numbers (a)  $Re_{\unicode[STIX]{x1D70F}}=180$ and (b)  $Re_{\unicode[STIX]{x1D70F}}=5200$ . The right vertical axis denotes the inner coordinates ( $y^{+}$ ).

Figure 3.

Force balance of the plane channel flow in (a) the outer layer and (b) the viscous wall region.

Second, the change of Reynolds number has little influence on the ratio $\overline{\unicode[STIX]{x1D6FC}}$ , as is shown in figure 1. We examine the profile of turbulent shear ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}$ ) and viscous shear ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{vis}$ ) in the channel in figure 2 for the two cases, $Re_{\unicode[STIX]{x1D70F}}=180$ and 5200. In most of the channel outside the viscous wall region, both forces (and thus the ratio) are fairly uniform. Nevertheless, in the viscous wall region, the two forces are of the same order of magnitude but with opposite signs. In contrast, outside the viscous region, the pressure gradient is the main driving force while the Reynolds shear stress is the resistance, with the viscous shear having negligible effects. The forces in the two distinct regions are illustrated schematically in figure 3. However, when calculating the ratio $\overline{\unicode[STIX]{x1D6FC}}=\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert /\Vert \boldsymbol{b}\Vert$ of the two norms, values within the viscous wall region clearly dominate the calculation of both norms, which involve integration of the functions squared. It can be seen that in both figure 2(a) ( $Re_{\unicode[STIX]{x1D70F}}=180$ ) and figure 2(b) ( $Re_{\unicode[STIX]{x1D70F}}=5200$ ) the areas enclosed by the blue solid curve and red dashed curve (with the vertical zero line) are similar. Squaring the function places even more weights on the regions of larger function values, i.e. the viscous wall region. This observation suggests that the ratio $\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert /\Vert \boldsymbol{b}\Vert$ is of $O(1)$ for both cases, as confirmed by examining figure 1. Consequently, the computed norm mostly reflects the values of the forces in the viscous region and not the outer layer. It is well known that the Reynolds-number effects are not pronounced within the viscous wall region. Increasing the Reynolds number merely extends the outer layer in terms of inner coordinates ( $y^{+}=y/y^{\ast }$ , where $y^{\ast }=\unicode[STIX]{x1D708}/\sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}}$ is the viscous unit and $\unicode[STIX]{x1D70F}_{w}$ is the wall shear stress). This explains why the ratio $\overline{\unicode[STIX]{x1D6FC}}$ does not vary significantly (much less than proportionally) with the Reynolds number as can be seen in figure 1. If the viscous region is neglected, the factor $\bar{\unicode[STIX]{x1D6FC}}$ will be higher for the larger Reynolds number, and thus it makes the global condition number a better indicator of model conditioning. However, the matrix condition number ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ depends on the mesh size and thus the global condition number is still not an ideal choice of evaluating model conditioning of RANS equations.

In summary, the matrix-based condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ as derived in (2.11) is not able to explain the increasing sensitivity of the velocities with respect to Reynolds stresses with increasing Reynolds number. In addition, the matrix condition number ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ has another critical drawback of being mesh-dependent. The mesh dependence is highly undesirable, as the condition number is to measure the conditioning property of turbulence models at the PDE level, not any particular numerical discretization thereof. These drawbacks clearly call for a better metric for measuring the conditioning property of Reynolds-stress-based turbulence models.

2.2 Proposed metric as a measure of model conditioning

In order to address the deficiency of the global condition number as presented in § 2.1, we derive a metric based on a local condition number function to measure the sensitivity of the solved mean velocity $\boldsymbol{u}$ at a given location $\boldsymbol{x}$ with respect to perturbation $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}$ on the Reynolds stress field  $\unicode[STIX]{x1D749}$ . Such a local condition number is formally defined as the following bound:

(2.14) $$\begin{eqnarray}\frac{|\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x})|}{U_{\infty }}\leqslant {\mathcal{K}}(\boldsymbol{x})\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}},\end{eqnarray}$$

where ${\mathcal{K}}(\boldsymbol{x})$ is the local condition number function defined as

(2.15) $$\begin{eqnarray}{\mathcal{K}}(\boldsymbol{x})=\frac{\Vert G(\boldsymbol{x},\unicode[STIX]{x1D743})\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}.\end{eqnarray}$$

Function $G$ is the Green’s function corresponding to the linear operator ${\mathcal{L}}$ (see e.g. Lanczos Reference Lanczos1996), such that the solution to the linearized RANS equation (2.5) can be written formally as

(2.16) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x})={\mathcal{L}}^{-1}[\boldsymbol{b}(\boldsymbol{x})]\equiv \int _{\unicode[STIX]{x1D6FA}}G(\boldsymbol{x};\unicode[STIX]{x1D743})\boldsymbol{b}(\unicode[STIX]{x1D743})\,\text{d}\unicode[STIX]{x1D743},\end{eqnarray}$$

where ${\mathcal{L}}^{-1}[\cdot ]$ is the inverse operator of ${\mathcal{L}}$ . Green’s function $G(\boldsymbol{x};\unicode[STIX]{x1D743})$ indicates the contribution of the source $\boldsymbol{b}(\unicode[STIX]{x1D743})$ at location $\unicode[STIX]{x1D743}$ to the solution $\boldsymbol{u}$ at location  $\boldsymbol{x}$ . The norm $\Vert \,f(\unicode[STIX]{x1D743})\Vert _{\unicode[STIX]{x1D6FA}}$ of function $f(\unicode[STIX]{x1D743})$ is an integration on domain $\unicode[STIX]{x1D6FA}$ defined as (Debnath & Mikusiński Reference Debnath and Mikusiński2005)

(2.17) $$\begin{eqnarray}\Vert \,f(\unicode[STIX]{x1D743})\Vert _{\unicode[STIX]{x1D6FA}}=\sqrt{\int _{\unicode[STIX]{x1D6FA}}|\,f(\unicode[STIX]{x1D743})|^{2}\,\text{d}\unicode[STIX]{x1D743}}.\end{eqnarray}$$

The detailed derivations to obtain (2.15) are presented in § A.2.

For functions discretized on a CFD mesh of $n$ cells (e.g. those in RANS simulations), the function norm $\Vert \cdot \Vert _{\unicode[STIX]{x1D6FA}}$ in (2.15) can be approximated by the norm of the discretized $n$ -vector through numerical quadrature. That is,

(2.18a,b ) $$\begin{eqnarray}\Vert G(\boldsymbol{x},\unicode[STIX]{x1D743})\Vert _{\unicode[STIX]{x1D6FA}}\approx \Vert \boldsymbol{r}_{j}\Vert _{n}\quad \text{and}\quad \Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}\approx \Vert [\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]\Vert _{n},\end{eqnarray}$$

where $\boldsymbol{r}_{j}$ is the $j$ th row of the matrix $\unicode[STIX]{x1D63C}^{-1}$ . Recall that $[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]$ indicates discretization of field $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}$ on the CFD mesh, but the bracket can be omitted inside a vector norm $\Vert \cdot \Vert _{n}$ without ambiguity. In this work, the norm $\Vert \cdot \Vert _{n}$ is defined by the $L^{2}$ -norm as

(2.19) $$\begin{eqnarray}\Vert \boldsymbol{v}\Vert _{n}=\left(\mathop{\sum }_{i}v_{i}^{2}\right)^{1/2}\end{eqnarray}$$

and the discretized condition number $n$ -vector corresponding to ${\mathcal{K}}(\boldsymbol{x})$ in (2.15) is thus

(2.20) $$\begin{eqnarray}{\mathcal{K}}_{j}=\frac{\Vert \boldsymbol{r}_{j}\Vert _{n}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{n}}{U_{\infty }}\quad \text{with}~j=1,2,\ldots ,n,\end{eqnarray}$$

which implies that the location $\boldsymbol{x}$ is the coordinate of the $j$ th cell in the CFD mesh.

The proposed local condition number function ${\mathcal{K}}(\boldsymbol{x})$ has two important merits compared to the global matrix-based condition number  ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ :

  1. (i) First, ${\mathcal{K}}(\boldsymbol{x})$ provides a tighter upper bound than the matrix-norm condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ . The main reason is that the upper bound of ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ can only be achieved when the following conditions are satisfied simultaneously: (1) the discretized mean velocity field vector $\boldsymbol{U}$ is aligned with the principal axis of the coefficient matrix $\unicode[STIX]{x1D63C}$ , and (2) the perturbation vector $\unicode[STIX]{x1D6FF}\boldsymbol{b}$ is aligned with the principal axis of $\unicode[STIX]{x1D63C}^{-1}$ . In contrast, the derivation of ${\mathcal{K}}(\boldsymbol{x})$ does not assume any conditions on the discretized mean velocity field  $\boldsymbol{U}$ . Consequently, the bound provided by ${\mathcal{K}}(\boldsymbol{x})$ is a more precise assessment of the sensitivity $\unicode[STIX]{x1D6FF}\boldsymbol{u}$ with respect to Reynolds stress perturbations.

  2. (ii) The discretization ${\mathcal{K}}_{j}$ of function ${\mathcal{K}}(\boldsymbol{x})$ is mesh-independent, which is an important property considering that this metric aims to measure the conditioning property of data-driven Reynolds stress models. Detailed analytical derivations to obtain (2.20) and numerical results (see figure 17 in the appendix) used to demonstrate the mesh independence of the local condition number $K_{j}$ are presented in appendix B.

While the local condition metrics proposed here are generally applicable to any linearized PDE and its discretization, it is important to emphasize that these conditioning metrics faithfully embody the physics described by the underlying PDEs. Taking the RANS equations as an example, the mean flow field contributes to the linearized differential operator and thus is reflected in the analysis of local conditioning. Therefore, these local condition metrics are flow-specific properties and thus reflect the mean flow physics. More detailed discussion can be found in §§ 3.2.1 and 3.2.2 on the complex, two-dimensional flows.

Finally, we comment briefly on the numerical implementation and computational complexity of conditioning metrics proposed below. The local condition number ${\mathcal{K}}_{j}$ does not require a full inversion of the matrix $\unicode[STIX]{x1D63C}$ but only needs the $j$ th row of $\unicode[STIX]{x1D63C}^{-1}$ to be computed. This can be achieved by solving the equation $\unicode[STIX]{x1D63C}^{\text{T}}\boldsymbol{r}_{j}=\unicode[STIX]{x1D644}_{j}$ based on the identity $\unicode[STIX]{x1D63C}^{-1}\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D644}$ , where $\unicode[STIX]{x1D644}_{j}$ denotes the $j$ th column of the identity matrix $\unicode[STIX]{x1D644}$ . Solving this equation is a standard, inexpensive routine available in CFD codes, e.g. with a computational complexity of $O(n\log n)$ if a multigrid linear solver is used, where $n$ indicates the number of elements in $\boldsymbol{r}$ . Therefore, it can be estimated that even obtaining the full condition number field only has a computational complexity of $O(n^{2}\log n)$ with a standard multigrid linear solver, which is much lower than the complexity of $O(n^{3})$ for typical algorithms of matrix inversions.

As the local condition number ${\mathcal{K}}(\boldsymbol{x})$ is a spatial function, and its discretization is an $n$ -vector, it is desirable to obtain a scalar quantity to provide an integral, more straightforward measure of model conditioning property similar to the global condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ . To this end, we define a volume-averaged local condition number $\overline{{\mathcal{K}}}_{\boldsymbol{x}}$ defined as

(2.21) $$\begin{eqnarray}\overline{{\mathcal{K}}}_{\boldsymbol{x}}=\frac{\displaystyle \mathop{\sum }_{j=1}^{n}[{\mathcal{K}}_{j}][\unicode[STIX]{x0394}V_{j}]}{V},\end{eqnarray}$$

where $\unicode[STIX]{x0394}V_{j}$ denotes the volume of the $j$ th cell in the CFD mesh, and $V$ is the total volume of the computational domain. This volume-averaged local condition number $\overline{{\mathcal{K}}}_{\boldsymbol{x}}$ has a similar interpretation to ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ but preserves the merits of ${\mathcal{K}}_{j}$ , i.e. tighter bounds and mesh independence.

In the derivations above, the Reynolds stress term is substituted directly into the RANS equation and is treated explicitly. When the Reynolds stress term is treated implicitly as in many practical implementations of Reynolds stress models (e.g. Basara & Jakirlic Reference Basara and Jakirlic2003; Maduta & Jakirlic Reference Maduta and Jakirlic2017), the corresponding local condition number of the model can be obtained similarly, except that the Green’s function is modified to account for the implicit modelling of the linear part of the Reynolds stress with the eddy-viscosity model. Specifically, the general form of implicit treatment of the Reynolds stress can be written as

(2.22) $$\begin{eqnarray}\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D708}_{t}\unicode[STIX]{x1D64E}+\unicode[STIX]{x1D749}^{\bot },\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{t}$ represents the eddy viscosity, $\unicode[STIX]{x1D64E}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}})$ denotes the strain-rate tensor and $\unicode[STIX]{x1D749}^{\bot }$ denotes the nonlinear part. In this work, the eddy viscosity $\unicode[STIX]{x1D708}_{t}$ is obtained from projecting DNS Reynolds stress onto DNS strain-rate tensor and not from RANS simulations. With such an optimal eddy viscosity  $\unicode[STIX]{x1D708}_{t}^{m}$ , equation (2.22) treats the linear part of the Reynolds stress tensor implicitly to enhance the conditioning of RANS equations. Consequently, the Green’s function $\widetilde{G}$ corresponding to the linear operator

(2.23) $$\begin{eqnarray}\widetilde{{\mathcal{L}}}(\boldsymbol{u})={\mathcal{L}}(\boldsymbol{u})-\unicode[STIX]{x1D708}_{t}^{m}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-(\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t}^{m})\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}\end{eqnarray}$$

should be used in (2.15) and (2.20), with $\widetilde{G}$ related to $\widetilde{{\mathcal{L}}}$ in a similar way as $G$ is to ${\mathcal{L}}$ in (2.16). The optimal eddy viscosity $\unicode[STIX]{x1D708}_{t}^{m}$ is computed by minimizing the discrepancy between the linear eddy-viscosity model and the DNS Reynolds stress data, i.e. 

(2.24) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{t}^{m}(\boldsymbol{x})=\mathop{\{}_{\unicode[STIX]{x1D708}_{t}}\Vert \unicode[STIX]{x1D749}^{DNS}-2\unicode[STIX]{x1D708}_{t}(\boldsymbol{x})\unicode[STIX]{x1D64E}^{DNS}\Vert ,\end{eqnarray}$$

where $\unicode[STIX]{x1D749}^{DNS}$ and $\unicode[STIX]{x1D64E}^{DNS}$ denote the Reynolds stress and the strain-rate tensor from the DNS database, respectively. Here we emphasize that the optimal eddy viscosity is location-dependent. The detailed derivations are presented in § A.3. It is noted that the eddy viscosity $\unicode[STIX]{x1D708}_{t}$ in (2.22) only quantifies the amount of Reynolds stress being treated implicitly and is not necessarily the optimal eddy viscosity  $\unicode[STIX]{x1D708}_{t}^{m}$ . By specifying different  $\unicode[STIX]{x1D708}_{t}$ , the amount of implicit treatment of Reynolds stress and the amount of the nonlinear part of the Reynolds stress vary accordingly, leading to different conditioning of RANS equations. Therefore, equation (2.22) provides a general form of the implicit treatment of Reynolds stresses. The optimal eddy viscosity is capped to be positive for numerical stability of solving the RANS equations. Therefore, the implicit treatment always leads to better conditioning of the system, but the difference between the implicit and explicit treatments becomes smaller in the regions where the linear part of the Reynolds stress is less dominant. In the extreme (albeit unlikely) situation where the Reynolds stress tensor is orthogonal to the strain-rate tensor (i.e. optimal eddy viscosity is zero across the flow domain), the conditioning with implicit and explicit treatments would be equivalent.

In summary, we proposed a local condition number to assess the sensitivity of local mean velocities with regard to data-driven Reynolds stress models. It has the following three forms: the spatial function ${\mathcal{K}}(\boldsymbol{x})$ (i.e. condition number function), an $n$ -vector ${\mathcal{K}}_{j}$ obtained by discretizing ${\mathcal{K}}(\boldsymbol{x})$ on the CFD mesh, and a scalar $\overline{{\mathcal{K}}}_{\boldsymbol{x}}$ obtained by integration of ${\mathcal{K}}(\boldsymbol{x})$ . This metric is applicable to different types of data-driven Reynolds stress models. The main contributions of this work are (1) highlighting the importance of the conditioning of PDE-governed systems with data-driven closure models, and (2) providing a quantitative metric for assessing such conditioning. Data-driven modelling is becoming an emerging area in the fluid mechanics community. However, all existing works of data-driven modelling focus on the accuracy of the model itself. Our work demonstrates that an accurate data-driven model does not necessarily guarantee satisfactory predictions of quantities of interest. Therefore, ensuring good conditioning of the problem formulation (i.e. PDEs with data-driven closure) is as important as improving the accuracy of data-driven model. We envision a standard practice in the future where all developed data-driven closures are presented along with corresponding conditioning analysis, in the same way in which experimental data reported nowadays are accompanied by their associated uncertainties.

3 Results

We first use the proposed local condition number to explain the model conditioning of RANS equations when using an explicit treatment with fixed Reynolds stress. Specifically, turbulent channel flows at different Reynolds numbers are studied and the results are shown in table 1. In addition, the implicit treatment of the Reynolds stress is studied to compare with the model conditioning of RANS equations by using an explicit treatment with fixed Reynolds stress. By studying these two types of RANS modelling, we show that the proposed local condition number can be used to assess the sensitivities of RANS simulations for data-driven modelling. We also extend the study of the model conditioning and the proposed local condition number to two more complex flows, including the flow over periodic hills and the flow in a square duct. The results show that the RANS equations can be ill-conditioned for other types of flows, which can be assessed by the proposed condition number metric.

The RANS simulations are performed in a finite-volume CFD platform OpenFOAM, using a modified flow solver that allows the explicit and implicit treatments of the Reynolds stress computed from DNS data. For numerical discretizations, the second-order central difference scheme is chosen for all terms except for the convection term, which is discretized with a second-order upwind scheme. The second-order upwind scheme was used to avoid the possible numerical instability when using the central difference scheme for the convection term. For the turbulent plane channel cases, the convection term disappears at steady state and thus the results presented in § 3.1 are not influenced by the discretization of the convection term thereof. In § 3.1, the mean velocity is obtained by directly solving (2.5) since the mean velocity and the pressure are decoupled for the RANS simulation of a fully developed plane channel flow. The convergence criterion of solving the momentum equations is set as  $10^{-8}$ in absolute error.

3.1 Turbulent channel flow

The fully developed turbulent plane channel flows are investigated by using the local condition number  ${\mathcal{K}}_{j}$ . In this work, we consider an ideal scenario in which the Reynolds stress $\unicode[STIX]{x1D749}$ is directly computed from a DNS database at various Reynolds numbers $Re_{\unicode[STIX]{x1D70F}}=180$ , 550, 1000, 2000 and 5200. The numbers of mesh cells $N$ are 36, 110, 200, 400 and 1040, respectively, for the channel flows at those Reynolds numbers. Non-uniform meshes are used and the expansion ratio is adjusted to ensure that the $y^{+}$ of the first cell centre is kept below 1. Therefore, the mesh size of the first cell in the wall-normal direction is below $h/N$ , where $h$ denotes the half channel height. The DNS data were obtained from the University of Texas, Austin, online database (Lee & Moser Reference Lee and Moser2015). The mean velocity field is then solved by substituting the computed Reynolds stress as the closure term of RANS equations.

In the practice of RANS modelling, iterations are involved in solving RANS equations and the modelling of the Reynolds stress is updated by the mean velocity field during the iterations. Therefore, it is possible that the mean velocity field and the Reynolds stress can adjust to each other during the iterations. We employed the ratio $\unicode[STIX]{x1D6FF}U_{rms}/U_{rms}^{DNS}$ to assess the error of the solved mean flow field at each iteration step. Specifically, the volume-averaged root-mean-squared (r.m.s.) error of the solved mean velocity is defined as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}U_{rms}=\sqrt{\frac{\displaystyle \mathop{\sum }_{j=1}^{n}([U]_{j}-[U^{DNS}]_{j})^{2}[\unicode[STIX]{x0394}V_{j}]}{V}}.\end{eqnarray}$$

The volume-averaged r.m.s. DNS velocity is defined as

(3.2) $$\begin{eqnarray}U_{rms}^{DNS}=\sqrt{\frac{\displaystyle \mathop{\sum }_{j=1}^{n}([U^{DNS}]_{j})^{2}[\unicode[STIX]{x0394}V_{j}]}{V}}.\end{eqnarray}$$

3.1.1 Reynolds stress models with explicit treatment

The Reynolds stress term is directly computed from DNS data and substituted into RANS as shown in algorithm 1. The purpose is to study most existing data-driven RANS models, in which the Reynolds stress is directly predicted by training on DNS data and then used to solve for mean velocity field. In this subsection, we only study the explicit treatment with fixed Reynolds stress, i.e. the dependence of the Reynolds stress on strain-rate tensor is not considered. It is because such a dependence has not been taken into consideration in many data-driven turbulence models, and thus we illustrate the corresponding issue here. More results of explicit treatment with dependence on strain-rate tensor can be found in appendix C, where we further show that the explicit coupling between Reynolds stress and mean velocity during the iterations can gradually amplify the errors.

The local condition numbers ${\mathcal{K}}_{j}$ of the explicit treatment of the Reynolds stress are shown in figure 4. With the increase of the Reynolds number, it can be seen that the magnitude of local condition number also increases. Specifically, the local condition number of the flow at $Re_{\unicode[STIX]{x1D70F}}=180$ is of $O(1)$ , while the local condition number of the flow at $Re_{\unicode[STIX]{x1D70F}}=5200$ is of $O(10^{2})$ . This rapid increase of local condition number agrees well with the increased error of solved mean velocity $U$ as summarized in table 1. In addition, the local condition number is greater near the channel centre than close to the wall, especially for the high-Reynolds-number cases. Such pattern of local condition number also agrees with the spatial pattern of the error of the solved mean velocity as illustrated in figure 9(b). As demonstrated by Thompson et al. (Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016), the error of solved channel flow mean velocity at a given point is an integration of Reynolds stress errors (which happen to be of the same sign in the entire domain, indicating a systematic nature) from the wall to that given point, i.e. the Reynolds stress errors in between are accumulated without cancellation. The local condition number in figure 4 faithfully reflects such an accumulative nature with deteriorated conditioning further away from the wall for all cases. Recall that the local condition number ${\mathcal{K}}_{j}$ assesses the relative error of solved mean velocity at a given point with respect to the error in the whole Reynolds stress field, and not with respect to the error of the Reynolds stress at the same given point.

Figure 4.

The profiles of local condition number ${\mathcal{K}}_{j}$ at different Reynolds numbers by using explicit treatment with fixed Reynolds stress, i.e. the fixed DNS Reynolds stress is substituted into RANS equations.

The averaged local condition number $\overline{{\mathcal{K}}}_{\boldsymbol{x}}$ increases with the Reynolds number by using explicit treatment of the Reynolds stress, which is clearly seen in figure 7. Such increase of averaged local condition number with Reynolds number reveals the potential shortcoming of explicit modelling of the Reynolds stress, i.e. a relatively accurate but explicit modelling of the Reynolds stress does not guarantee the satisfactory mean velocity by solving RANS equations, especially for high-Reynolds-number flows. This observation has been reported in the work of Thompson et al. (Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean2016), and the proposed averaged local condition number can be used as an integral indicator to estimate the extent of error propagation from the modelled Reynolds stress to the solved mean velocity field. The error propagation with iterations is presented in figure 5(a) together with the averaged local condition number in each iteration. It can be seen that the error in the solved mean velocity stays constant in every iteration step. This observation is consistent with our expectation since the convection term disappears in the channel flows and the RANS equations become linear. Therefore, the error within the solved mean velocity in a given iteration does not influence the model conditioning and thus has no effect upon the error in the solved mean velocity at the next iteration.

Figure 5.

The error propagation analysis of solving streamwise velocity iteratively by using explicit treatment with fixed Reynolds stress, including (a) relative error of mean velocity and (b) volume-averaged local condition number.

3.1.2 Reynolds stress models with implicit treatment

The eddy viscosity is directly computed from DNS data and substituted into RANS equations to study the ideal situation of data-driven Reynolds stress models with implicit treatment as shown in algorithm 2.

Compared to the Reynolds stress models, it is well known that eddy-viscosity models can enhance the stability and conditioning of RANS equations with turbulence closures. In the practice of traditional RSM, the modelled Reynolds stress is empirically blended with the Reynolds stress from eddy-viscosity models to achieve better convergence and conditioning (Basara & Jakirlic Reference Basara and Jakirlic2003; Maduta & Jakirlic Reference Maduta and Jakirlic2017). In this work, we demonstrate that the local condition number ${\mathcal{K}}_{j}$ can quantitatively explain the improved conditioning of the implicit treatment of the Reynolds stress by introducing an eddy viscosity. It can be seen in figure 6 that the local condition number ${\mathcal{K}}_{j}$ is significantly reduced compared with the results of explicit treatment of the Reynolds stress as shown in figure 4, especially for high-Reynolds-number flows. Although the local condition number of high-Reynolds-number flow is still greater than that for low Reynolds number, they are of the same order of magnitude for different Reynolds numbers. The volume-averaged local condition number in figure 7 is also significantly reduced by using an implicit treatment of the Reynolds stress, demonstrating the merit of using an implicit treatment of the Reynolds stress in improving the conditioning when solving RANS equations for mean velocity field. The conditioning can be further improved by adjusting $\unicode[STIX]{x1D70F}_{DNS}^{\bot }$ to enhance the implicit treatment part, e.g. increasing $\unicode[STIX]{x1D708}_{t}^{m}$ by $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{t}$ and adjusting the nonlinear part of the Reynolds stress as $\unicode[STIX]{x1D749}_{DNS}^{\bot }-\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{t}(\unicode[STIX]{x1D735}\overline{u}^{(i)}+(\unicode[STIX]{x1D735}\overline{u}^{(i)})^{\text{T}})$ accordingly. However, it should be noted that such a purely numerical enhancement may introduce excessive errors to iterative solvers when the chosen $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{t}$ is too large. The proposed scheme aims to strike a balance between accuracy and conditioning.

Figure 6.

The local condition number at different Reynolds numbers by using implicit treatment of Reynolds stress, i.e. the linear part of Reynolds stress is implicitly treated by introducing an optimal eddy viscosity.

Figure 7.

The volume-averaged local condition number at different Reynolds numbers for explicit and implicit treatments of Reynolds stress.

We further show that the relative error of mean velocity is much smaller by using an implicit treatment of the Reynolds stress in RANS simulations. It can be seen in figure 8(a) that the relative error of the solved mean velocity is much smaller than the one shown in figure 5(a). In addition, the volume-averaged local condition number stays at $O(1)$ as shown in figure 8(b), which explains the better convergence of solving for mean velocity field by using an implicit treatment of the Reynolds stress.

Figure 8.

The error propagation analysis of solving streamwise velocity iteratively by using implicit treatment of Reynolds stress, including (a) relative error of mean velocity and (b) volume-averaged local condition number.

Figure 9.

The comparison of solved mean velocity by using explicit and implicit treatments of Reynolds stress, including (a) mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=180$ , (b) mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=5200$ , (c) percentage error of mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=180$ and (d) percentage error of mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=5200$ .

The mean velocity $U$ is solved and presented in figure 9 at Reynolds numbers $Re_{\unicode[STIX]{x1D70F}}=180$ and $5200$ by using explicit and implicit treatments of the Reynolds stress. At Reynolds number $Re_{\unicode[STIX]{x1D70F}}=180$ , it can be seen in figure 9(a) that the solved mean velocity $U$ by using both kinds of treatments has a good agreement with DNS data. This demonstrates that the error propagation from Reynolds stress to mean velocity is not severe at low Reynolds number, and the percentage error of the mean velocity is comparable by using either an explicit or implicit treatment of the Reynolds stress as shown in figure 9(c). These results show good agreement with the local condition number presented in figures 4 and 6, which shows that the local condition number is of the same order for the flow at Reynolds number $Re_{\unicode[STIX]{x1D70F}}=180$ by using both types of treatments. However, the solved mean velocity fields are noticeably different at high Reynolds number ( $Re_{\unicode[STIX]{x1D70F}}=5200$ ) as shown in figure 9(b). Specifically, the solved mean velocity by using an explicit treatment of the Reynolds stress shows a significant difference from the DNS data, while the solved mean velocity by using an implicit treatment of the Reynolds stress still agrees well with the DNS data at $Re_{\unicode[STIX]{x1D70F}}=5200$ . The subtle differences between the explicit and implicit treatments are further discussed in appendix D. The percentage error of solved mean velocity at $Re_{\unicode[STIX]{x1D70F}}=5200$ in figure 9(d) also confirms that the error in mean velocity by using an explicit treatment of the Reynolds stress is orders of magnitude higher than the error of using an implicit treatment of the Reynolds stress. Such comparison of solved mean velocity fields agrees well with the differences in local condition number  ${\mathcal{K}}_{j}$ , demonstrating that the proposed local condition number can be used to quantitatively assess the error propagation from Reynolds stress to mean velocity when solving RANS equations. At a high Reynolds number ( $Re_{\unicode[STIX]{x1D70F}}=5200$ ), the profile of the error is dominated by the local condition number, which is much larger near the channel centre as shown in figure 4. At a low Reynolds number ( $Re_{\unicode[STIX]{x1D70F}}=180$ ), the local condition number is of the same order of magnitude in the whole domain, and thus the profile of the error is dominated by the error within the DNS Reynolds stress. Therefore, the noisy pattern of the profile of the error can only be seen in figure 9(c) for $Re_{\unicode[STIX]{x1D70F}}=180$ but not in figure 9(d) for $Re_{\unicode[STIX]{x1D70F}}=5200$ . Inspired by the comparison of conditioning with explicit and implicit treatments of Reynolds stress closures, Wu et al. (Reference Wu, Xiao and Paterson2018) further proposed an implicit treatment of the Reynolds stress for machine-learning-assisted RANS modelling to improve the conditioning when solving for mean velocity field.

3.2 Model conditioning of RANS equations for more complex flows

We further study the model conditioning of RANS equations for more complex flows including the flow in a square duct at $Re_{b}=3500$ (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010) and the flow over periodic hills at $Re_{b}=5600$ (Breuer et al. Reference Breuer, Peller, Rapp and Manhart2009), where $Re_{b}$ is defined by the bulk velocity $U_{b}$ at the inlet. The flow configurations are shown in figure 10. In this work, we show that the RANS equations of more complex flows can also be ill-conditioned. In addition, we demonstrate that the proposed local condition number can be used to assess the model conditioning of RANS equations in these more complex flows.

Figure 10.

The configuration of (a) the flow in a square duct at $Re_{b}=3500$ and (b) the flow over periodic hills at $Re_{b}=5600$ .

3.2.1 Flow in a square duct

We first study the solved mean secondary velocity $U_{z}$ by using an explicit treatment with fixed Reynolds stress and an implicit treatment of the Reynolds stress. The comparison of mean velocity profiles demonstrates that the implicit treatment of the Reynolds stress leads to solved mean velocity that shows better agreement with DNS data in figure 11. We then focus on the analysis of the errors in this work. The error is quantified by the ratio $\Vert U_{z}-U_{z}^{DNS}\Vert /U_{z,rms}^{DNS}$ , where $U_{z,rms}^{DNS}$ is a volume-averaged velocity of $U_{z}^{DNS}$ as defined in (3.2). It can be seen in figure 12(a) that some large errors exist in the region of the vertical symmetry plane and around the diagonal within the cross-plane. Compared to the errors as shown in figure 12(a), noticeable reduction of errors can be observed in figure 12(b), where the implicit treatment of the Reynolds stress is used.

Figure 11.

The comparison of solved mean velocities by using explicit treatment with fixed Reynolds stress and implicit treatment of Reynolds stress. The computational domain covers a quarter of the cross-section of the physical domain, i.e.  $h=D/2$ . This is due to the symmetry of the mean flow in both $y$ and $z$ directions as shown in figure 10(a). It should be noted that the Reynolds stress is obtained from a DNS database (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010).

Figure 12.

The normalized error of the solved secondary flow velocity $U_{z}$ of the flow in a square duct by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

The proposed local condition number can be used to analyse the model conditioning for the flow in a square duct as demonstrated in figure 13. It can be seen in figure 13(a) that the local condition number is also large in the region of the vertical symmetry plane and around the diagonal within the cross-plane, which is consistent with the error of mean velocity as shown in figure 12(a). In addition, figure 12(b) shows that the local condition number is generally smaller across the whole domain by using an implicit treatment of the Reynolds stress. Such a reduction of local condition number also correlates well with the comparison of mean velocity error in figure 12. It should be noted that the mean velocity error is determined by both the local condition number and the error in Reynolds stress. Therefore, the spatial pattern of mean velocity error in figure 12 cannot be solely explained by the local condition number in figure 13. However, the analysis of local condition number can still provide some information about whether the solved mean velocity is reliable. In practical applications, the error of the Reynolds stress is usually unknown, and caution should be exercised when regions with large local condition number exist. Note that the results shown in figures 12 and 13 demonstrate that the vertical velocity $U_{z}$ is not symmetric with itself about the diagonal of the domain. Rather, the diagonal symmetry of this flow is such that the velocity components $U_{z}$ and $U_{y}$ are symmetrical to each other about the diagonal. That is, it is expected that $U_{y}(\boldsymbol{x})=U_{z}(\boldsymbol{x}^{\prime })$ and not $U_{z}(\boldsymbol{x})=U_{z}(\boldsymbol{x}^{\prime })$ , where $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime }$ denote two points symmetric about the diagonal.

Figure 13.

The local condition number of the flow in a square duct by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

3.2.2 Flow over periodic hills

Unlike the flow in a square duct, the error of mean velocity can be large across the whole domain for the flow over periodic hills. We first show the comparison of the Reynolds stress from different DNS datasets, including the data by Breuer et al. (Reference Breuer, Peller, Rapp and Manhart2009) (dataset 1) and two datasets obtained by using Incompact3d (Laizet & Lamballais Reference Laizet and Lamballais2009; Laizet & Li Reference Laizet and Li2011) by using two different mesh resolutions (datasets 2 and 3). The simulation details are summarized in table 2. It can be seen in figure 14(a) that the Reynolds stresses obtained from these three datasets are very close to each other. According to the comparison of these Reynolds stresses, we might intuitively expect that the solved mean velocity fields are similar to each other by substituting the fixed DNS Reynolds stresses into RANS equations. However, this is not the case, as shown in figure 14(b), where we compare the solved mean velocity with the DNS mean velocity field. It can be seen that the solved mean velocity field by using fixed DNS Reynolds stress from dataset 2 show noticeable differences across the whole domain. In contrast, the solved mean velocity fields from datasets 1 and 3 have better agreement with the DNS mean velocity field, but they are still different from each other. It is expected that implicit treatment leads to better agreement of solved mean velocity with DNS data, which has been demonstrated in a related work (Wu et al. Reference Wu, Xiao and Paterson2018).

Figure 14.

The comparison of (a) DNS Reynolds stress and (b) solved mean velocity between different DNS datasets. The solved mean velocity is obtained by substituting the corresponding DNS Reynolds stress and solving RANS equations. The solid black lines in panel (b) denote the DNS mean velocity as a benchmark. It should be noted that the DNS mean velocities from different datasets have no noticeable difference, and here we present the DNS mean velocity from Breuer et al. (Reference Breuer, Peller, Rapp and Manhart2009).

Table 2.

Summary of the datasets of the flow over periodic hills at Reynolds number $Re=5600$ , including the numerical methods, treatment of solid boundaries, accuracy of numerical discretization, mesh sizes and type. Here $N_{x}$ , $N_{y}$ , and $N_{z}$ indicate number of grids in streamwise, wall-normal and spanwise directions, respectively.

The comparison in figure 14 can be explained by studying the model conditioning of RANS equations with specified Reynolds stress. Specifically, it can be seen in figure 15(a) that the local condition number by using fixed Reynolds stress is of $O(10^{2})$ in most areas, indicating that the RANS equations are ill-conditioned in these regions. On the other hand, the local condition number is smaller in the near-wall and the recirculation regions. The large local condition number can explain the comparison in figure 14, i.e. the similar Reynolds stress fields lead to dramatically different mean velocity fields by solving the RANS equations. The physical justification of such a pattern of the local condition number field is that the mean velocity error is strongly correlated along the streamlines, and the periodic boundary condition of the inlet and outlet exacerbates the model conditioning for this flow. We also study the local condition number by using an implicit treatment of the Reynolds stress. It can be seen in figure 15(b) that the local condition number is much smaller than those where fixed Reynolds stress is used in figure 15(a). Therefore, it can be expected that the solved mean velocity by using an implicit treatment of DNS Reynolds stress have a better agreement with the DNS mean velocity, which has been confirmed by Wu et al. (Reference Wu, Xiao and Paterson2018).

In this work, the local condition number assesses the relative error of the solved mean velocity at a given point with regard to the errors in the Reynolds stress field. Therefore, the global effect of error can be captured by the local condition number. Specifically, such an effect depends on the mean flow pattern, which is embodied in the differential operator of linearized RANS equations and influences the local condition number via the Green’s function $G(\boldsymbol{x},\unicode[STIX]{x1D743})$ as defined in (2.15). For instance, the errors in the solved mean velocity of the periodic hill flow is generally large across the upper channel region as shown in figure 14. This is largely due to the fact that the cyclic boundary conditions of the inlet and outlet introduce a strong correlation among errors in the upper channel region. It can be seen in figure 15 that the local condition number ${\mathcal{K}}_{j}$ is also large across the upper channel region, demonstrating that the local condition number takes into account the mean flow pattern and thus truthfully reflects the potentially large error in the solved mean velocities therein.

Figure 15.

The local condition number ${\mathcal{K}}_{j}$ of the flow over periodic hills at $Re=5600$ by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

4 Discussion

While this work primarily focuses on the conditioning of a particular class of data-driven Reynolds stress models, the conditioning issue is an equally important challenge for traditional Reynolds stress models based on Reynolds stress transport equations. Although solving a monolithic system of Reynolds stress transport equations and RANS equations is the most effective way to improve conditioning, it is uncommon due to increased computational costs. Moreover, monolithic coupling is by no means a panacea that guarantees well-conditioning and stability. The conditioning and stability ultimately depend on the characteristics of the turbulence model itself. For example, the popularity of the SA model in external aerodynamics is largely attributed to its excellent robustness in terms of both model conditioning and numerical stability, which many other models do not have. As of now, most commercial and open-source general-purpose CFD packages (e.g. Weller et al. Reference Weller, Tabor, Jasak and Fureby1998) still solve the turbulence transport equations and the RANS mean flow equations in a segregated manner, even in solvers where velocity and pressure are solved concurrently. In such segregated solvers, the modelled Reynolds stress is often updated with the mean velocity field at every iteration step, in the hope that the mean velocity field and the Reynolds stress can consistently adjust to each other during the iterations. However, if the RANS equations are ill-conditioned and the Reynolds stresses are treated explicitly as source terms, the error can be amplified within each iteration. Consequently, a small error in the modelled Reynolds stress can lead to large errors in the solved mean velocity field, which is carried over to the Reynolds stresses in the next iteration step and further amplified. Such an error amplification destabilizes the solution procedure and can potentially lead to divergence.

For the reasons outlined above, RANS simulations with Reynolds-stress-based turbulence models need to be stabilized to increase the robustness of the solvers. Examples of stabilization include (1) using the velocity solved with eddy-viscosity models as an initial condition for iterations in the RSM-based solver and (2) partial implicit treatment of the Reynolds stress, among others. In the latter category, researchers introduced a hybrid scheme of computing Reynolds stress by blending the RSM-modelled Reynolds stress with that computed from eddy-viscosity models, with the latter stabilizing the solution (Basara & Jakirlic Reference Basara and Jakirlic2003; Maduta & Jakirlic Reference Maduta and Jakirlic2017). However, the choice of the blending factor is largely ad hoc due to the lack of a quantitative method to evaluate the model conditioning. A large blending factor improves the conditioning and stabilizes the solution of the RANS equations, but it impairs the accuracy of the solved mean velocity, since the linear eddy-viscosity assumption would be increasingly dominant. The metric proposed in this work can assess the model conditioning with any given blending factor, and thus it is possible to choose a minimum blending factor that maintains good conditioning.

A monolithic coupling for data-driven turbulence models is more challenging than for traditional PDE-based models, if it is possible at all. For example, for a neural-network-based data-driven turbulence model (e.g. Ling et al. Reference Ling, Kurzawski and Templeton2016; Zhu et al. Reference Zhu, Zhang, Kou and Liu2019), a monolithic coupling is possible, because neural network models are differentiable. However, for non-differentiable models, e.g. those based on random forests or other tree-based models (e.g. Wang et al. Reference Wang, Wu and Xiao2017), a monolithic coupling is not straightforwardly viable.

5 Conclusion

Recently, several researchers have employed DNS Reynolds stress data as the closure term and solved the RANS equations for mean velocities on turbulent channel flows. They reported unexpected results that the obtained mean velocities deviated significantly (up to 35 %) from the DNS data at high Reynolds numbers. In this work, we aim to identify a metric to quantitatively assess the conditioning of RANS equations with data-driven Reynolds stress closures, i.e. how a small error in Reynolds stress can lead to large errors in the mean velocity by solving RANS equations. The turbulent channel flow is studied to evaluate the candidate metrics. Our analysis shows that the global matrix-based condition number is not able to distinguish the different sensitivity of solved mean velocities at different Reynolds numbers. A local condition number function is then derived as a more precise indicator of model conditioning. We demonstrate that such a local condition number explains the error propagation from the modelled Reynolds stress to the solved mean velocity in RANS simulations for turbulent channel flows at different Reynolds numbers. Two more complex flows are also studied to further demonstrate the capability of the proposed local condition number in evaluating the conditioning of RANS equations with data-driven Reynolds stress closures. The proposed condition number provides a quantitative metric to assess the model conditioning of RANS equations, facilitating the development of conditioning-oriented schemes in data-driven turbulence modelling.

Acknowledgements

The authors would like to thank Dr S. Laizet of Imperial College London and Professor M. Breuer of Helmut-Schmidt-Universität, Hamburg, for providing us with their high-fidelity simulation data of the flow over periodic hills. The authors would also like to thank Dr G. N. Coleman of NASA Langley and Dr F. Menter of ANSYS for helpful discussions during this research. We gratefully acknowledge Dr S. Murman of NASA Ames and Professor S. Poroseva of the University of New Mexico for their constructive comments on the first draft of the manuscript. The authors would also like to thank the reviewers for their constructive and valuable comments, which helped to improve the quality and clarity of this paper.

Appendix A. Derivations of condition numbers

A.1 Derivation of global matrix-based condition number

The global matrix-based condition number ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ is defined as follows:

(A 1) $$\begin{eqnarray}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{U}\Vert }{\Vert \boldsymbol{U}\Vert }\leqslant {\mathcal{K}}_{\unicode[STIX]{x1D70F}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert }{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert },\end{eqnarray}$$

where ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ measures the sensitivity of the solved mean velocity field due to the perturbation of the Reynolds stress field, and $|\cdot |$ indicates the Euclidean norm of a vector (of all values in the discretized velocity field). To derive the formulation of  ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ , the perturbation $\unicode[STIX]{x1D6FF}\boldsymbol{b}$ in (2.8) is further written as

(A 2) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\boldsymbol{b}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}-\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D735}p).\end{eqnarray}$$

For the purpose of the sensitivity study here, it is assumed that a constant pressure gradient is imposed to drive the flow, i.e.  $\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D735}p)=0$ , and thus we have

(A 3) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\boldsymbol{b}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}.\end{eqnarray}$$

Hence,

(A 4) $$\begin{eqnarray}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{U}\Vert }{\Vert \boldsymbol{U}\Vert }\leqslant {\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D6FF}\boldsymbol{b}\Vert }{\Vert \boldsymbol{b}\Vert }={\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert }{\Vert \boldsymbol{b}\Vert }=\underbrace{{\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert }{\Vert \boldsymbol{b}\Vert }}_{{\mathcal{K}}_{\unicode[STIX]{x1D70F}}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert }{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert }.\end{eqnarray}$$

Comparing the forms of (A 4) with the definition of ${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$ in (A 1), the matrix-norm-based condition number for Reynolds-stress-based turbulence models is thus

(A 5) $$\begin{eqnarray}{\mathcal{K}}_{\unicode[STIX]{x1D70F}}={\mathcal{K}}_{\unicode[STIX]{x1D63C}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert }{\Vert \boldsymbol{b}\Vert }.\end{eqnarray}$$

A.2 Derivation of local condition number function

The continuous local condition number ${\mathcal{K}}(\boldsymbol{x})$ is defined as follows:

(A 6) $$\begin{eqnarray}\frac{|\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x})|}{U_{\infty }}\leqslant {\mathcal{K}}(\boldsymbol{x})\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}},\end{eqnarray}$$

where ${\mathcal{K}}(\boldsymbol{x})$ measures the sensitivity of the solved mean velocity at any given location $\boldsymbol{x}$ due to the perturbation of the Reynolds stress field, and $U_{\infty }$ is a constant representative velocity magnitude for normalization. The function norm $\Vert \cdot \Vert _{\unicode[STIX]{x1D6FA}}$ of function $f(\unicode[STIX]{x1D743})$ on domain $\unicode[STIX]{x1D6FA}$ is defined as in (2.17).

To derive the formulation of this local condition number, we first consider the solution $\boldsymbol{u}$ at a particular location  $\boldsymbol{x}^{\prime }$ :

(A 7) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x}^{\prime })=\int _{\unicode[STIX]{x1D6FA}}G(\boldsymbol{x}^{\prime };\unicode[STIX]{x1D743})\boldsymbol{b}(\unicode[STIX]{x1D743})\,\text{d}\unicode[STIX]{x1D743},\end{eqnarray}$$

where $G$ represents the Green’s function of the linear differential operator ${\mathcal{L}}$ in the linearized RANS equations as defined in (2.6). Denoting $G_{\boldsymbol{x}^{\prime }}=G(\boldsymbol{x}^{\prime };\unicode[STIX]{x1D743})$ , the perturbation of the solution is thus

(A 8) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x}^{\prime })=\int _{\unicode[STIX]{x1D6FA}}G(\boldsymbol{x}^{\prime };\unicode[STIX]{x1D743})\,\unicode[STIX]{x1D6FF}\boldsymbol{b}(\unicode[STIX]{x1D743})\,\text{d}\unicode[STIX]{x1D743}=\langle G_{\boldsymbol{ x}^{\prime }},\unicode[STIX]{x1D6FF}\boldsymbol{b}\rangle _{\unicode[STIX]{x1D6FA}},\end{eqnarray}$$

where $\langle \cdot \rangle _{\unicode[STIX]{x1D6FA}}$ is the inner product of functions defined on domain  $\unicode[STIX]{x1D6FA}$ .

Using the Schwartz inequality (Steele Reference Steele2004; Debnath & Mikusiński Reference Debnath and Mikusiński2005) leads to

(A 9) $$\begin{eqnarray}\displaystyle |\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x}^{\prime })| & {\leqslant} & \displaystyle \Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D6FF}\boldsymbol{b}\Vert _{\unicode[STIX]{x1D6FA}}\end{eqnarray}$$
(A 10) $$\begin{eqnarray}\displaystyle & = & \displaystyle \Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}.\end{eqnarray}$$

As in § A.1, the pressure gradient is assumed constant and thus $\unicode[STIX]{x1D6FF}\boldsymbol{b}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}$ . Finally, the sensitivity of mean velocity $\boldsymbol{u}$ with respect to the Reynolds stress $\unicode[STIX]{x1D749}$ perturbations is derived as follows:

(A 11) $$\begin{eqnarray}\displaystyle \frac{|\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x}^{\prime })|}{U_{\infty }} & {\leqslant} & \displaystyle \frac{\Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}\end{eqnarray}$$
(A 12) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}\,\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}.\end{eqnarray}$$

Therefore, by comparing (A 12) and (A 6), we define a local condition number function ${\mathcal{K}}$ of spatial location $\boldsymbol{x}$ as

(A 13) $$\begin{eqnarray}{\mathcal{K}}(\boldsymbol{x})=\frac{\Vert G_{\boldsymbol{x}}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}=\frac{\Vert G(\boldsymbol{x},\unicode[STIX]{x1D743})\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}.\end{eqnarray}$$

Without causing ambiguity, we have dropped the subscript of $\boldsymbol{x}^{\prime }$ in the equation above and in the text for simplicity of notation.

A.3 Local condition number for implicit treatment of Reynolds stress

In the practice of RANS modelling, eddy-viscosity models are widely used, and the modelled eddy viscosity influences the differential operator ${\mathcal{L}}$ associated with RANS equations. Therefore, we extend the derivation of (2.15) to make it compatible with the implicit treatment of the Reynolds stress. According to the general form of the implicit treatment of the Reynolds stress (Pope Reference Pope1975) in (2.22), the linearized RANS equations in (2.5) can be rearranged as follows:

(A 14) $$\begin{eqnarray}\widetilde{{\mathcal{L}}}(\boldsymbol{u})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{\bot }-\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where $\widetilde{{\mathcal{L}}}={\mathcal{L}}-\unicode[STIX]{x1D708}_{t}^{m}\unicode[STIX]{x1D6FB}^{2}$ is the modified linear differential operator by using an implicit treatment of the Reynolds stress. Examples of optimal eddy viscosity $\unicode[STIX]{x1D708}_{t}^{m}$ for the flow over periodic hills and the flow in a square duct are shown in figure 16. Here we only study the perturbation on the nonlinear term $\unicode[STIX]{x1D749}^{\bot }$ of the Reynolds stress  $\unicode[STIX]{x1D749}$ , i.e. 

(A 15) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}^{\bot }.\end{eqnarray}$$

Finally, we have the local condition number ${\mathcal{K}}(\boldsymbol{x}^{\prime })$ in (2.15) rederived as follows for the implicit treatment of the Reynolds stress:

(A 16) $$\begin{eqnarray}\displaystyle \frac{|\unicode[STIX]{x1D6FF}\boldsymbol{u}(\boldsymbol{x}^{\prime })|}{U_{\infty }} & {\leqslant} & \displaystyle \frac{\Vert \widetilde{G}_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}^{\bot }\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}\end{eqnarray}$$
(A 17) $$\begin{eqnarray}\displaystyle & = & \displaystyle \underbrace{\frac{\Vert \widetilde{G}_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}\,\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{U_{\infty }}}_{{\mathcal{K}}_{j}}\frac{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}}{\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}},\end{eqnarray}$$

where the kernel function $\Vert \widetilde{G}_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}$ represents the Green’s function that corresponds to the linear differential operator $\widetilde{{\mathcal{L}}}$ defined in (2.23), taking into account the implicit modelling of the linear part of the Reynolds stress by introducing an eddy viscosity.

Figure 16.

The optimal eddy-viscosity fields for (a) the flow in a square duct at $Re_{b}=3500$ and (b) the flow over periodic hills at $Re_{b}=5600$ .

Appendix B. Mesh independence of the local condition number function

We present the numerical discretization of the proposed local condition number on a CFD mesh and show that the discretized local condition number is mesh-independent. First, the function norms $\Vert \cdot \Vert _{\unicode[STIX]{x1D6FA}}$ are approximated in vector norms $\Vert \cdot \Vert _{n}$ through numerical integration on a CFD mesh of $n$ cells. This is derived as follows:

(B 1) $$\begin{eqnarray}\displaystyle \Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}} & = & \displaystyle \sqrt{\int _{\unicode[STIX]{x1D6FA}}|G(\boldsymbol{x}^{\prime };\unicode[STIX]{x1D743})|^{2}\,\text{d}\unicode[STIX]{x1D743}}\end{eqnarray}$$
(B 2) $$\begin{eqnarray}\displaystyle & {\approx} & \displaystyle \sqrt{\mathop{\sum }_{i=1}^{n}([\boldsymbol{r}_{j,i}]^{2}\,\unicode[STIX]{x0394}V_{i})}\end{eqnarray}$$
(B 3) $$\begin{eqnarray}\displaystyle & = & \displaystyle \sqrt{\mathop{\sum }_{i=1}^{n}\left([\boldsymbol{r}_{j,i}]\,\sqrt{[\unicode[STIX]{x0394}V_{i}]}\right)^{2}}\end{eqnarray}$$
(B 4) $$\begin{eqnarray}\displaystyle & = & \displaystyle \Vert \boldsymbol{r}_{j}\sqrt{[\unicode[STIX]{x0394}\boldsymbol{V}]}\Vert _{n}=\Vert \tilde{\boldsymbol{r}}_{j}\Vert _{n},\end{eqnarray}$$

with $\tilde{\boldsymbol{r}}_{j}=\boldsymbol{r}_{j}\sqrt{\unicode[STIX]{x0394}\boldsymbol{V}}$ and ${\approx}$ indicating numerical discretization of the integral involved in the function norm in (B 1); $\unicode[STIX]{x0394}V_{i}$ is the volume for the $i$ th cell; $\unicode[STIX]{x0394}\boldsymbol{V}$ is the $n$ -vector consisting of volumes of cells in the mesh; and $\boldsymbol{r}_{j}$ is the $j$ th row of the matrix $\unicode[STIX]{x1D63C}^{-1}$ , with $\unicode[STIX]{x1D63C}$ being the discretization of the operator ${\mathcal{L}}$ as seen in (2.6). The numbering implies that the location $\boldsymbol{x}^{\prime }$ is the coordinate of the $j$ th cell.

Similarly, the function norm of the Reynolds stress divergence is approximated on the CFD mesh as follows:

(B 5) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}\approx \Vert [\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]\sqrt{[\unicode[STIX]{x0394}\boldsymbol{V}]}\Vert _{n}.\end{eqnarray}$$

It is clear that the function norm $\Vert G_{\boldsymbol{x}^{\prime }}\Vert _{\unicode[STIX]{x1D6FA}}$ is mesh-independent as its definition does not involve any discretization mesh (quadrature), so its numerical approximation $\Vert \boldsymbol{r}_{j}\sqrt{[\unicode[STIX]{x0394}\boldsymbol{V}]}\Vert _{n}$ should also be mesh-independent on a sufficiently fine mesh. In the same way, both the function norm $\Vert \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\Vert _{\unicode[STIX]{x1D6FA}}$ and its numerical approximation $\Vert [\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]\sqrt{[\unicode[STIX]{x0394}\boldsymbol{V}]}\Vert _{n}$ are mesh-independent. The mesh independence can be further verified in the special case, where the divergence field $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}$ is a non-zero constant $\unicode[STIX]{x1D6FD}$ and the mesh consists of $n$ uniformly sized cells. In this case we have $\unicode[STIX]{x0394}V=|\unicode[STIX]{x1D6FA}|/n$ , where $|\unicode[STIX]{x1D6FA}|$ is the total volume of the computational domain $\unicode[STIX]{x1D6FA}$ (independent of the discretization mesh). Therefore, the vector norm, which is a numerical approximation of the function form, is as follows:

(B 6) $$\begin{eqnarray}\displaystyle \Vert [\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]\sqrt{[\unicode[STIX]{x0394}\boldsymbol{V}]}\Vert _{n} & = & \displaystyle \sqrt{\mathop{\sum }_{i=1}^{n}\left(\unicode[STIX]{x1D6FD}\sqrt{\frac{|\unicode[STIX]{x1D6FA}|}{n}}\right)^{2}}\end{eqnarray}$$
(B 7) $$\begin{eqnarray}\displaystyle & = & \displaystyle \sqrt{n\unicode[STIX]{x1D6FD}^{2}\frac{|\unicode[STIX]{x1D6FA}|}{n}}\end{eqnarray}$$
(B 8) $$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FD}\sqrt{|\unicode[STIX]{x1D6FA}|},\end{eqnarray}$$

which is clearly independent of the number of cells in the mesh. In order to complement and validate the derivations above, a numerical study of mesh convergence is presented in figure 17, where the local condition number ${\mathcal{K}}_{j}$ at $Re=5200$ is calculated by using different mesh resolutions. It can be seen that asymptotical convergence of local condition number ${\mathcal{K}}_{j}$ can be achieved by gradually refining the mesh resolution. Moreover, even on a coarser mesh the overall magnitude and spatial distribution of the condition number ${\mathcal{K}}_{j}$ are correctly captured.

Figure 17.

Mesh convergence study of local condition number ${\mathcal{K}}_{j}$ at $Re=5200$ by using explicit treatment with fixed Reynolds stress, i.e. the fixed DNS Reynolds stress is substituted into RANS equations.

Appendix C. Explicit treatment of Reynolds stress with dependence on strain rate

The conditioning of RANS equations with explicit treatment of the Reynolds stress has been discussed in § 3.1.1, where the Reynolds stress is obtained from a DNS database and is fixed when solving RANS equations. If the dependence of the Reynolds stress on strain rate is considered, it can be seen in figure 18 that the conditioning of RANS equations in the first iteration is the same as the conditioning with fixed Reynolds stress as shown in figure 5. This is because the detailed explicit treatment of the Reynolds stress does not influence the model conditioning at a given state. Here we further show that the explicit coupling between Reynolds stress and mean velocity during the iterations can gradually amplify the errors and lead to divergence. Such explicit coupling is often used in the Reynolds stress transport models. Specifically, the Reynolds stress is obtained by solving its transport equations with the mean velocity field at the previous iteration step. In the following, we use a simplified example with an iterative solver as shown in algorithm 3 to illustrate the convergence issue of Reynolds stress transport models. In addition, we demonstrate that the proposed local condition number can be used to detect and explain the corresponding ill-conditioning issue during iterations.

Figure 18.

The error propagation analysis of solving streamwise velocity iteratively by using explicit treatment of Reynolds stress that depends on the strain rate, including (a) relative error of mean velocity and (b) volume-averaged local condition number.

The Reynolds stress at the $i\text{th}$ iteration step is explicitly treated by using DNS data according to algorithm 3. Unlike the data-driven Reynolds stress modelling as shown in algorithm 1, this simplified explicit Reynolds stress treatment allows the update of the Reynolds stress at each iteration based on the solved mean velocity field at the previous iteration step. Compared to the implicit treatment of the Reynolds stress as shown in algorithm 2, the only difference of this explicit treatment is the computing of the Reynolds stress with the mean velocity at the previous iteration step, which is indicated by the superscript $i-1$ at line 3 of algorithm 3.

The errors of solved mean velocity field by using an explicit treatment of the Reynolds stress is presented in figure 18(a). The DNS mean velocity is used as the initial field in RANS simulations, and thus the initial value of $\unicode[STIX]{x1D6FF}U_{rms}/U_{rms}^{DNS}$ is small. However, the value of $\unicode[STIX]{x1D6FF}U_{rms}/U_{rms}^{DNS}$ increases rapidly during the first several iteration steps. It demonstrates that the conditioning issue within each iteration can lead to error amplification, i.e. a small error in the modelled Reynolds stress can lead to noticeable errors in the solved mean velocity field and thus influence the modelled Reynolds stress in the next iteration step. Owing to such coupling of error amplification, even a small error of modelled Reynolds stress can lead to divergence of the simulation eventually. It can be seen in figure 18(b) that the volume-averaged local condition number is of $O(10^{2})$ within the first three iteration steps, explaining the rapid growth of error in the solved mean velocity. The error of the solved mean velocity grows rapidly and eventually leads to divergence of the simulation. Therefore, the solved mean velocity is not presented in this work since a converged solution was not achieved. Therefore, the proposed local condition number is still of importance in traditional RANS modelling since it provides a quantitative assessment of model conditioning at every iteration step. The divergence of the mean velocity is because the Reynolds stress is updated according to the mean velocity at each iteration step. With the rapidly increased error in the mean velocity at the first several iteration steps, the error in the Reynolds stress also increases accordingly. The decrease of the condition number as shown in figure 18(b) does not guarantee the decrease of the error in the mean velocity in such a scenario, considering that the error of mean velocity is also influenced by the error in the modelled Reynolds stress. Therefore, the small condition number needs to be interpreted with caution when the source term in RANS equations changes during the simulation. With large condition number in some previous iteration steps, it is possible that the error of the modelled Reynolds stress becomes large and the mean velocity error keeps increasing even though the condition number decreases.

Appendix D. Discrepancies in velocities obtained with different treatments

It was shown in figure 9 that the solved mean velocity can be significantly different depending on whether the DNS Reynolds stress used to solve (2.4) is treated explicitly or implicitly. In other words, solving the following two equations

(D 1) $$\begin{eqnarray}\displaystyle & {\mathcal{L}}(\boldsymbol{u}^{exp})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{exp}-\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$
(D 2) $$\begin{eqnarray}\displaystyle & {\mathcal{L}}(\boldsymbol{u}^{imp})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{imp}-\unicode[STIX]{x1D735}p & \displaystyle\end{eqnarray}$$

yields drastically different velocities. The superscript ‘imp’ indicates the implicit treatment of the Reynolds stress and ‘exp’ denotes the explicit treatment, i.e.  $\unicode[STIX]{x1D749}^{imp}=\unicode[STIX]{x1D708}_{t}(\unicode[STIX]{x1D735}\boldsymbol{u}^{imp}+(\unicode[STIX]{x1D735}\boldsymbol{u}^{imp})^{\text{T}})+\unicode[STIX]{x1D749}^{\bot ,DNS}$ and $\unicode[STIX]{x1D749}^{exp}=\unicode[STIX]{x1D749}^{DNS}$ . This finding apparently contradicts the common understanding in traditional CFD practice that the converged solution of the mean velocity should be the same regardless of how the Reynolds stress is treated. Indeed, the Reynolds stresses used in the two formulations in (D 1) and (D 2) are approximately equal, since $\unicode[STIX]{x1D708}_{t}(\unicode[STIX]{x1D735}\boldsymbol{u}^{imp}+(\unicode[STIX]{x1D735}\boldsymbol{u}^{imp})^{\text{T}})+\unicode[STIX]{x1D749}^{\bot ,DNS}\approx \unicode[STIX]{x1D749}^{DNS}$ . More precisely, the difference between $\boldsymbol{u}^{imp}$ and $\boldsymbol{u}^{DNS}$ is approximately 0.1 %, and the difference between $\unicode[STIX]{x1D749}^{imp}$ and $\unicode[STIX]{x1D749}^{exp}$ should be at a similar level. However, the condition number with regard to the nonlinear differential operator ${\mathcal{L}}$ is large for the flows at high Reynolds numbers, and thus a small difference between $\unicode[STIX]{x1D749}^{imp}$ and $\unicode[STIX]{x1D749}^{exp}$ can lead to a large difference between the solved mean velocities $\boldsymbol{u}^{imp}$ and  $\boldsymbol{u}^{exp}$ .

In addition, the better solution of $\boldsymbol{u}^{imp}$ with the implicit treatment of the Reynolds stress can be explained by the improved model conditioning, i.e. the condition number is smaller with regard to the linear differential operator $\widetilde{{\mathcal{L}}}={\mathcal{L}}-\unicode[STIX]{x1D708}_{t}^{m}\unicode[STIX]{x1D6FB}^{2}$ for the implicit treatment of the Reynolds stress. Specifically, we first define an optimal Reynolds stress $\unicode[STIX]{x1D749}^{op}$ that can lead to $\boldsymbol{u}^{DNS}$ by solving the RANS equations

(D 3) $$\begin{eqnarray}{\mathcal{L}}(\boldsymbol{u}^{DNS})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{op}-\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where $\unicode[STIX]{x1D749}^{op}$ denotes the true Reynolds stress that provides $\boldsymbol{u}^{DNS}$ by solving the RANS equations. The errors $\Vert \unicode[STIX]{x1D749}^{DNS}-\unicode[STIX]{x1D749}^{op}\Vert$ and $\Vert \unicode[STIX]{x1D749}^{\bot ,DNS}-\unicode[STIX]{x1D749}^{\bot ,op}\Vert$ are of the same order of magnitude. Therefore, $\Vert \boldsymbol{u}^{imp}-\boldsymbol{u}^{DNS}\Vert$ is smaller than $\Vert \boldsymbol{u}^{exp}-\boldsymbol{u}^{DNS}\Vert$ due to the smaller sensitivity of solving the mean velocity by using the implicit treatment of the Reynolds stress.

Notation

We use $[\unicode[STIX]{x1D719}]$ to indicate the $n$ -vector obtained by discretizing the field $\unicode[STIX]{x1D719}$ on the mesh, where $n$ is the number of cells/grid points; $\Vert [\unicode[STIX]{x1D719}]\Vert$ denotes the norm of vector $[\unicode[STIX]{x1D719}]$ resulting from the discretization. Since the norm is always taken on the discretized $n$ -vector, we abbreviated $\Vert [\unicode[STIX]{x1D719}]\Vert$ as $\Vert \unicode[STIX]{x1D719}\Vert$ without ambiguity. When mentioning function norm and $n$ -vector norm simultaneously, we used $\Vert \cdot \Vert _{\unicode[STIX]{x1D6FA}}$ and $\Vert \cdot \Vert _{n}$ to distinguish them, with $\unicode[STIX]{x1D6FA}$ denoting the domain on which the norm is defined.

$\boldsymbol{u}$

mean velocity field

$\boldsymbol{U}$

discretized mean velocity ( $n$ -vector)

$\unicode[STIX]{x1D749}$

Reynolds stress tensor

$\unicode[STIX]{x1D64E}$

rate-of-strain tensor

$\boldsymbol{b}$

imbalance vector between Reynolds stress divergence and pressure gradient

$\unicode[STIX]{x1D63C}$

$n\times n$ coefficients matrix of discretized RANS equations

${\mathcal{N}}$

nonlinear differential operator

${\mathcal{L}}$

linear differential operator

$G$

Green’s function corresponding to ${\mathcal{L}}$

${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$

condition number of matrix $\unicode[STIX]{x1D63C}$

$\overline{\unicode[STIX]{x1D6FC}}$

ratio between Reynolds stress divergence and the total source term

${\mathcal{K}}_{\unicode[STIX]{x1D70F}}$

matrix-norm condition number associated with Reynolds stress perturbation

${\mathcal{K}}_{j}$

local condition number vector approximated on a CFD mesh ( $n$ -vector)

$\overline{{\mathcal{K}}}_{x}$

volume-averaged local condition number (scalar)

$\unicode[STIX]{x1D6FF}U_{rms}$

volume-averaged r.m.s. error of solved mean velocity

$U_{rms}^{DNS}$

volume-averaged r.m.s. DNS mean velocity

$\bot$

superscript indicating the nonlinear part of Reynolds stress

References

Basara, B. & Jakirlic, S. 2003 A new hybrid turbulence modelling strategy for industrial CFD. Intl J. Numer. Meth. Fluids 42 (1), 89–116.Google Scholar
Bernardini, M., Pirozzoli, S. & Orlandi, P. 2014 Velocity statistics in turbulent channel flow up to Re = 4000. J. Fluid Mech. 742, 171–191.Google Scholar
Breuer, M., Peller, N., Rapp, C. & Manhart, M. 2009 Flow over periodic hills: numerical and experimental study in a wide range of Reynolds numbers. Comput. Fluids 38 (2), 433–457.Google Scholar
Chandrasekaran, S. & Ipsen, I. C. F. 1995 On the sensitivity of solution components in linear systems of equations. SIAM J. Matrix Anal. Applics. 16 (1), 93–112.Google Scholar
Debnath, L. & Mikusiński, P. 2005 Hilbert Spaces with Applications. Academic Press.Google Scholar
Del Alamo, J. C. & Jiménez, J. 2003 Spectra of the very large anisotropic scales in turbulent channels. Phys. Fluids 15 (6), L41–L44.Google Scholar
Gamahara, M. & Hattori, Y. 2017 Searching for turbulence models by artificial neural network. Phys. Rev. Fluids 2 (5), 054604.Google Scholar
Geneva, N. & Zabaras, N. 2019 Quantifying model form uncertainty in Reynolds-averaged turbulence models with Bayesian deep neural networks. J. Comput. Phys. 383, 125–147.Google Scholar
Hamlington, P. E. & Dahm, W. J. A. 2008 Reynolds stress closure for nonequilibrium effects in turbulent flows. Phys. Fluids 20 (11), 115101.Google Scholar
Hamlington, P. E. & Ihme, M. 2014 Modeling of non-equilibrium homogeneous turbulence in rapidly compressed flows. Flow Turbul. Combust. 93 (1), 93–124.Google Scholar
King, R. N., Hamlington, P. E. & Dahm, W. J. A. 2016 Autonomic closure for turbulence simulations. Phys. Rev. E 93 (3), 031301.Google Scholar
Laizet, S. & Lamballais, E. 2009 High-order compact schemes for incompressible flows: a simple and efficient method with quasi-spectral accuracy. J. Comput. Phys. 228 (16), 5989–6015.Google Scholar
Laizet, S. & Li, N. 2011 Incompact3d: a powerful tool to tackle turbulence problems with up to O (105) computational cores. Intl J. Numer. Meth. Fluids 67 (11), 1735–1757.Google Scholar
Lanczos, C. 1996 Linear Differential Operators. SIAM.Google Scholar
Launder, B. E. & Sharma, B. I. 1974 Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transfer 1 (2), 131–137.Google Scholar
Lee, M. & Moser, R. D. 2015 Direct numerical simulation of turbulent channel flow up to Re = 5200. J. Fluid Mech. 774, 395–415.Google Scholar
Ling, J., Kurzawski, A. & Templeton, J. 2016 Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. J. Fluid Mech. 807, 155–166.Google Scholar
Maduta, R. & Jakirlic, S. 2017 Improved RANS computations of flow over the 25° -slant-angle Ahmed body. SAE Intl J. Passenger Cars – Mech. Syst. 10 (2), 649–661.Google Scholar
Maulik, R. & San, O. 2017 A neural network approach for the blind deconvolution of turbulent flows. J. Fluid Mech. 831, 151–181.Google Scholar
Maulik, R., San, O., Rasheed, A. & Vedula, P. 2019 Sub-grid modelling for two-dimensional turbulence using neural networks. J. Fluid Mech. 858, 122–144.Google Scholar
Menter, F. R. 1994 Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32 (8), 1598–1605.Google Scholar
Parish, E. J. & Duraisamy, K. 2016 A paradigm for data-driven predictive modeling using field inversion and machine learning. J. Comput. Phys. 305, 758–774.Google Scholar
Pinelli, A., Uhlmann, M., Sekimoto, A. & Kawahara, G. 2010 Reynolds number dependence of mean flow structure in square duct turbulence. J. Fluid Mech. 644, 107–122.Google Scholar
Pope, S. B. 1975 A more general effective-viscosity hypothesis. J. Fluid Mech. 72 (2), 331–340.Google Scholar
Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.Google Scholar
Poroseva, S. V., Colmenares, F., Juan, D. & Murman, S. M. 2016 On the accuracy of RANS simulations with DNS data. Phys. Fluids 28 (11), 115102.Google Scholar
Singh, A. P. & Duraisamy, K. 2016 Using field inversion to quantify functional errors in turbulence closures. Phys. Fluids 28, 045110.Google Scholar
Singh, A. P., Medida, S. & Duraisamy, K. 2017 Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils. AIAA J. 55 (7), 2215–2227.Google Scholar
Slotnick, J., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W., Lurie, E. & Mavriplis, D.2014 CFD Vision 2030 Study: a path to revolutionary computational aerosciences. Tech. Rep. National Aeronautics and Space Administration, Langley Research Center, Hampton, VA.Google Scholar
Spalart, P. R. & Allmaras, S. R. 1994 A one-equation turbulence model for aerodynamic flows. Rech. Aerosp, 1, 5–21.Google Scholar
Speziale, C. G. & Xu, X.-H. 1996 Towards the development of second-order closure models for nonequilibrium turbulent flows. Intl J. Heat Fluid Flow 17 (3), 238–244.Google Scholar
Steele, J. M. 2004 The Cauchy–Schwarz Master Class: An Introduction to the Art of Mathematical Inequalities. Cambridge University Press.Google Scholar
Strang, G. 2016 Introduction to Linear Algebra, 5th edn. Wellesley–Cambridge Press.Google Scholar
Thompson, R. L., Sampaio, L. E. B., de Bragança Alves, F. A. V., Thais, L. & Mompean, G. 2016 A methodology to evaluate statistical errors in DNS data of plane channel flows. Comput. Fluids 130, 1–7.Google Scholar
Vollant, A., Balarac, G. & Corre, C. 2017 Subgrid-scale scalar flux modelling based on optimal estimation theory and machine-learning procedures. J. Turbul. 18 (9), 854–878.Google Scholar
Wang, J.-X., Wu, J.-L. & Xiao, H. 2017 Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data. Phys. Rev. Fluids 2 (3), 034603.Google Scholar
Weatheritt, J. & Sandberg, R. 2016 A novel evolutionary algorithm applied to algebraic modifications of the RANS stress–strain relationship. J. Comput. Phys. 325, 22–37.Google Scholar
Weller, H. G., Tabor, G., Jasak, H. & Fureby, C. 1998 A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 12 (6), 620–631.Google Scholar
Wilcox, D. C. 1988 Reassessment of the scale-determining equation for advanced turbulence models. AIAA J. 26 (11), 1299–1310.Google Scholar
Wu, J.-L., Xiao, H. & Paterson, E. 2018 Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework. Phys. Rev. Fluids 3 (7), 074602.Google Scholar
Zhu, L., Zhang, W., Kou, J. & Liu, Y. 2019 Machine learning methods for turbulence modeling in subsonic flows around airfoils. Phys. Fluids 31 (1), 015105.Google Scholar
Figure 0

Table 1.Summary of the results of the channel flow test, showing percentage errors in the turbulent shear stresses and the propagated mean velocities. The large errors in the high-Reynolds-number case $Re_{\unicode[STIX]{x1D70F}}=5200$ are highlighted in bold. The true Reynolds shear stresses are obtained by analytically integrating the RANS equation with the DNS mean velocities (Thompson et al.2016).

Figure 1

Figure 1.Conditioning measure of Reynolds-stress-based turbulence models based on ${\mathcal{K}}_{\unicode[STIX]{x1D63C}}$ and the ratio $\overline{\unicode[STIX]{x1D6FC}}$ as defined in (2.12). The Reynolds stress is computed from DNS data to study the ideal scenario of the RANS modelling.

Figure 2

Figure 2.The balance among forces due to turbulent shear stress ($\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}$), viscous shear stress ($\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{vis}$) and pressure gradient ($\unicode[STIX]{x1D735}p$) for two plane channel flows at frictional Reynolds numbers (a) $Re_{\unicode[STIX]{x1D70F}}=180$ and (b) $Re_{\unicode[STIX]{x1D70F}}=5200$. The right vertical axis denotes the inner coordinates ($y^{+}$).

Figure 3

Figure 3.Force balance of the plane channel flow in (a) the outer layer and (b) the viscous wall region.

Figure 4

Figure 4.The profiles of local condition number ${\mathcal{K}}_{j}$ at different Reynolds numbers by using explicit treatment with fixed Reynolds stress, i.e. the fixed DNS Reynolds stress is substituted into RANS equations.

Figure 5

Figure 5.The error propagation analysis of solving streamwise velocity iteratively by using explicit treatment with fixed Reynolds stress, including (a) relative error of mean velocity and (b) volume-averaged local condition number.

Figure 6

Figure 6.The local condition number at different Reynolds numbers by using implicit treatment of Reynolds stress, i.e. the linear part of Reynolds stress is implicitly treated by introducing an optimal eddy viscosity.

Figure 7

Figure 7.The volume-averaged local condition number at different Reynolds numbers for explicit and implicit treatments of Reynolds stress.

Figure 8

Figure 8.The error propagation analysis of solving streamwise velocity iteratively by using implicit treatment of Reynolds stress, including (a) relative error of mean velocity and (b) volume-averaged local condition number.

Figure 9

Figure 9.The comparison of solved mean velocity by using explicit and implicit treatments of Reynolds stress, including (a) mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=180$, (b) mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=5200$, (c) percentage error of mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=180$ and (d) percentage error of mean velocity $U$ at $Re_{\unicode[STIX]{x1D70F}}=5200$.

Figure 10

Figure 10.The configuration of (a) the flow in a square duct at $Re_{b}=3500$ and (b) the flow over periodic hills at $Re_{b}=5600$.

Figure 11

Figure 11.The comparison of solved mean velocities by using explicit treatment with fixed Reynolds stress and implicit treatment of Reynolds stress. The computational domain covers a quarter of the cross-section of the physical domain, i.e. $h=D/2$. This is due to the symmetry of the mean flow in both $y$ and $z$ directions as shown in figure 10(a). It should be noted that the Reynolds stress is obtained from a DNS database (Pinelli et al.2010).

Figure 12

Figure 12.The normalized error of the solved secondary flow velocity $U_{z}$ of the flow in a square duct by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

Figure 13

Figure 13.The local condition number of the flow in a square duct by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

Figure 14

Figure 14.The comparison of (a) DNS Reynolds stress and (b) solved mean velocity between different DNS datasets. The solved mean velocity is obtained by substituting the corresponding DNS Reynolds stress and solving RANS equations. The solid black lines in panel (b) denote the DNS mean velocity as a benchmark. It should be noted that the DNS mean velocities from different datasets have no noticeable difference, and here we present the DNS mean velocity from Breuer et al. (2009).

Figure 15

Table 2.Summary of the datasets of the flow over periodic hills at Reynolds number $Re=5600$, including the numerical methods, treatment of solid boundaries, accuracy of numerical discretization, mesh sizes and type. Here $N_{x}$, $N_{y}$, and $N_{z}$ indicate number of grids in streamwise, wall-normal and spanwise directions, respectively.

Figure 16

Figure 15.The local condition number ${\mathcal{K}}_{j}$ of the flow over periodic hills at $Re=5600$ by using (a) explicit treatment with fixed Reynolds stress and (b) implicit treatment of Reynolds stress.

Figure 17

Figure 16.The optimal eddy-viscosity fields for (a) the flow in a square duct at $Re_{b}=3500$ and (b) the flow over periodic hills at $Re_{b}=5600$.

Figure 18

Figure 17.Mesh convergence study of local condition number ${\mathcal{K}}_{j}$ at $Re=5200$ by using explicit treatment with fixed Reynolds stress, i.e. the fixed DNS Reynolds stress is substituted into RANS equations.

Figure 19

Figure 18.The error propagation analysis of solving streamwise velocity iteratively by using explicit treatment of Reynolds stress that depends on the strain rate, including (a) relative error of mean velocity and (b) volume-averaged local condition number.