2.1. Turbulence modelling from a kinetic perspective

The study and simulation of turbulence can be categorised into three distinct levels: the microscopic level, the hydrodynamic level, and the average flow field level. Each of these levels corresponds to a specific characteristic length (Wilcox Reference Wilcox1998; Yakhot Reference Yakhot2007): the molecular mean free path scale $\lambda$ , the Kolmogorov length scale ${\eta _K} \approx L{{{Re}} ^{ - 3/4}}$ , Re is the Reynolds number and the macroscopic length scale $L$ , as depicted in figure 1.

Figure 1.

Schematic diagram of scales of length and velocity fields defined in turbulence.

At the microscopic level ( $r \approx \lambda$ ), our focus begins with molecules, where the velocity fields emphasise molecular velocities – specifically, the molecules’ instantaneous velocities $\boldsymbol{c}$ , and their thermal velocities $\boldsymbol{c}'$ . The thermal velocity can be determined by $\boldsymbol{c}' = \boldsymbol{c} - \boldsymbol{v}$ , where $\boldsymbol{v}$ represents the instantaneous velocity of a fluid particle. At the hydrodynamic level, fluid particles serve as the fundamental elements. Hydrodynamic approximations can be obtained by eliminating small-scale modes from the interval $r \leqslant {\eta _K}$ , resulting in the velocities of fluid particles – specifically, their instantaneous velocities $\boldsymbol{v}$ , and their fluctuating velocities $\boldsymbol{v}'$ . The fluctuating velocity can be determined by $\boldsymbol{v}' = \boldsymbol{v} - \bar{\boldsymbol{v}}$ , where $\bar{\boldsymbol{v}}$ represents the mean velocity at the average flow field level. Subsequently, turbulence models can be derived by averaging over the interval ${\eta _K} \leqslant r \leqslant L$ , focusing solely on mean velocities.

In principle, kinetic equations, such as the Boltzmann equation and its simplified models, can provide a comprehensive description of turbulence at the microscopic level (Gallis et al. Reference Gallis, Bitter, Koehler, Torczynski, Plimpton and Papadakis2017; McMullen et al. Reference McMullen, Krygier, Torczynski and Gallis2022). Based on this foundation, there are two fundamental approaches for deriving turbulence models at the average flow field level from kinetic equations (Succi Reference Succi2020). As depicted in figure 2, the first approach involves initially projecting out the velocity degrees of freedom from the kinetic equations to derive conservation equations – namely, the mass and momentum conservation equations – at the hydrodynamic level. The appearance of molecular viscous stress makes these equations unclosed. Assuming a weak departure from local equilibrium, we can apply the Chapman–Enskog expansion (Chapman & Cowling Reference Chapman and Cowling1990) to link the molecular viscous stress with hydrodynamic variables. Truncating this expansion to the first order yields the NS equations. By further coarse-graining the NS equations through ensemble averaging, we can derive the RANS equations with Reynolds stress, which forms the foundation of mainstream turbulence modelling.

The second approach begins by constructing a kinetic equation for turbulence at the hydrodynamic level from the molecular kinetic equations, either through coarse-graining or by analogy. Then by projecting out the velocity degrees of freedom from this kinetic equation, the RANS equations can be derived. To close these equations, the Chapman–Enskog expansion can be used to derive a constitutive relationship between the Reynolds stress and the average flow field variables. This methodology is adopted in the current study.

Figure 2.

Two basic approaches to coarse-graining the kinetic equations to obtain turbulence models at the average flow field level.

Symbolically, the first approach (the hydrodynamic approach) can be expressed as

(2.1) \begin{equation} {U_h}\left ( {x;t} \right ) = {P_x} \, {P_v} \, f\left ( {\boldsymbol{x},\boldsymbol{v};t} \right ), \end{equation}

and the second approach (the kinetic approach) can be expressed as

(2.2) \begin{equation} {U_k}\left ( {x;t} \right ) = {P_v} \, {P_x} \, f\left ( {\boldsymbol{x},\boldsymbol{v};t} \right ). \end{equation}

In the above, $P_x$ denotes a space-filtering projection operator, while $P_v$ denotes the velocity projector associated with Chapman–Enskog expansion (Succi Reference Succi2020). For turbulence modelling, these two projectors do not commute, which indicates the potential of the kinetic approach to yield new turbulence models.

2.2. The coarse-grained Fokker–Planck equation for turbulence

2.2.1. The microscopic Fokker–Planck equation

At the microscopic level, the positions $\boldsymbol{x^*}{ ( t )}$ and velocities $\boldsymbol{c^*}{ ( t )}$ of molecules can be described using the Boltzmann equation. Alternatively, they can also be considered as a stochastic diffusion process (Heinz Reference Heinz2003). This process can be represented by deriving an equation for the motion of a single molecule colliding with all other molecules, resulting in the Langevin equation (Langevin Reference Langevin1908). Assuming the neglect of possible contributions from external forces, the Langevin equation can be expressed as

(2.3) \begin{equation} \frac {{\mathrm{d}x_i^*}}{{\mathrm{d}t}} = c_i^*,\quad \frac {{\mathrm{d}c_i^*}}{{\mathrm{d}t}} = - \frac {1}{{{\tau _m}}}\left ( {c_i^* - {u_i}} \right ) + \sqrt {\frac {{4e}}{{3{\tau _m}}}} \frac {{\mathrm{d}{W_i}}}{{\mathrm{d}t}}. \end{equation}

Notably, the mean velocity $u_i$ is equivalent to the fluid velocity measured on the hydrodynamic scale, in accordance with the ergodic theorem. According to (2.3), the molecular velocity is driven by a drift towards the hydrodynamic velocity and a stochastic noise term governed by the derivative of a Wiener process $W ( t )$ . A Wiener process is a Gaussian process with $\mathrm{d}{W_i} ( t ) = {W_i} ( {t + {\rm d}t} ) - {W_i} ( t )$ , $\left \langle {\mathrm{d}{W_i} ( t )} \right \rangle \equiv 0$ and $\left \langle {\mathrm{d}{W_i} ( t )\,\mathrm{d}{W_j} ( t )} \right \rangle \equiv \mathrm{d}t\,{\delta _{ij}}$ , where the angle brackets denote the expectation of the stochastic processes inside. Here, $\tau _m$ is a molecular relaxation time characterising the molecular velocity relaxation and is related to the molecular viscosity, while $e$ is the sensible energy relevant to the intensity of molecular fluctuations. All molecular statistics can be obtained by directly solving (2.3). For an Eulerian description, (2.3) can be reformulated as an equivalent Fokker–Planck equation for the molecular velocity distribution function $f ( {{c_i},{x_i},t} )$ , i.e.

