2.1 ML-RANS framework

The ML-RANS framework consists of two independent phases: training phase and prediction. Throughout the ML-RANS modeling process, neural networks, input features and target variables play integral roles. Liu et al. [Reference Liu, Fang, Rolfo, Moulinec and Emerson29] employed the TensorFlow platform [Reference Abadi, Barham, Chen, Chen, Davis, Dean, Devin, Ghemawat, Irving, Isard, Kudlur, Levenberg, Monga, Moore, Murray, Steiner, Tucker, Vasudevan, Warden, Wicke, Yu and Zheng53] to construct an artificial neural network. In the training phase, M2 method described in 1 is first adopted for the construction of machine learning dataset. The machine learning regression algorithms, along with prescribed input features and target variables, are used to determine the weights and biases of the neural network, resulting in a machine learning model. During the prediction phase, this machine learning model is coupled with a baseline turbulence model. It first computes the input features using numerical solutions of both transport equations in the baseline turbulence model and RANS equations. The comparison of turbulent quantities ( $k$ , $\omega $ , …) obtained with DNS/LES with those of RANS formulation can refer to previous studies, for example Figs. 1 and 2 in Ref. [Reference Fang, Bao, Xu, Zhou, Du and Jin35]. Subsequently, these input features are mapped onto the target variable to generate turbulence viscosity. This framework exhibits two notable attributes: consistency in the computation of input features across both phases and the evolution of the turbulence viscosity within the prediction phase.

Figure 1. Comparisons of $\left( {{q_1},{q_2}} \right)$ in the transitional channel flow in transition region (TR), non-equilibrium region (NE) and equilibrium region (EQ) respectively [Reference Fang, Bao, Xu, Zhou, Du and Jin35].

In the ML-RANS framework, the target variable refers to the object of turbulence modeling, which involves leveraging machine learning techniques to repesent Reynolds stresses, either as mathematical expressions or through artificial neural networks, derived from high-fidelity flow data. However, certain researchers [Reference Thompson, Sampaio, de Bragança Alves, Thais and Mompean54, Reference Poroseva, Colmenares and Murman55] have observed discrepancies between the resolved mean velocity when employing the Reynolds stresses from high-fidelity data to close RANS equations and the mean velocity in the original high-fidelity data. Wu et al. [Reference Wu, Xiao, Sun and Wang28] highlighted that the RANS equations coupled with Reynolds stress closure models could be ill-conditioned, where minor a-priori errors in modeled terms, such as the Reynolds stresses, can escalate into significant a-posteriori errors during the solution of the RANS equations. To tackle this issue, Wu et al. [Reference Wu, Xiao, Sun and Wang28] partitioned the Reynolds stresses into an implicitly treatment (or linear) component involving ${\nu _t}$ and a nonlinear part, i.e.,

(1) \begin{align}\boldsymbol{{a}} = \boldsymbol\tau - \frac{1}{3}tr( \boldsymbol{\tau})\boldsymbol{{I}} = 2{\nu _t}\boldsymbol{{S}} + {\boldsymbol{{a}}_ \bot }.\end{align}

The linear portion corresponds to one aspect of Pope’s proposed mathematical representation of the Reynolds stress [Reference Pope12], which adheres to the Boussinesq eddy viscosity hypothesis. From the high-fidelity dataset, the Reynolds stress tensor τ _h and the mean shear stress tensor S _h are extracted to calculate ${\nu _t}$ by minimising the nonlinear part, resulting in

(2) \begin{align}{\nu _t} = \left| {\frac{{{\tau _{{h_{ij}}}}{S_{{h_{ij}}}}}}{{2{S_{{h_{ij}}}}{S_{{h_{ij}}}}}}} \right|.\end{align}