(2.4) \begin{equation} \frac {{\partial f}}{{\partial t}} + {c_i}\frac {{\partial f}}{{\partial {x_i}}} = \frac {1}{{{\tau _m}}}\frac {\partial }{{\partial {c_k}}} [ {f( {{c_k} - {u_k}})} ] + \frac {{2e}}{{3{\tau _m}}}\frac {{{\partial ^2}f}}{{\partial {c_k}\,\partial {c_k}}}. \end{equation}

By multiplying (2.4) by the conservation variables and integrating the resulting equations over the velocity space, one obtains the conservation equations at the hydrodynamic level (Heinz Reference Heinz2004; Jenny et al. Reference Jenny, Torrilhon and Heinz2010; Singh & Ansumali Reference Singh and Ansumali2015).

2.2.2. The Fokker–Planck equation for turbulence

Equation (2.3) or (2.4) holds at the finest level of resolution, namely from the macroscopic length scale $L$ all the way down to the molecular mean free path scale $\lambda$ . However, such a microscale description is computationally prohibitive for most turbulent flows. Thus a kinetic model for turbulence at the hydrodynamic level is highly required – one that omits the intricate details of molecular motions.

Although differences between kinetic theory and statistical turbulence theory exist, as Taylor (Reference Taylor1935) pointed out, one can still leverage the fundamental similarities underlying both statistical phenomena. In Brownian stochastic processes, molecular dynamics is modelled using two characteristic times: a long characteristic time for the non-equilibrium degree of freedom, embedded in a drift term, and a short characteristic time for other degrees of freedom that are in equilibrium, embedded in the Wiener process of the diffusion term. In turbulent flows, eddies of various sizes intermingle, and the statistical properties of one fluid particle result from the combined influence of these different-sized eddies. This reveals certain similarities between turbulence and Brownian motion. In addition, it has been established that two statistically independent characteristic times exist in turbulence – a long one for lower wavenumbers, and a short one for higher wavenumbers. According to Kolmogorov’s hypotheses (Kolmogorov Reference Kolmogorov1962), these two characteristic times are sufficiently separated when the Reynolds number is large, enabling the application of Brownian statistics to describe turbulent flow.

Recently, the experiments by Chong et al. (Reference Chong, Shi, Ding, Ding, Lu, Zhong and Xia2020) further demonstrate that vortices can be treated as Brownian particles in a rotating turbulent convective flow. Based on these observations, we propose the following Langevin equation for turbulent flow at the hydrodynamic level:

(2.5) \begin{align} \frac {{\mathrm{d}x_i^*}}{{\mathrm{d}t}} = u_i^*, \quad \frac {{\mathrm{d}u_i^*}}{{\mathrm{d}t}} = \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} - \frac {1}{{{\tau _T}}}( {u_i^* - {{\bar u}_i}} ) + \sqrt {\frac {{4k}}{{3{\tau _T}}}} \frac {{\mathrm{d}{W_i}}}{{\mathrm{d}t}}. \end{align}

And the corresponding Fokker–Planck equation can be obtained as

(2.6) \begin{equation} \frac {{\partial f}}{{\partial t}} + {u_i}\frac {{\partial f}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial f}}{{\partial {u_i}}} = \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_k}}} [ {f( {{u_k} - {{\bar u}_k}} )} ] + \frac {{2k}}{{3{\tau _T}}}\frac {{{\partial ^2}f}}{{\partial {u_k}\,\partial {u_k}}}. \end{equation}

The analogy between molecular and turbulent motions provides the basis for this model. Unlike (2.3) and (2.4), the velocities in (2.5) and (2.6) are at the hydrodynamic scale. Here, $\bar u_i$ is the mean velocity at the average flow field scale according to the ergodic theorem. The distribution function $f ( {{u_i},{x_i},t} )$ describes the occupancy probability of the fluid particles in the phase space. Fluctuations are characterised by energy and a characteristic relaxation time. It should be noted that for modelling turbulent motions, the molecular energy $e$ and relaxation time scale $\tau _m$ in (2.3) and (2.4) are replaced by the turbulent kinetic energy $k$ and the relaxation time $\tau _T$ for turbulence in (2.5) and (2.6). Notably, a term relevant to molecular viscous stress, defined as ${\sigma _{ij}} = - \overline {\rho {c'_i}{c'_j}}$ , is introduced in (2.5) and (2.6) to ensure the correct transport of mean velocities. As demonstrated in Appendix A, by taking moments of (2.6), one obtains the RANS equations and a partially modelled Reynolds stress equation at the average flow field level.

By focusing on the probabilistic evolution of velocity fluctuations, the Fokker–Planck model, as shown in (2.6), effectively captures the statistical nature of turbulence. The drift and diffusion terms account for the integrated effects of eddy interactions. Specifically, the drift term reflects the mean flow behaviour or large-scale structures, corresponding to coherent eddies, and is modelled as a relaxation towards the mean velocity, with $\tau _T$ characterising the relaxation time scale. As will be demonstrated in §§ 2.3 and 2.4, $\tau _T$ can be related to the eddy viscosity with an appropriate choice. The diffusion term captures the random, small-scale velocity fluctuations that arise from the chaotic breakdown of larger structures. While coherent structures are dependent on flow conditions, Kolmogorov’s hypotheses suggest that small-scale chaotic motions exhibit a surprising degree of independence from the larger scales. Thus these small-scale fluctuations are modelled using a Wiener process, with turbulent energy $k$ controlling their intensity.

2.2.3. Connection with the PDF method

Pope (Reference Pope1985, Reference Pope2001) derives the exact transport equation for $f ( {{u_i},{x_i},t} )$ as a hydrodynamic description from the NS equations. He then proposes a generalised Langevin model (GLM) to close this transport equation, which is the well-known PDF method. The form of GLM can be expressed as

(2.7) \begin{equation} \frac {{\mathrm{d}x_i^*}}{{\mathrm{d}t}} = u_i^*, \quad \frac {{\mathrm{d}u_i^*}}{{\mathrm{d}t}} = \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + {G_{ij}}( {u_j^* - {{\bar u}_j}} ) + \sqrt {{C_0}\varepsilon } \frac {{\mathrm{d}{W_i}}}{{\mathrm{d}t}}. \end{equation}

And the corresponding PDF transport equation is

(2.8) \begin{equation} \frac {{\partial f}}{{\partial t}} + {u_i}\frac {{\partial f}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial f}}{{\partial {u_i}}} = - \frac {\partial }{{\partial {u_i}}}\left [ {f{G_{ij}}\left ( {{u_j} - {{\bar u}_j}} \right )} \right ] + \frac {1}{2}{C_0}\varepsilon \frac {{{\partial ^2}f}}{{\partial {u_k}\,\partial {u_k}}}. \end{equation}

In fact, the GLM is a family of models, where ${G_{ij}} ( {\boldsymbol{x},t} )$ and ${C_0} ( {\boldsymbol{x},t} )$ are coefficients that define each specific model. By calibrating these coefficients against a broad range of experiments, Haworth & Pope (Reference Haworth and Pope1986, Reference Haworth and Pope1987) develop a model that exhibits satisfactory performance in both homogeneous shear flow and free shear flows. Furthermore, it can be contracted to the isotropisation of production model (IPM). According to Fu & Pope (Reference Fu and Pope1994) and Pope (Reference Pope2001), the GLM can be considered a nonlinear extension of Reynolds stress transport models, such as Rotta’s model (Rotta Reference Rotta1951) and the Launder–Reece–Rodi model (Launder et al. Reference Launder, Reece and Rodi1975). Fu & Pope (Reference Fu and Pope1994) derive the Reynolds stress transport model corresponding to the Haworth–Pope model, and assess its performance in complex flow conditions induced by strong swirl. The results are compared with those from IPM and the $k\text{-}\varepsilon$ model, revealing that the IPM demonstrates the best agreement with experimental data regarding the size of the recirculation, while the Haworth–Pope model tends to produce an excessively strong recirculation. However, when it comes to the profiles of mean flow and turbulence quantities, none of the models outperforms the others significantly.

Notably, our proposed model exhibits a form similar to the GLM. If we specify $G_{ij}$ as isotropic and the coefficient $C_0$ as

(2.9) \begin{align} {G_{ij}} &= - {\delta _{ij}}/{\tau _T},\nonumber \\[3pt] {C_0} &= \frac {4}{{3{\tau _T}}}\frac {k}{\varepsilon },\end{align}

the GLM ((2.7) and (2.8)) recovers (2.5) and (2.6). This implies that our proposed Fokker–Planck equation for turbulence can also be regarded as a specific instance of the GLM. In the next two subsections, the connection between the turbulent Fokker–Planck equation and the Reynolds stress constitutive relationship will be established.

2.3. First-order expansion: linear EVMs

Building on the derivation by Chen et al. (Reference Chen, Orszag, Staroselsky and Succi2004) for the BGK model, we employ a Chapman–Enskog expansion technique (Chapman & Cowling Reference Chapman and Cowling1990) to derive turbulence models from the Fokker–Planck equation for turbulence, as illustrated in (2.6). We first use a reference velocity $u_r$ , a reference collision time $\tau _r$ and a reference length $L$ to non-dimensionalise (2.6) as

(2.10) \begin{equation} \frac {{\partial \hat f}}{{\partial \hat t}} + {\hat u_j}\frac {{\partial \hat f}}{{\partial {{\hat x}_j}}} + \frac {1}{{\hat \rho }}\frac {{\partial {{\hat \sigma }_{ij}}}}{{\partial {{\hat x}_i}}}\frac {{\partial \hat f}}{{\partial {{\hat u}_j}}} = \frac {1}{{{{\hat \tau }_T}}}\frac {\partial }{{\partial {{\hat u}_i}}}\left ( {{{\hat u}_i}'\hat f + \frac {2}{3}\hat k\frac {{\partial \hat f}}{{\partial {{\hat u}_i}}}} \right ) \frac {1}{\xi }. \end{equation}

In (2.10), $u'$ is the fluctuation velocity, and $\xi$ is defined as

(2.11) \begin{equation} \xi = \frac {{{u_r}{\tau _r}}}{L}. \end{equation}

In essence, $\xi$ represents the ratio of the unresolved scale of motion to the resolved one, and can therefore be interpreted as a Knudsen number based on turbulent fluctuations (Righi Reference Righi2016). Notably, it vanishes when returning to the dimensional form. For simplicity, we will omit the caret placed over the symbol in the following context. If we assume that $\xi$ is a small quantity less than unity, then we can expand $f$ as the power series

(2.12) \begin{equation} f = {f^{\left ( 0 \right )}} + \xi {f^{\left ( 1 \right )}} + {\xi ^2}{f^{\left ( 2 \right )}} + \cdots . \end{equation}

Substituting (2.12) into (2.10), we obtain

(2.13) \begin{equation} \begin{aligned} &- \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_i}}}\left [ {{u'_i}{f^{\left ( 0 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}}} \right ] \\[6pt] &{}+ \xi \left [ {\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}} - \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 1 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}}} \right )} \right ]\\[6pt] &{}+ {\xi ^2}\left [ {\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}} - \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 2 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 2 \right )}}}}{{\partial {u_i}}}} \right )} \right ] = o( {{\xi ^3}} ). \end{aligned} \end{equation}

Note that $f^{ ( 0 )}$ is an equilibrium distribution that does not include the non-trivial flow-induced fluctuations. Following the previous work of Chen et al. (Reference Chen, Kandasamy, Orszag, Shock, Succi and Yakhot2003, Reference Chen, Orszag, Staroselsky and Succi2004), it is reasonable to use a local Gaussian form, specifically the Maxwell distribution

(2.14) \begin{equation} {f^{\left ( 0 \right )}} = \frac {n}{{{{\left ( {4/3\pi k} \right )}^{3/2}}}} \exp \left ( { - \frac {{\boldsymbol{u}{'^2}}}{{4k/3}}} \right ), \end{equation}

to represent $f^{ ( 0 )}$ . Notably, measurements from experiments and DNS in homogeneous turbulence (Tavoularis & Corrsin Reference Tavoularis and Corrsin1981; Yeung & Pope Reference Yeung and Pope1989) indicate that the one-point velocity distribution function is indeed jointly normal. Consequently, we find that the first term of (2.13) vanishes.

We further expand the time derivation as

(2.15) \begin{equation} \left \{ \begin{aligned} \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial t}} = \frac {{{\partial _0}{f^{\left ( 0 \right )}}}}{{\partial t}} + \xi \frac {{{\partial _1}{f^{\left ( 0 \right )}}}}{{\partial t}} + {\xi ^2}\frac {{{\partial _2}{f^{\left ( 0 \right )}}}}{{\partial t}} + \cdots, \\[6pt] \frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial t}} = \frac {{{\partial _0}{f^{\left ( 1 \right )}}}}{{\partial t}} + \xi \frac {{{\partial _1}{f^{\left ( 1 \right )}}}}{{\partial t}} + {\xi ^2}\frac {{{\partial _2}{f^{\left ( 1 \right )}}}}{{\partial t}} + \cdots. \end{aligned} \right . \end{equation}

Substituting (2.15) into (2.13), after rearrangement we obtain

(2.16) \begin{align} &\left [ {\frac {{{\partial _0}{f^{\left ( 0 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}} - \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 1 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}}} \right )} \right ]\nonumber\\[6pt] & + \xi \left [ {\frac {{{\partial _1}{f^{\left ( 0 \right )}}}}{{\partial t}} + \frac {{{\partial _0}{f^{\left ( 1 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}} - \frac {1}{{{\tau _T}}}\frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 2 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 2 \right )}}}}{{\partial {u_i}}}} \right )} \right ] \nonumber\\[6pt] & = o ( {{\xi ^2}} ). \end{align}