Their analysis showed that enhancing the implicit treatment of the Reynolds stress, for instance, by adding an increment to the optimal ${\nu _t}$ while correspondingly subtracting from the nonlinear part, can further improve the conditioning of the model, thereby reducing the sensitivity of the RANS equation solution to the unclosed terms [Reference Wu, Xiao, Sun and Wang28]. In earlier work, Wu et al. [Reference Wu, Xiao and Paterson19] had independently trained machine learning models for the linear and nonlinear components of the Reynolds stress, culminating in a well-conditioned and accurate Reynolds stress model. Liu et al. [Reference Liu, Fang, Rolfo, Moulinec and Emerson29, Reference Liu, Song and Fang51, Reference Liu, Fang, Rolfo, Moulinec and Emerson52], building upon this, took a step further by focusing solely on the turbulence viscosity as the target variable, bypassing the nonlinear part of the Reynolds stress altogether. To enhance the effectiveness of the machine learning regression, they first nondimensionalised the turbulence viscosity using the traditional definition $k/\omega $ from the $k - \omega $ model [Reference Wilcox3, Reference Wilcox8], followed by normalisation to arrive at the actual form of the target variable presented in Table 1.

Table 1. The input features and the target variable for the practical ML-RANS model [Reference Liu, Song and Fang51, Reference Liu, Fang, Rolfo, Moulinec and Emerson52]

Figure 2. Mean velocity profiles by use of different high-fidelity methods.

As for the input features, Liu et al. [Reference Liu, Fang, Rolfo, Moulinec and Emerson29] initially derived from Pope’s analysis [Reference Pope12] an essential input feature (i.e., $k/\nu \omega $ ) related to the target variable. Based on empirical evidence, Liu et al. bypassed tensor invariances in the input features to stabilise simulations and introduced five supplementary input features that included turbulence intensity, wall distance, cross-diffusion of $k$ and $\omega $ , along with two SST model variables specifically capturing properties of the viscous sublayer and the turbulent region [Reference Liu, Fang, Rolfo, Moulinec and Emerson29]. In subsequent research [Reference Liu, Song and Fang51, Reference Liu, Fang, Rolfo, Moulinec and Emerson52], they streamlined their approach by retaining only the first two input features and subjected them to normalisation based on their initial design, as detailed in Table 1. Notably, our previous study [Reference Fang, Bao, Xu, Zhou, Du and Jin35] has demonstrated the ability to identify distinct regimes of turbulence development exclusively using these two dimensionless and normalised input features $\left( {{q_1},{q_2}} \right)$ , as illustrated in Fig. 1. In previous practices, DNS datasets [Reference Abe, Kawamura and Matsuo33, Reference Moser, Kim and Mansour34] and RANS calculations [Reference Liu, Fang, Rolfo, Moulinec and Emerson29, Reference Liu, Song and Fang51, Reference Liu, Fang, Rolfo, Moulinec and Emerson52] were used, respectively.

2.2 Non-equilibrium flows

The property of multi-scale energy transfer is a good starting point for understanding turbulence. The well-known Richardson-Kolmogorov energy cascade indicates an equilibrium dissipation law, that can be depicted as a balance of high Reynolds turbulence between the forward energy transfer and the dissipation.

(3) \begin{align}\epsilon = \frac{{{C_\epsilon }{{\mathcal U}^3}}}{{\mathcal L}}.\end{align}

where ${{\mathcal U}^2}$ is turbulent kinetic energy, ${\mathcal L}$ is the integral length scale and the dissipation coefficient ${C_\epsilon }$ is constant. In fact, in real flows there exist various dissipation laws, including the equilibrium or quasi-equilibrium laws and the non-equilibrium laws, as Vassilicos [Reference Vassilicos56] introduced. It was discovered that there exists a non-equilibrium dissipation law between ${C_\epsilon }$ and Taylor-scale Reynolds number ${\rm{Re}}_{\lambda}$ , which is approximately ${C_\epsilon }\sim {\rm{Re}}_\lambda ^{ - 1}$ , or rigorously by theoretical analysis [Reference Bos and Rubinstein57],

(4) \begin{align}{C_\epsilon }\sim {\rm{Re}}_\lambda ^{ - \frac{{15}}{{14}}}.\end{align}

It denotes the imbalance between energy production and dissipation, which mainly exists with a spatially inhomogeneous geometry, even when the flow is a fully developed turbulence. The non-equilibrium canonical flows discussed in sec:intro largely conform to this dissipation rate. Researchers have actively explored the integration of non-equilibrium phenomena within turbulence models by leveraging the ratio between the production and dissipation of turbulence kinetic energy [Reference Wu and Zhang42, Reference Li, Zhang and Chen58]. Recently, it was found that for flows whose transfer spectrum is out-of-equilibrium, there also exists another non-equilibrium dissipation law

(5) \begin{align}{C_\epsilon }\sim {\rm{Re}}_\lambda ^{ - 2},\end{align}

which is usually rather evident in short-time evolution by comparing to the $ - 15/14$ scaling [Reference Liu, Lu, Bos and Fang59–Reference Araki and Bos62] and is closely related to the self-organisation procedure in a flow that is (temporally or spacially) not fully developed. During the evolution of turbulent flows, this novel scaling law typically implies the presence of turbulence energy backscatter (i.e., the transfer of energy from smaller to larger scales). Despite their appearance prior to the establishment of such a scaling law, the conventional turbulence modeling efforts [Reference Liu, Lu, Fang and Gao43, Reference Liu, Tang, Scillitoe and Tucker44] described in sec:intro, aiming at enhancing the Spalart-Allmaras and SST models by incorporating helicity, can be regarded as early attempts to account for this out-of-equilibrium in third-order statistical quantities. However, until now, there has not yet been a data-driven turbulence modeling that specifically addresses the $ - 2$ scaling behaviour associated with energy backscatter phenomena.

In previous studies we showed that a spatially transitional channel flow can perfectly involve all the known states of energy transfer of turbulence evolution, including both the $ - 15/14$ and $ - 2$ scalings of non-equilibrium flows [Reference Liu, Lu, Bos and Fang59]. On another hand, we also illustrated that the flows in turbomachinery are almost non-equilibrium everywhere [Reference Fang, Zhao, Lu, Liu and Yan7]. In the present paper we then use this dataset, aiming at involving various flows in the training procedure, and hoping to improve the performance of turbulence models in turbomachinery.

3.1 Training set

A training set based on the spatially transitional channel flow [Reference Liu, Lu and Fang48] will be introduced and evaluated in this subsection. The channel has a $64{\pi}h \times 2h \times 2{\pi}h$ domain with $h$ channel half-width with a spatial resolution $1,024 \times 96 \times 64$ in $x$ , $y$ , $z$ directions, including a PSB (periodic suction and blowing disturbance) region in the streamwise interval ${L_x} \in \left[ {{\pi}h,2{\pi}h} \right]$ and a fringe region [Reference Ma, Van Doorne, Zhang and Nieuwstadt63] in the streamwise interval ${L_x} \in \left[ {56{\pi}h,64{\pi}h} \right]$ . The existence of two regions contributes to the generation of the spatial transition. The commonly used dynamic Smagorinsky model was selected as the subgrid-scale model for LES of the channel flow. No-slip boundary conditions are imposed at wall (excluding the PSB region), whereas periodic boundary conditions are utilised along the spanwise direction with the LES. Laminar flow distribution is set at inlet and zero-gradient is set at outlet. Furthermore, the grid spacings and Reynolds numbers are documented in Table 2 for more details. More details refer to Refs [Reference Liu, Lu and Fang48, Reference Xu, Zhang, den Toonder and Nieuwstadt64, Reference Liu65].

Table 2. Information of LES on the spatially transitional channel flow [Reference Liu, Lu and Fang48]

The mean flow profile deviates from its parabolic shape within the range ${L_x} \gt 12\pi $ with the maximum value of the mean velocity reaching approximately 1.14 at ${L_x} \approx 23.4\pi $ . This confirms the typical morphology of a turbulent velocity profile beyond ${L_x} \gt 23.4\pi $ , which is in good agreement with DNS results from Moser et al. [Reference Moser, Kim and Mansour34] and those from Abe et al. [Reference Abe, Kawamura and Matsuo33] at the same friction Reynolds number, ${\rm{R}}{{\rm{e}}_\tau }$ . However, the flow had not yet been fully developed until nearly ${L_x} = 31\pi $ , at which the profile of the total shear stress was a horizontal straight line. In the interval, ${L_x} \in \left[ {23.4{\pi}h,31{\pi}h} \right]$ , the dissipation coefficient ${C_\epsilon }$ increased weakly to constant, and the skewness of longitudinal velocity derivative ${S_k}$ decreased slightly to constant. In addition, ${C_\epsilon }$ was inversely proportional to the integral scale Reynolds number at a low Reynolds number. The phenomena were also observed in the non-equilibrium study of grid turbulence [Reference Hearst and Lavoie66]. Consequently, the segment $\left[ {23.4{\pi}h,31{\pi}h} \right]$ along the streamwise direction was classified as NE, while ${L_x} \lt 23.4\pi $ was characterised as TR and ${L_x} \gt 31\pi $ was considered as EQ. The near-wall mean velocity profile at EQ phase is compared with DNS results of Abe et al. [Reference Abe, Kawamura and Matsuo33], as well as the experiment of Hussain and Reynolds [Reference Hussain and Reynolds67]. Pope reviewed the classical wall law, including the linear relation ${U^ + } = {d^ + }$ in the viscous sublayer ( ${d^ + } \lt 5$ ) and the logarithmic law