Consequently, (2.10) is transformed into an infinite hierarchy of equations based of the order of $\xi$ . In particular,

(2.17) \begin{equation} {f^{\left ( 0 \right )}} = {f_M} = \frac {n}{{{{\left ( {4/3\pi k} \right )}^{3/2}}}} \exp \left ( { - \frac {{\boldsymbol{u}{'^2}}}{{4k/3}}} \right ), \end{equation}

(2.18) \begin{equation} \frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 1 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}}} \right ) = {\tau _T}\left ( {\frac {{{\partial _0}{f^{\left ( 0 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}}} \right ), \end{equation}

(2.19) \begin{equation} \frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 2 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 2 \right )}}}}{{\partial {u_i}}}} \right ) = {\tau _T}\left ( {\frac {{{\partial _1}{f^{\left ( 0 \right )}}}}{{\partial t}} + \frac {{{\partial _0}{f^{\left ( 1 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}}} \right ). \end{equation}

Next, we will determine the solution corresponding to the first-order expansion. Since the spatial and time dependence of $f^{ ( 0 )}$ is only through $\rho$ , $\bar u_i$ and $k$ , we can organise the right-hand side of (2.18) into

(2.20) \begin{align} & {\tau _T}\left ( \frac {{{\partial _0}{f^{\left ( 0 \right )}}}}{{\partial t}} \right .+\left . {u_i}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}} \right ) \nonumber \\[6pt] & \qquad = {\tau _T} \left [ \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial \rho }} \left ( {\frac {{{\partial _0}\rho }}{{\partial t}} + {u_j}\frac {{\partial \rho }}{{\partial {x_j}}}} \right ) + \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {{\bar u}_i}}}\left ( {\frac {{{\partial _0}{{\bar u}_i}}}{{\partial t}} + {u_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) \right . \nonumber \\[6pt] & \qquad \quad+ \left . \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial k}}\left ( {\frac {{{\partial _0}k}}{{\partial t}} + {u_j}\frac {{\partial k}}{{\partial {x_j}}}} \right ) - \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} \times {f^{\left ( 0 \right )}}\frac {{{u'_i}}}{{2k/3}} \right ]. \end{align}

Using the expression for $f^{ ( 0 )}$ , it can be demonstrated directly that

(2.21) \begin{equation} \left \{ \begin{aligned} \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial \rho }} &= \frac {{{f^{\left ( 0 \right )}}}}{\rho },\\[6pt] \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {{\bar u}_i}}} &= {f^{\left ( 0 \right )}}\frac {{{u'_i}}}{{2k/3}},\\[6pt] \frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial k}} &= {f^{\left ( 0 \right )}}\left [ {\frac {{{u'_i}{u'_i}}}{{4{k^2}/3}} - \frac {3}{{2k}}} \right ]. \end{aligned} \right . \end{equation}

For the time derivation terms involved in (2.20), we will use the RANS equation and the governing equation of $k$ , i.e.

(2.22) \begin{equation} \left \{ \begin{aligned} \frac {{\partial \rho }}{{\partial t}} &= - \rho \frac {{\partial {{\bar u}_j}}}{{\partial {x_j}}} - {{\bar u}_j}\frac {{\partial \rho }}{{\partial {x_j}}},\\[6pt] \frac {{\partial {{\bar u}_i}}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}^T}}}{{\partial {x_j}}},\\[6pt] \frac {{\partial k}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial k}}{{\partial {x_j}}} + \frac {1}{\rho }{\sigma _{ij}^T}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac {\partial }{{\partial {x_j}}}\left ( {\nu _{eff}^k\frac {{\partial k}}{{\partial {x_j}}}} \right ) - \varepsilon . \end{aligned} \right . \end{equation}

In the above equations, $\nu _{eff}^k = \nu + {\nu _T}/{\sigma _k}$ is the effective diffusivity for $k$ , and ${\sigma _{ij}^T} = - \overline {\rho {u'_i}{u'_j}}$ is the Reynolds stress. We can expand $\sigma _{ij}^T$ into a power series as

(2.23) \begin{equation} {\sigma _{ij}^T} = - \frac {2}{3}\rho k \delta _{ij} + \sigma _{ij}^{T\left ( 0 \right )} + \xi \sigma _{ij}^{T\left ( 1 \right )} + {\xi ^2}\sigma _{ij}^{T\left ( 2 \right )} + \cdots, \end{equation}

where $\sigma _{ij}^{T ( \alpha )} = - \int _{{R^3}} {mu{'_{\left \langle i \right .}}u{'_{\left . j \right \rangle }}} {f^{ ( \alpha )}}\,\mathrm{d}\boldsymbol{u}$ . Notably, to enhance clarity and conciseness, we use angle brackets to denote the symmetric and deviatoric (or trace-free) part of the tensor. Its complete form is given by

(2.24) \begin{equation} u{'_{\left \langle i \right .}}u{'_{\left . j \right \rangle }} = \frac {1}{2} ( {u{'_i}u{'_j} + u{'_j}u{'_i}} ) - \frac {1}{3}u{'_k}u{'_k}{\delta _{ij}}. \end{equation}

Additionally, the Einstein summation convention is applied here and throughout the text. It is easy to demonstrate that

(2.25) \begin{equation} \sigma _{ij}^{T\left ( 0 \right )} = - \int _{{R^3}} {mu{'_{\left \langle i \right .}}u{'_{\left . j \right \rangle }}} {f^{\left ( 0 \right )}}\,\mathrm{d}\boldsymbol{u} = 0. \end{equation}

Thus

(2.26) \begin{equation} {\sigma _{ij}^T} = - \frac {2}{3}\rho k \delta _{ij} + \xi \sigma _{ij}^{T\left ( 1 \right )} + {\xi ^2}\sigma _{ij}^{T\left ( 2 \right )} + \cdots . \end{equation}

Substituting (2.26) into (2.22), we obtain

(2.27) \begin{equation} \left \{ \begin{aligned} \frac {{\partial \rho }}{{\partial t}} &= - \rho \frac {{\partial {{\bar u}_j}}}{{\partial {x_j}}} - {{\bar u}_j}\frac {{\partial \rho }}{{\partial {x_j}}},\\[6pt] \frac {{\partial {{\bar u}_i}}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} - \frac {1}{\rho }\,\frac {2}{3}\frac {{\partial \rho k{\delta _{ij}}}}{{\partial {x_j}}} + \frac {1}{\rho }\sum \limits _\alpha {{\xi ^\alpha }\frac {{\partial \sigma _{ij}^{T\left ( \alpha \right )}}}{{\partial {x_j}}}}, \\[6pt] \frac {{\partial k}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial k}}{{\partial {x_j}}} - \frac {1}{\rho } \cdot \frac {2}{3}\rho k\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}{\delta _{ij}} + \frac {1}{\rho }\sum \limits _\alpha {{\xi ^\alpha }\sigma _{ij}^{T\left ( \alpha \right )}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} - \nu _{eff}^k\frac {{{\partial ^2}k}}{{\partial {x_j}^2}} - \varepsilon . \end{aligned} \right . \end{equation}

In this manner, the time derivation of $\rho$ , $\bar u_i$ and $k$ can be expanded in a power series as

(2.28) \begin{equation} \left \{ \begin{aligned} \frac {{{\partial _0}\rho }}{{\partial t}} &= - \rho \frac {{\partial {{\bar u}_j}}}{{\partial {x_j}}} - {{\bar u}_j}\frac {{\partial \rho }}{{\partial {x_j}}},\\[6pt] \frac {{{\partial _0}{{\bar u}_i}}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}} - \frac {1}{\rho }\,\frac {2}{3}\frac {{\partial \rho k}}{{\partial {x_i}}},\\[6pt] \frac {{{\partial _0}k}}{{\partial t}} &= - {{\bar u}_j}\frac {{\partial k}}{{\partial {x_j}}} - \frac {1}{\rho } \, \frac {2}{3}\rho k\frac {{\partial {{\bar u}_i}}}{{\partial {x_i}}} - \nu _{eff}^k\frac {{{\partial ^2}k}}{{\partial {x_j}^2}} - \varepsilon, \end{aligned} \right .\left ( {\alpha = 0} \right ), \end{equation}

(2.29) \begin{equation} \left \{ \begin{aligned} \frac {{{\partial _\alpha }\rho }}{{\partial t}} &= 0,\\[6pt] \frac {{{\partial _\alpha }{{\bar u}_i}}}{{\partial t}} &= \frac {1}{\rho }{\xi ^\alpha }\frac {{\partial \sigma _{ij}^{T\left ( \alpha \right )}}}{{\partial {x_j}}},\\[6pt] \frac {{{\partial _\alpha }k}}{{\partial t}} &= \frac {1}{\rho }{\xi ^\alpha }\sigma _{ij}^{T\left ( \alpha \right )}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}. \end{aligned} \right .\left ( {\alpha \geqslant 1} \right ). \end{equation}

Substituting (2.21) and (2.28) into (2.20), and omitting terms that are associated with the second derivative of $k$ and $\varepsilon$ for simplicity, we then obtain

(2.30) \begin{align} {\tau _T}&\left ( {\frac {{{\partial _0}{f^{\left ( 0 \right )}}}}{{\partial t}} + {u_i}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {x_i}}} + \frac {1}{\rho }\frac {{\partial {\sigma _{ij}}}}{{\partial {x_j}}}\frac {{\partial {f^{\left ( 0 \right )}}}}{{\partial {u_i}}}} \right ) \nonumber \\[6pt] &= \frac {{{\tau _T}}}{{2k/3}}{f^{\left ( 0 \right )}}\left [ {\left ( {{u'_i}{u'_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac {1}{3}{u'_i}{u'_i}\frac {{\partial {{\bar u}_j}}}{{\partial {x_j}}}} \right ) + {u'_j}\left ( {\frac {{\boldsymbol{u}{'^2}}}{2} - \frac {5}{2} \cdot \frac {2}{3}k} \right )\frac {{\partial \ln k}}{{\partial {x_j}}}} \right ]. \end{align}