(6) \begin{align}{U^ + } = \frac{1}{\kappa }{\rm{ln}}({d^+}) + C,\end{align}

due to von Karman for ${d^ + } \gt 30$ , where $\kappa \approx 0.41$ represents the Karman constant and $C$ was taken as 5.2 despite dependance on Reynolds number [Reference Abe, Kawamura and Matsuo33, Reference Pope68]. As shown in Fig. 2, all three simulations and the experiment agree with the linear relation in the viscous sublayer. In the logarithmic region, the DNS result at ${\rm{R}}{{\rm{e}}_\tau } = 640$ closely matches both the theoretical formula and corresponding experimental findings at the same Reynolds number, thereby substantiating the accuracy of the DNS programme employed by Abe et al. At ${\rm{R}}{{\rm{e}}_\tau } = 395$ , the LES result aligns closely with the DNS result, which indicates that the LES result remains highly credible. For in-depth information, see the paper published by Liu et al. [Reference Liu, Lu and Fang48].

In TR, positive ${S_k}$ denotes the energy backscatter signal [Reference Zhao, Liu, Shao, Fang and Dong69], which also occurs after altering the sign of velocities in time-reversed turbulence [Reference Liu, Lu, Bos and Fang59] where the large-scale dynamics cannot feel the effect of energy backscatter as concurrently as small-scale dynamics does. This non-equilibrium behaviour, characterised by a $ - 2$ scaling law, continues to impact the channel flow until the large-scale structures within NE are notably affected. Subsequently, the flows begins to conform to $ - 15/14$ scaling law. The flow eventually transforms into equilibrium turbulence in EQ after the rapid and slow non-equilibrium procedures. The physical information of transition, non-equilibrium and equilibrium turbulence is retrieved from the flow data in the two-dimensional mid-span section of the spatially transitional channel flow in the current study.

3.2 Modeling process

During the training phase, Google’s TensorFlow [Reference Abadi, Barham, Chen, Chen, Davis, Dean, Devin, Ghemawat, Irving, Isard, Kudlur, Levenberg, Monga, Moore, Murray, Steiner, Tucker, Vasudevan, Warden, Wicke, Yu and Zheng53] was adopted as the platform for machine learning. The training dataset was separated randomly into the training part and the validation one in a ratio of four to one. According to training effect, the hyperparameters in Table 3 were adopted to construct a feedforward artificial neural network (i.e., multi-layer perception). For sake of prevention from overfitting [Reference Tibshirani70], L1 regularisation was adopted to minimise the weight of each layer in the neural network. The total loss was thus defined as