We can transform the left-hand side of (2.18) into

(2.31) \begin{equation} \frac {\partial }{{\partial {u_i}}}\left ( {{u'_i}{f^{\left ( 1 \right )}} + \frac {2}{3}k\frac {{\partial {f^{\left ( 1 \right )}}}}{{\partial {u_i}}}} \right ) = \frac {\partial }{{\partial {u_i}}}\left [ {\frac {2}{3}k \, {f_M}\frac {\partial }{{\partial {u_i}}}\left ( {\frac {{{f^{\left ( 1 \right )}}}}{{{f_M}}}} \right )} \right ]. \end{equation}

Combining with (2.30), (2.18) can be written as

(2.32) \begin{align} \frac {1}{{{f_M}}}&\frac {\partial }{{\partial {u_i}}}\left [ {\frac {2}{3}k \, {f_M}\frac {\partial }{{\partial {u_i}}}\left ( {\frac {{{f^{\left ( 1 \right )}}}}{{{f_M}}}} \right )} \right ] \nonumber \\[6pt] &= \frac {{{\tau _T}}}{{2k/3}}\left [ {\left ( {{u'_i}{u'_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac {1}{3}{u'_i}{u'_i}\frac {{\partial {{\bar u}_j}}}{{\partial {x_j}}}} \right ) + {u'_j}\left ( {\frac {{\boldsymbol{u}{'^2}}}{2} - \frac {5}{2} \, \frac {2}{3}k} \right )\frac {{\partial \ln k}}{{\partial {x_j}}}} \right ]. \end{align}

The left-hand side of (2.32) is a linear operator of ${f^{ ( 1 )}}/{f_M}$ , so we can denote it as $L [ {{f^{ ( 1 )}}/{f_M}} ]$ . If we define

(2.33) \begin{align} {A_{ij}} &= {u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}, \end{align}

(2.34) \begin{align} {B_j} &= {u'_j}\left ( {\frac {\boldsymbol{u}^{\prime 2}}{2} - \frac {5}{2} \, \frac {2}{3}k} \right ), \end{align}

then (2.32) can be simplified as

(2.35) \begin{equation} L\left [ {\frac {{{f^{\left ( 1 \right )}}}}{{{f_M}}}} \right ] = \frac {{{\tau _T}}}{{2k/3}}\left [ {{A_{ij}}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + {B_j}\frac {{\partial \ln k}}{{\partial {x_j}}}} \right ]. \end{equation}

Notably, it is easy to demonstrate that $A_{ij}$ and $B_j$ are eigenfunctions of operator $L [ {\ } ]$ , specifically,

(2.36) \begin{align} L\left [ {{A_{ij}}} \right ] &= - 2{A_{ij}},\nonumber \\[6pt] L\left [ {{B_j}} \right ] &= - 3{B_j}. \end{align}

So $f^{ ( 1 )}$ must be a linear combination of $A_{ij}$ and $B_j$ . We can assume the form of $f^{ ( 1 )}$ as

(2.37) \begin{equation} \frac {{{f^{\left ( 1 \right )}}}}{{{f_M}}} = \frac {{{\tau _T}}}{{2k/3}}\left ( {a \, {A_{ij}} + b \, {B_j}} \right ), \end{equation}

and the parameters can be determined easily, i.e.

(2.38) \begin{align} a &= - \frac {1}{2}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}, \nonumber \\[6pt] b &= - \frac {1}{3}\frac {{\partial \ln k}}{{\partial {x_j}}}. \end{align}

Consequently, we obtain the expression for $f^{ ( 1 )}$ as

(2.39) \begin{equation} {f^{\left ( 1 \right )}} = {f_M} \, \frac {{{\tau _T}}}{{2k/3}}\left [ { - \frac {1}{2}\left ( {{u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}} \right )\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac {1}{3}{u'_j}\left ( {\frac {{\boldsymbol{u}{'^2}}}{2} - \frac {5}{2} \, \frac {2}{3}k} \right )\frac {{\partial \ln k}}{{\partial {x_j}}}} \right ]. \end{equation}

The first-order approximation of Reynolds stress can be determined by taking moments of $f^{ ( 1 )}$ , i.e.

(2.40) \begin{equation} \sigma _{ij}^{T\left ( 1 \right )} = - \int _{{R^3}} {mu{'_{\left \langle i \right .}}u{'_{\left . j \right \rangle }}{f^{\left ( 1 \right )}}\,\mathrm{d}\boldsymbol{u}'} = \frac {{{\tau _T}}}{2} \, \rho \frac {2}{3}k\left [ {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac {{\partial {{\bar u}_j}}}{{\partial {x_i}}} - \frac {2}{3}\frac {{\partial {{\bar u}_k}}}{{\partial {x_k}}}{\delta _{ij}}} \right ]. \end{equation}

If $\tau _T$ is determined by $3{\nu _T}/k$ , then (2.40) together with the isotropic part $-2/3\rho k$ is exactly the linear EVM, i.e.

(2.41) \begin{equation} {\sigma _{ij}^T} = - \frac {2}{3}\rho k{\delta _{ij}} + 2\rho {\nu _T}S_{ij}^d, \end{equation}

where $S_{ij}^d$ is the strain tensor of the mean flow, i.e.

(2.42) \begin{equation} S_{ij}^d = \frac {1}{2}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac {1}{2}\frac {{\partial {{\bar u}_j}}}{{\partial {x_i}}} - \frac {1}{3}\frac {{\partial {{\bar u}_k}}}{{\partial {x_k}}}{\delta _{ij}}. \end{equation}

2.4. Second-order expansion: nonlinear EVMs

As noted in the previous subsection, the parameter $\xi$ can be regarded as a Knudsen number for turbulence, with the Boussinesq approximation’s applicability relying on it being small. The linear EVM offers a limited representation of anisotropy and is effective only for simple shear flows. To some extent, a mean turbulent flow more closely resembles a flow with a finite Knudsen number (Wilcox Reference Wilcox1998). For example, when the Reynolds number exceeds 20,000, an overlap region exists between the inner layer ( $y \lt 0.1\delta$ ) and the outer layer ( $y \gt 50{\delta _v}$ ) of the boundary layer (Pope Reference Pope2001), where $y$ denotes the distance from the wall. Here, $\delta$ represents the boundary layer thickness, while $\delta _v$ is the viscosity scale, defined as ${\delta _v} \equiv \nu \sqrt {\rho /{\tau _w}}$ , $\nu$ denotes the kinematic viscosity, $\rho$ is the density, and $\tau _w$ is the wall shear stress. In this region, the characteristic length scale of large-scale turbulent eddies $L$ , given by $L \equiv {k^{3/2}}/\varepsilon$ , and the mixing length $l_{mix}$ both vary linearly with $y$ (Wilcox Reference Wilcox1998; Pope Reference Pope2001). The proportionality constant for $L/y$ is approximately 2.5, while that for ${l_{mix}}/y$ corresponds to the Kármán constant, $\kappa \approx 0.41$ . Consequently, the turbulent Knudsen number is typically of order 1, i.e.

(2.43) \begin{equation} Kn = \frac {{{l_{mix}}}}{L} = \frac {{0.41y}}{{2.5y}} \approx 0.164. \end{equation}

This implies that the contributions from higher-order terms become non-negligible. To enhance the representation of Reynolds stress anisotropy, we next expand the turbulent Fokker–Planck equation to second order.

Starting with (2.19) and applying calculations similar to those in the previous subsection, we can reorganise (2.19) as

(2.44) \begin{equation} \frac {1}{{{f_M}}}\frac {\partial }{{\partial {u_i}}}\left [ {\frac {2}{3}k \, {f_M}\frac {\partial }{{\partial {u_i}}}\left ( {\frac {{{f^{\left ( 2 \right )}}}}{{{f_M}}}} \right )} \right ] = - \frac {1}{2}\frac {{\tau _T^2}}{{2k/3}}\left [ { - \left ( {{E_1} + {E_2}} \right ) + A + B + C + D} \right ], \end{equation}

where

(2.45) \begin{eqnarray} {E_1} &=& 2 \times \frac {2}{3}k \times {u'_i}\frac {{\partial S_{ij}^d}}{{\partial {x_j}}}, \nonumber\\[6pt] {E_2} &=& 2 \times \frac {2}{3}k \times \frac {1}{3}\left ( {\frac {{{u'_i}{u'_i}}}{{4k/3}} - \frac {3}{2}} \right ) \times 2S_{ij}^d\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}},\nonumber\\[6pt] A &=& \frac {1}{{2k/3}}{\left [ {\left ( {{u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}} \right )\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ]^2},\nonumber\\[6pt] B &=& \left [ {\frac {{{\partial _0}}}{{\partial t}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) + {{\bar u}_l}\frac {\partial }{{\partial {x_l}}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right )} \right ] \times \left ( {{u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}} \right )\nonumber\\[4pt] && {}+ {u'_l}\frac {\partial }{{\partial {x_l}}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) \times \left ( {{u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}} \right ),\nonumber\\[4pt] C &=& - \left [ {{u'_l}{u'_j}\frac {{\partial {{\bar u}_i}}}{{\partial {x_l}}} + {u'_l}{u'_i}\frac {{\partial {{\bar u}_j}}}{{\partial {x_l}}} - \frac {2}{3}{u'_k}{u'_l}\frac {{\partial {{\bar u}_k}}}{{\partial {x_l}}}{\delta _{ij}}} \right ]\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}},\nonumber\\[6pt] D &=& - \frac {1}{\rho }\frac {{\partial {\sigma _{lj}}}}{{\partial {x_j}}}\left [ {\frac {{{u_l}'}}{{2k/3}} \times \left ( {{u'_i}{u'_j} - \frac {1}{3}{u'_k}{u'_k}{\delta _{ij}}} \right )\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right. \nonumber\\[4pt] && \left. {+ \left ( {{u'_j}\frac {{\partial {{\bar u}_l}}}{{\partial {x_j}}} + {u'_i}\frac {{\partial {{\bar u}_i}}}{{\partial {x_l}}} - \frac {2}{3}{u'_l}{\delta _{ij}}\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right )} \right ]. \end{eqnarray}

We further express (2.44) in terms of the eigenfunctions, with the expressions of the eigenfunctions provided in Appendix B, as

(2.46) \begin{align} L\left [ {\frac {{{f^{\left ( 2 \right )}}}}{{{f_M}}}} \right ] &= - \frac {1}{2}\frac {{\tau _T^2}}{{2k/3}} \left \{ - \frac {4}{3}k \, {Q_{{E_1}}} - \frac {8}{9}k \, S_{ij}^d\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} \, {Q_{{E_2}}} + \frac {1}{{2k/3}} \, {Q_A} \right . \nonumber \\[4pt] &\quad\left . {}+ \left [ {\frac {{{\partial _0}}}{{\partial t}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) + {{\bar u}_l}\frac {\partial }{{\partial {x_l}}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right )} \right ] {Q_{{B_1}}} + {Q_{{B_2}}} - {Q_C} - \frac {1}{{2k/3}}{Q_{{D_1}}} - {Q_{{D_2}}} \right \}. \\[-24pt] \nonumber \end{align}

We can assume that $f^{ ( 2 )}$ is the combination of eigenfunctions of the Fokker–Planck operator, i.e.

(2.47) \begin{align} \frac {{{f^{\left ( 2 \right )}}}}{{{f_M}}} & = - \frac {1}{2}\frac {{\tau _T^2}}{{2k/3}} \left ( p {Q_{{E_1}}} + q {Q_{{E_2}}} + x {Q_A} + y {Q_{{B_1}}} + z {Q_{{B_2}}} \right . \nonumber \\[2pt] & \quad \left . {}+ u {Q_C} + v {Q_{{D_1}}} + w {Q_{{D_2}}} + m D + n M \right ), \end{align}