(7) \begin{align}{\rm{Loss}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N {\left(\widehat{\nu _t^{0,i}} - \nu _t^{0,i}\right)^2} + \lambda \mathop \sum \limits_{i,j}^{m,n} | {{w_{i,j}}} |.\end{align}

where $\nu _t^{0,i}$ is the target variable in the training data, $\widehat{\nu _t^{0,i}}$ is the value predicted by the neural network, and ${w_{i,j}}$ is the weight in the neural network. $\lambda $ represents the L1 regularisation coefficient. $N$ is the number of training data, $m$ and $n$ are respectively the numbers of layers and nodes in the neural network. The first part of the loss is the mean square error (MSE) from the training data, and the second part is the L1 regularisation penalty. As seen in Fig. 3, the evolution of errors indicates that our model is rapidly established once the training process starts. To estimate visually the model result, we plot the profiles of the target variable $\nu _t^0$ at different streamwise locations in Fig. 4, and the relation between the target variable and both input features is shown in Fig. 5. It can be observed that the training data and ML predictions are all in good agreement.

Table 3. Hyperparameters of the neural network for the ML-RANS TR-NE-EQ model

Figure 3. Evolution of loss in the training process. The x-axis represents the ratio of current training step and total steps.

Figure 4. Profiles of the target variable $\nu _t^0$ at $x/( {{\pi}h} ) = 20,30,40,50$ in TR (blue), NE (red) and EQ (green) regions: comparison of ML predictions (lines) and a snapshot of training dataset (points).

Figure 5. Relation between the target variable $\nu _t^0$ and two input features $\left( {{q_1},{q_2}} \right)$ at four different streamwise locations: comparison of ML predictions (line) and a snapshot of training dataset (point).

Following training, a comprehensive data-driven turbulence model is created by combining the machine learning model with the predictor library [Reference Liu, Song and Fang51]. Figure 6 illustrates that the target variable $\nu _t^0$ is pretty independent of the input features where ${q_1}$ is large for the ML-RANS EQ model, while it decreases obviously as ${q_2}$ increases in the same region for the ML-RANS TR-NE-EQ model.

Figure 6. Mapping from $\left( {{q_1},{q_2}} \right)$ to $\nu _t^0$ in two ML-RANS models.

3.3 A-posteriori test

Having implemented the novel ML-RANS TR-NE-EQ model in OpenFOAM, we applied it into the spatially transitional channel flow used during training, comparing its performance against the earlier ML-RANS EQ model and the baseline model. Inspired by Zhao’s doctoral dissertation [Reference Zhao71] where several turbulence models were utilised for RANS simulations of this transitional channel flow, this study followed suit by removing the fringe area close to the outlet. A computationally economical mesh configuration of $440 \times 50 \times 32$ was established across the computational domain of $55{\pi}h \times 2h \times 2{\pi}h$ , ensuring effective refinement near the top and bottom walls, particularly with the first cell layer meeting the requirement of ${d^ + } \lt 1$ .

Consistent with LES conducted by Liu et al. [Reference Liu, Lu and Fang48], the transition is facilitated by the PSB region near the inlet; thus, for the wall segment where the streamwise coordinates ${L_x}$ lie within $\left[ {{\pi}h,2{\pi}h} \right]$ , we also designed the PSB boundary condition for the velocity field while imposing a zero-gradient condition for other physical quantities. For the remaining walls, we set the no-slip condition with low-Reynolds number corrections. In the meantime, laminar flow inlet and zero-gradient outlet were configured accordingly. Considering that the velocity boundary values vary over time in the PSB region, we conducted unsteady RANS simulations using the baseline turbulence model, ML-RANS EQ model, and the ML-RANS TR-NE-EQ model, each assisted by the PISO algorithm. To maintain numerical stability, the time step was constrained to ensure a Courant number below 1.

For incompressible flows, the bulk mean velocity remains constant, and the transition from laminar to turbulence can be reflected through variations in the maximum value of the mean velocity. Upon comparing the simulation outputs generated by the three turbulence models against the maximum value of the mean velocity from the LES data, it is observed that the ML-RANS TR-NE-EQ model delivers superior performance among three RANS turbulence models, as shown in Fig. 7. Firstly, both ML-RANS models effectively avoid the occurrence of the local velocity peak within the PSB region simulated by the baseline model. Subsequently, following the PSB region, the maximum value of the mean velocity computed by all three models begin to decline, which indicates that the flow starts to deviate from the laminar. However, the ML-RANS TR-NE-EQ model exhibits a notably delayed initiation of flow transition, as its calculated mean velocity maximum doesn’t experience a sharp drop until reaching ${L_x} = 10{\pi}h$ , aligning closely with the LES-derived transitional point ${L_x} = 12{\pi}h$ . The RANS simulation adopting the ML-RANS TR-NE-EQ model presents remarkable agreement with the LES results in the NE region. Within the EQ region where the large-scale non-equilibrium is virtually absent, the maximum value of the mean velocity calculated by all four numerical simulations converge. Among three RANS models, the ML-RANS TR-NE-EQ model’s prediction shows a closer resemblance to the LES result than that of the baseline model. In summary, the ML-RANS TR-NE-EQ model provides the best simulation of the transitional channel flow among the three RANS turbulence models.

Figure 7. Evolution of maximum value of the mean velocity along streamwise direction.