where $p,q,x,y,z,u,v,w,m,n$ are parameters that need to be determined. By substituting (2.47) into the Fokker–Planck operator and using the relationships provided in Appendix B, we find that

(2.48) \begin{equation} \begin{aligned} \!L\left [ {\frac {{{f^{\left ( 2 \right )}}}}{{{f_M}}}} \right ] = - \frac {1}{2}\frac {{\tau _T^2}}{{2k/3}}&\left [\!\vphantom{\left ( {\frac {4}{3}k \cdot u + \frac {8}{3}k \cdot n} \right )} - p \cdot {Q_{{E_1}}} - 2q \cdot {Q_{{E_2}}} - 4x \cdot {Q_A} - 2y \cdot {Q_{{B_1}}} - 3z \cdot {Q_{{B_2}}} - 2u \cdot {Q_C} \right . \\[6pt] & \,\ - 3v \cdot {Q_{{D_1}}} + \left ( {\frac {4}{3}k \cdot v - w} \right ){Q_{{D_2}}} + \left ( {\frac {4}{3}k \cdot z - m} \right ) \cdot D \\[6pt] & \,\ \left .+ \left ( {\frac {4}{3}k \cdot x - 2n} \right ) \cdot M + \left ( {\frac {4}{3}k \cdot u + \frac {8}{3}k \cdot n} \right ) \cdot N \right ]. \end{aligned} \end{equation}

By comparing (2.48) with (2.46), we obtain eleven equations for the ten unknown parameters. Evidently, these equations are consistent, and the parameters can be determined as

(2.49) \begin{equation} \begin{aligned} p &= \frac {4}{3}k,\quad q = \frac {4}{9}k \cdot S_{ij}^d\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}},\quad x = - \frac {1}{{8k/3}},\quad y = - \frac {1}{2}\left [ {\frac {{{\partial _0}}}{{\partial t}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) + {{\bar u}_l}\frac {\partial }{{\partial {x_l}}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right )} \right ],\\[4pt] z &= - \frac {1}{3},\quad u = \frac {1}{2},\quad v = \frac {1}{{2k}},\quad w = \frac {5}{3},\quad m = - \frac {4}{9}k,\quad n = - \frac {1}{4}. \end{aligned} \end{equation}

So far, we have derived the expression for $f^{ ( 2 )}$ , i.e.

(2.50) \begin{align} {f^{\left ( 2 \right )}} = \frac {1}{2}\frac {{\tau _T^2}}{{2k/3}} \cdot {f_M}\left ( - \frac {4}{3}k \cdot {Q_{{E_1}}} - \frac {4}{9}k \cdot S_{ij}^d\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} \cdot {Q_{{E_2}}} + \frac {1}{{8k/3}} \cdot {Q_A} \right . \nonumber \\[4pt] + \frac {1}{2}\left [ {\frac {{{\partial _0}}}{{\partial t}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right ) + {{\bar u}_l}\frac {\partial }{{\partial {x_l}}}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}}} \right )} \right ] \cdot {Q_{{B_1}}} + \frac {1}{3} \cdot {Q_{{B_2}}} \nonumber \\[4pt] \left . - \frac {1}{2} \cdot {Q_C} - \frac {1}{{2k}}{Q_{{D_1}}} - \frac {5}{3}{Q_{{D_2}}} + \frac {4}{9}k \cdot D + \frac {1}{4} \cdot M \right ). \end{align}

The second-order expansion of Reynolds stress can be obtained by taking moments of $f^{ ( 2 )}$ term by term, with the relevant integrals presented in Appendix C. After reorganisation, the result can be expressed as

(2.51) \begin{equation} \begin{aligned} \sigma _{ij}^{T\left ( 2 \right )} &= - \frac {1}{2}\rho \tau _T^2 \cdot \left ( {\frac {2}{3}k} \right ) \cdot \left ( {\frac {{{\partial _0}}}{{\partial t}} + {{\bar u}_k}\frac {\partial }{{\partial {x_k}}}} \right )S_{ij}^d\\[4pt] &\quad {}- \rho \tau _T^2 \cdot \left ( {\frac {2}{3}k} \right ) \cdot \left ( {S_{ik}^dS_{kj}^d - \frac {1}{3}S_{mn}^dS_{mn}^d{\delta _{ij}}} \right ) + \frac {1}{2}\rho \tau _T^2 \cdot \left ( {\frac {2}{3}k} \right ) \cdot ( {S_{ik}^d{\Omega _{kj}} + S_{jk}^d{\Omega _{ki}}} ). \end{aligned} \end{equation}

The turbulent relaxation time $\tau _T$ can be related to the eddy viscosity $\nu _T$ according to the relationship

(2.52) \begin{equation} {\nu _T} = \frac {{{\tau _T}}}{2} \cdot \frac {2}{3}k, \end{equation}

which is consistent with the first-order expansion. In this way, (2.51) can be expressed in a more formal manner, i.e.

(2.53) \begin{align} \sigma _{ij}^{T\left ( 2 \right )} &= - 2\rho \frac {{\nu _T^2}}{{2k/3}} \cdot \left ( {\frac {{{\partial _0}}}{{\partial t}} + {{\bar u}_k}\frac {\partial }{{\partial {x_k}}}} \right )S_{ij}^d \nonumber \\[4pt] &\quad{}- 4\rho \frac {{\nu _T^2}}{{2k/3}} \cdot \left ( {S_{ik}^dS_{kj}^d - \frac {1}{3}S_{mn}^dS_{mn}^d{\delta _{ij}}} \right ) + 2\rho \frac {{\nu _T^2}}{{2k/3}} \cdot ( {S_{ik}^d{\Omega _{kj}} + S_{jk}^d{\Omega _{ki}}} ). \end{align}

In (2.53), $\Omega _{ij}$ is the vorticity tensor of the mean velocity field, which is defined as

(2.54) \begin{equation} {\Omega _{ij}} \equiv \frac {1}{2}\left ( {\frac {{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac {{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right ). \end{equation}

This second-order expansion, as shown in (2.53), provides a more accurate representation of the Reynolds stress by including the effects of higher-order velocity gradients, which are essential for capturing complex turbulent behaviours such as secondary flow in a square duct.